I have recently handled several help requests for solving differential equations in MATLAB. All of the cases I worked on boil down to how to transform the higher-order equation(s) given to a system of first order equations. In this post I will outline how to accomplish this task and solve the equations in question. What I am not going to talk about is details of ODE solving algorithms (maybe another time).
Let's start with a simple, linear, second order system, a mass-spring system with motion constrained to one-dimension, a horizontal one (and therefore no gravity). If is the mass and is the spring constant, the equations of motion for the system are:
To solve this system, we need to know , , , and initial conditions, e.g., (also known as position and velocity). Let's simplify things and set , i.e., no external forces. Let's also set some initial conditions, , in other words, start with the spring unstretched and the mass moving. I end up with this system:
Looking in the help, I need to set up an system of equations to enable me to use one of the numerical ODE solvers in MATLAB.
To start the transformation, let me define a new variable that I will substitute in the system.
I can derive
and now rewrite my ODE system in terms of .
with initial conditions as above.
Now let me reorganize these 2 equations in a vector/matrix equation where
I now write me equation solely in terms of , the new vector (consisting of position and velocity). With that in mind, I will reorganize the existing equations first so I have on the left-hand sides.
or, in terms of ,
Let's try it. Set values for m and k
m = 1; k = 10;
I can create the ODE code in a file, or I can set up the equations as an anonymous function which is what I'll do here.
springmass = @(t,z)[z(2); -k * z(1)/ m];
Set up the initial conditions.
ic = [0; 1];
Solve between .
tspan = [0 10];
Call an ODE solver and plot the results.
[t,y] = ode23(springmass, tspan, ic); plot(t,y(:,1)), title('Position vs. Time')
And the velocity.
plot(t,y(:,2)), title('Velocity vs. Time')
You can see the left-most points on the plots match the initial conditions specified.
Do you need to recast problems to fit into a particular formulation such as for solving ODEs? Let's hear about them here.