MATH2071: LAB #3: Simple ODE Methods

Introduction
Euler's Method
Euler's Method Exercise
ODE Systems
Higher Order ODE's
ASSIGNMENT

Introduction

A very simple ordinary differential equation (ODE) is the explicit scalar first-order initial value problem:

y'(x) = f(x,y),
y(x₀) = y₀.

Here y' is shorthand for the derivative dy/dx. The equation is explicit because y' can be written explicitly as a function of x and y. It is scalar because we assume that y(x) is a scalar quantity. It is first-order because the highest derivative that appears is the first derivative y'. It is an initial value problem because we are given the value of the solution at some time or location x₀ and are asked to produce a formula for the solution at later times.

An analytic solution of an ODE is a formula y(x), which we can evaluate, differentiate, or analyze in any way we want. However, analytic solutions can only be determined for a very small class of ODE's.

A numerical solution of an ODE is simply a table of abscissas and approximate values (x_i,y_i) which are believed to approximate the value of an analytic solution. In general, a numerical solution is always wrong; the important question is, how wrong is it? One way to pose this question is to determine how close the computed values (x_i,y_i) are to the analytic solution, which we might write as (x_i,y(x_i)).

Euler's Method

The simplest method for producing a numerical solution of an ODE is known as Euler's method. Given a solution value (x_k,y_k), we estimate the next solution by:

x_k+1 = x_k + dx,
y_k+1 = y_k + dx * y'(x_k,y_k).

We can take as many steps as we want with this method, using the approximate answer from one step as the starting point for the next step. Each step may be of a different size, although we aren't ready to say why that can be useful.

Let's assume for the moment that we have the solution at x=0 and that we want the solution at x=1. We could take a single giant Euler step and get an approximate answer at x=1. We now have a numerical solution to the problem. But here are some important questions:

Is there any way to determine if this approximation is accurate?
If we're not happy with this approximation, what can we do?
If we take two steps of size 0.5, we get a different answer at x=1. Is the new answer better? Worse? Would 4 steps of size 0.25 be even better?

You may know, or have guessed, that we're going to want to take many small steps, and that usually, the smaller the steps, the more accurate the answer will be.

Euler's Method Exercise

Problem: Consider the initial value problem:

y'(x) = y
y(0) = 1.

M FILE: Write a MATLAB M file called euler.m, which has the form

function [ x, y ] = euler ( f, x_range, y_initial, nstep )

x(1) = x_range(1);
dx = ( x_range(2) - x_range(1) ) / nstep;
y(1) = y_initial;

for i = 1 : nstep
x(i+1) = x(i) + dx;
y(i+1) = y(i) + dx * feval ( f, x(i), y(i) );
end

The input to this function is:

f, the name of the derivative function, which will have to be in quotes;
x_range, a row vector of two values, the starting and stopping values of the independent variable;
y_initial, the value of the dependent variable at the initial time;
nstep, the number of Euler steps to take.

The function returns two vectors:

x, the value of the independent variable at nstep+1 equally spaced points;
y, the value of the dependent variable at each x.

Notice the use of the MATLAB routine feval. Because we're going to pass the name of the derivative function as an argument, we're not allowed to write:

y(i+1) = y(i) + dx * f(x(i),y(i));

but instead have to use the helper routine feval in this way.

M FILE: We will need a very simple M file to evaluate the derivative function. Write a MATLAB M file called model_ode.m, which has the form

function yprime = model_ode ( x, y )
yprime = y;

INTERACTIVE: Now we can use our implementation of Euler's method to "march" from x=0:

y_init = 1.0;
[ x, y ] = euler ( 'model_ode', [ 0.0, 2.0 ], y_init, 20 );
y_true = exp ( x );

You can use the plot(x,y) command to examine your approximate solution. You should also compare it to the true solution with the command plot(x,y_true).

COMPUTATION: Compute the numerical solution of the model ODE, from x = 0.0 to x = 2.0, using Euler's method and stepsizes of 0.2, 0.1, 0.05 and 0.025. (You have to translate the stepsize into the appropriate number of steps nstep). For each case, record the value of the numerical solution at x = 2.0, and the error, that is, the difference between the numerical solution and the true solution at the end point.

        Stepsize DX    Computed Y(2.0)      Error

        0.2            __________           __________
        0.1            __________           __________
        0.05           __________           __________
        0.025          __________           __________

When your stepsize is divided by 2, what happens to the error?

ODE Systems

Euler's method is easily extended to a system of differential equations, that is, a set of equations for which the unknown function y(x) is a vector, say of dimension n. In this case, we might write an individual component of the vector as y_i or y_i(x). Assuming that the right hand side function f returns a vector, then the system of equations has the same form as the scalar equation, and we can apply Euler's method componentwise. The only complication (which doesn't affect Euler's method itself) is that each component derivative y'_i can depend on any or all of the components y_iof the solution function.

Here's an implementation issue: we've already been using the variable y as a vector, namely, the list of its values at the corresponding x values. Since we've done this as a row vector, it's possible for us to handle the new case by assuming the individual solution values are column vectors. In this setting, a numerical solution will be a table of values, an initial condition specifies the first column of that table, and each row specifies the behavior of one of the components of the solution. If we think about all this in the right way, MATLAB will make it easy to extend Euler's method.

M FILE: Modify euler.m so that it can handle ODE systems. (It will still solve scalar problems too).

function [ x, y ] = euler ( f, x_range, y_initial, nstep )

x(1) = x_range(1);
dx = ( x_range(2) - x_range(1) ) / nstep;
y(:,1) = y_initial;

for i = 1 : nstep
x(i+1) = x(i) + dx;
y(:,i+1) = y(:,i) + dx * feval ( f, x(i), y(:,i) );
end

Notice that we don't tell MATLAB how the dimension of the system. It figures that out. The interesting thing here is the use of the "colon" index. We use two indices to access parts of the numerical solution table. Using a colon for the first index, and specifying a value for the second (which is the column) means we are going to examine or modify all the data in a particular column of the table.

A simple and interesting example of an ODE system is a predator-prey model. The variables are:

y'₁, the rate of change of a prey population, subject to a natural growth rate, from which is subtracted a loss factor depending on the number of predators;
y'₂, the rate of change of the predator population, subject to a natural death rate, to which is added a positive growth factor depending on the availability of prey.

A simple form of such a model might be:

y'₁ = alpha - beta * y₂
y'₂ = - gamma + delta * y₁

For a predator-prey model, assume the parameter values are

alpha = 1.0
beta = 0.5
gamma = 0.75
delta = 0.25
x(0) = 0
y₁(0) = 5.0
y₂(0) = 1.0

M FILE: Write a MATLAB M file called predator_ode.m, which returns the derivative function as a column vector:

function yprime = predator_ode ( x, y )
yprime = [ 1.00 - 0.50 * y(2);
- 0.75 + 0.25 * y(1) ];

(The first semicolon is different from the second one. It separates rows of the vector!)

INTERACTIVE: Now use your implementation of Euler's method to integrate to x = 50. Start by using 10 steps. Double the number of steps, and compare solutions. Keep doubling until the solution seems "reasonable". If either variable becomes negative, you are definitely not using enough steps!

y_init = [ 5.0; 1.0 ];
[ x, y ] = euler ( 'predator_ode', [ 0.0, 50.0 ], y_init, 10 );

Repeat after me: MATLAB makes vectors easy! It's worth a little trouble to use them! Don't believe me? Try the command plot(x,y) now!

Higher Order ODE's

Euler's method can also be extended to higher order differential equation, that is an explicit equation involving a derivative like y''. In order to represent derivatives of arbitrary order, we may use the notation yⁿ or y^(n-1) to represent a derivative of order n or n-1, say.

For clarity, we'll suppose we have a second order initial value problem:

y'' = f(x,y,y')
y(x₀) = y₀
y'(x₀) = y₁

We make this look like a system by making up names for the solution y and its derivative:

z₁ = y
z₂ = y'

We can now write an equivalent set of first order equations:

z'₁ = z₂
z'₂ = f(x,z₁,z₂)

with initial conditions:

z₁ = y₀
z₂ = y₁

But now this is a problem we know how to solve! It's just like the predator-prey example. And this same idea of converting a higher order system will work for systems of order 3 and so on.

Exercise - a pendulum swings through a circular arc. At any time, the angle THETA(T) the pendulum makes with the downward vertical is described by:

THETA'' + 1.5 * THETA = 0

with initial conditions

T₀ = 0
THETA(T₀) = 1
THETA'(T₀) = 0

Rewrite this as a first order system, write the appropriate M file pendulum_ode.m, and solve for THETA up to T = 25. Why do we expect the absolute value of THETA never to exceed 1?

Assignment

Your lab instructor is standing 1000 feet away from you. He is about to expel you from the class. You decide to change his mind by firing an exploding cannonball, which will kill him if it is no more than 10 feet above him when it explodes.

The equations of motion of a cannonball under the influence of gravity are:

x'' = 0
y'' = - 32

with initial conditions

t₀ = 0
x(0) = 0
x'(0) = 100
y(0) = 0
y'(0) = ALPHA

where ALPHA is unspecified. This is a combination higher-order vector system. Define the variables:

z₁ = x
z₂ = x'
z₃ = y
z₄ = y'

Determine the form of the first order ODE system. Write an M file called cannon_ode.m. Since the cannonball explodes at time t = 10, we want to know the value of the solution at that time.

By trial and error, determine a value of ALPHA so that the cannon ball kills your lab instructor. You may find it very useful to plot the path (x,y) of your cannonball; use the command plot ( z(1,:), z(3,:) ).

Once you are done, you will need to report:

the value of ALPHA that you used;
the value of NSTEP, the number of Euler steps you took;
the position of the cannon ball at t=10, which had better be between (1000,0) and (1000,10).

Mail your results to me.

Last revised on 20 January 2000.