# Solving Ordinary Differential Equations11

Posted by Loren Shure,

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).

### Simple System

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:

### Conditions, Initial and Otherwise

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:

### Transform Equation

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 ,

### Try It!

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.

### Transforming Equations

Do you need to recast problems to fit into a particular formulation such as for solving ODEs? Let's hear about them here.

Get the MATLAB code

Published with MATLAB® 7.10

### Note

Jason Merrill replied on : 1 of 11

Matlab has recently grown some symbolic capabilities, right? How ’bout an interface that does this kind of mechanical manipulation automatically?

Aslak Grinsted replied on : 2 of 11

The flow of ice (i.e. as in glaciers) can be modelled as diffusion of ice-thickness. There is also an additional flux accounting for snowfall and melt. I’ve experimented with solving that using the ode solvers. I have a youtube video of that here:

http://www.glaciology.net/Home/Miscellaneous-Debris/ashallowiceglaciermodel

The main difficulty is the fact that the ice-thickness cannot be less than zero. I’ve not found a workable solution to that problem.

We have used the ode solvers in teaching a course on simple climate models. One of the models is the Stommel 1961 model of the Atlantic overturning circulation.

Loren replied on : 3 of 11

Jason-

Good idea, I’ll show that soon since it is possible to use the Symbolic Math Toolbox to manipulate the equations and generate the code for you.

–Loren

Loren replied on : 4 of 11

JuanPi-

That would depend on the format expected by the solvers so I don’t have a clear answer for you.

–Loren

JuanPi replied on : 5 of 11

Hi,
Good post!

When dealing with PDEs I usually use a finite method approach for solving the spatial part of the problem.
Time integration I usually do it by hand (explicit methods), because I haven’t really managed to formulate the problem as an implicit ODE, for example.
What is the usual/most-useful approach (apart from using the PDE toolbox)?

Thanks!

Loren replied on : 6 of 11

JuanPi-

I do not. I don’t know why you feel you can’t use a for loop instead. It might be more efficient since some ICs may require fewer steps than other ones. A joint solver will have them all take the most number of steps that the worst-case solution needs. Of course, if you use a for-loop, be sure to preallocate results.

–Loren

JuanPi replied on : 7 of 11

Yes, I guess it depends on the solver.

So another related question would be.

Let say I have the ODE of the mass-spring you showed. But I am trying to solve multiple initial conditions one one call to the integrator.

The solution I am used to use is to write a function of the form

function dxdt = multiSpring(t,x);
k=(2*pi)^2; m=1;
dxdt = zeros(size(x));
dxdt(1:2:end) = x(2:2:end);
dxdt(2:2:end) = -k * x(1:2:end)/ m;


and the solve it

[t y]= ode45(@multiSpring,tsan, [0 0, 0 -.1, 0 .1]);


where I have to send each initial value as a vector of initial conditions.

This works but makes everything a little bit more complicated to understand. Now every second column of the solution represent the same variable. So if I want to plot the value of x for the different initial conditions I have to do

plot(t,y(:,1:2:end))


Is there a way to use the columns of the differential equations defined in multiSpring?

That is, when I need to evaluate the differential equations for different values of the states of the system (useful for quiver plots), I write the system as

function dxdt = multiSpring2(t,x);
k=(2*pi)^2;m=1;
dxdt = zeros(size(x));
dxdt(1,:) = x(2,:);
dxdt(2,:) = -k * x(1,:)/ m;


(note that now the meaning of the rows and columns is clear).

Then

[X V]=meshgrid(linspace(-1,1,20));
J = multiSpring2(0,[X(:)';V(:)']);
JJ=reshape(J,2,20,20);
quiver(X,V,JJ(1,:,:),JJ(2,:,:)/10);


gives a sketch of the phase flow.

Do you know a trick to integrate multiple initial conditions similarly?

Thank you very much

Loren replied on : 8 of 11

JuanPi-

I don’t know which would be faster. You should time it. It may depend on how many points you compute.

–Loren

JuanPi replied on : 9 of 11

Loren-

mmm… That is a good point!

Nevertheless using a for loop (for the same reason you point out) forces you to store all the time vectors that the ode generates for each IC (since now easy solution will have less points). To do that, I have found no other solution but to use cells, something like,

[t{i} sol{i}]=ode45(...


where t and sol are preallocated cell arrays.
Another solution I tried is to ask for the solution at given times (by giving an array in the time span argument of the ode solver)… but in these case I could have all the solutions calculated at once! as I did before.

Noiw, what would be faster? A for-loop with preallocated cells (non-vectorized ode)or a big array of solutions at defined time values (vectorized ode)?

Thanks again!

Sandeep replied on : 10 of 11

Hi Loren,

Thanks a lot for that!! I’d had a lot of trouble with understanding odes until your easy to follow example!

Cheers,
Sandeep

Amit Kumar replied on : 11 of 11

Firstly , it looks your concepts are clear but some where you have not explained it in a well manner .But your explanation is quit similar to tutorvista.com .I have seen it .You must also visit that explanation see this link How to solve equations