Ordinary Differential Equation (ODE) solvers

Introduction

The ODE solvers in MATLAB solve initial value problems with a variety of properties. In this activity we will learn how to use the most common ODE solver to solve first and second order ODEs.

Before starting

Use the MATLAB Live Editor to edit and run this Live Script in your browser or your desktop.

Read each section carefully.
Run the code within a section before starting to read the next.
To run the code from each section, click into the code section and then click on the Run Section button (from the toolstrip) or click on the blue stripe on left side of that section as shown below:

Remark: Run the code of each section from top to bottom, otherwise you may get an error.

- The end of a section is indicated with a thin line, like the next one -

1. Basic implementation

In order to solve an Initial Value Problem (IVP) in MATLAB we can use the following MATLAB code:

[t,y] = ode45(odefun, tspan, y0)

where

[t,y]: Defines two column vectors to store values of t and . Each row in the solution array y corresponds to a value at the time given by the same row in t.
odefun: Ordinary differential equation function
tspan: Interval of integration
y0: Initial conditions
ode45: MATLAB function to solve ODEs

Remark: In MATLAB there are several ODE solvers that we could use (e.g. ode23, ode113, ode23s, ode23t, and ode23tb). However, in the following activities we will use only the solver ode45. You can learn more details about these commands here Choose an ODE Solver, which provides general guidelines on when to use each of the different solvers.

For example, consider the IVP problem

where

We can define

as an anonymous function:

f = @(t,y) 1 - exp(-4*t) - 2*y;

Then we solve the IVP with the command ode45:

[t,y] = ode45(f, [0 2], 1);

Finally we plot the solution:

plot(t, y,'r-', t, y, 'b+')

How would you change these commands to solve the IVP for the same ODE but with the condition

on the interval

? Try it and re-run this section to see the result.

2. Plotting solutions of IVP problems and slope fields

Consider now the following IVP

where

We can use the previous method to plot a numerically approximated solution to the IVP over the top of a slope field.

First, we need defined the ODE and use the command ode45:

g = @(t,y) 2*t*exp(-1/t) + y/t^2;
[gt, gy] = ode45(g, [1 2], 0.1);

Then we plot the slope field including the approximated solution. Run this section to see the result:

% Interval bounds
tmin = 1;
tmax = 2;
ymin = 0;
ymax = 2;
% Grid spacing
spacing = 0.15;
 
% Define x and y coordinates
[pt, py] = meshgrid(tmin:spacing:tmax,ymin:spacing:ymax);
 
% Calculate slopes
ydash = 2*pt.*exp(-1./pt) + py./pt.^2;
 
% Define slope vector components
dt = (spacing/2)./sqrt(1+ydash.^2);
dy = (spacing/2)*ydash./sqrt(1+ydash.^2);
 
% Clear the figure
clf;
 
% Create slope field
q = quiver(pt, py, dt, dy,'b','AutoScale','off');
set(q,'ShowArrowHead','off','LineWidth',1)
 
% Include approximated solution
hold on
 
plot(gt, gy,'r-', gt, gy, 'k+')
 
hold off
 
% Set axis and labels
axis([tmin tmax ymin ymax])
xlabel('x')
ylabel('y')

Note: In the previous code we used the command set() to set graphics object properties. For more details see: Set Graphics Object.

3. Systems of differential equations

Anonymous function can also be used to define systems of ODEs. The easiest way to think of this is that the function will take a vector of inputs Y and puts the inputs into a vector of equations.

For example, a second order ODE of the form

can be expressed as the coupled system of equations. Let

. Then we have that

. If we substitute in equation (1) we obtain the system:

Define the vector Y as

The derivative of Y is a vector whose entries correspond to the coupled equations:

Now we have two first order equations that we can store as a vector of anonymous functions with vector input Y. We do this by writting:

secODE = @(t,Y) [Y(2); -Y(2) + Y(1)+ cos(t)];

3.1 Numerical solution of systems of ODEs

Now consider the IVP

which can be re-written as

To find and plot a numerical solution of this system over the interval

, first we define the initial and final time; and the initial conditions at 0:

t0 = 0;   % Initial time
tf = 10;  % Final time, here we have chosen 10
z0 = 1;   % Initial condition z(0)=1
y0 = 2;   % Initial condition z'(0)=2

The anonymous function representing the system is then defined as:

secODE = @(t,Y) [Y(2); -Y(2) + Y(1) + cos(t)];

Now, we solve and plot the solution:

[T,Y] = ode45(secODE, [t0 tf], [z0; y0]);
plot(T,Y(:,1),'ro-');

Run this section to see the result. Change the initial conditions of the IVP, or modify the equation of the system. Re-run this section to see the new output.

Remark: Since we solved a coupled system of 2 equations there are 2 columns in Y. The first column Y(:,1) gives the values of

, while the second column gives the values of

Finally, if we wish to compute the approximate value for z at

, we just need to use the command

Y(end,1)

Run this section to see the output. The result should be approximately 879.

4. Hands on Practice

Let's practice what we just learned.

Activity 1

Consider the following IVP

Plot a numerically approximated solution to the IVP over the top of a slope field.
Estimate the value of .

Write your code here:

Activity 2

The figure below shows a pendulum with length L and the angle θ for the vertical to the pendulum.

The motion of this pendulum is determined by the θ as a function of time, which satisfies the nonlinear differential equation

where g is the acceleration due to gravity.

In your notebook, re-write this second order ODE as a system of differential equations. Then compute an approximated solution to the IVP

Write your code here: