Introduction

This lab is concerned with solution of ordinary differential equations (ODEs) using a Matlab function for the solution. You will see it applied first to a simple scalar equation, then to a system of equations, and then to a higher order equation converted into a system. In later labs, you will be writing your own ODE solver routines so you can understand the underlying theory.

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

Here is shorthand for the derivative . The equation is

An *analytic solution* of an ODE is a formula
,
that we can evaluate, differentiate, or analyze in any way we want.
However, analytic solutions can only be determined for a limited
class of ODE's.

A *numerical solution* of an ODE is simply a table of
abscissæ and approximate values
that approximate the value of an analytic solution. The values
are often called ``steps,'' especially for initial value
problems.
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
are to the analytic solution, which we might write as
.

This lab will take three sessions. If you print this lab, you may prefer to use the pdf version.

Matlab ODE solvers

Matlab has a number of built-in ODE solvers. These include:

Matlab ODE solvers | |

ode23 | non-stiff, low order |

ode113 | non-stiff, variable order |

ode15s | stiff, variable order, includes DAE |

ode23s | stiff, low order |

ode23t | trapezoid rule |

ode23tb | stiff, low order |

ode45 | non-stiff, medium order (Runge-Kutta) |

Do not worry about the meaning of ``stiff'' in this context, you will see it later this term. All of these functions use the very best methods, are highly reliable, use adaptive step size control, and allow close control of errors and other parameters. As a bonus, they can be ``told'' to precisely locate interesting events such as zero crossings and will even allow user-written functions to be called when certain types of events occur. There is very little reason to write your own ODE solvers unless you are actively researching new methods. In the following exercises you will see how to use these solvers.

In the following exercise, you will solve the following initial value problem (IVP).

**Exercise 1**:- Copy the following code to a function m-file named
`ex1_ode.m`.function yprime=ex1_ode(x,y) % yprime=ex1_ode(x,y) % computes the right side of the ODE y'=sin(x)-y % x,y are scalars, yprime is y' % your name and the date yprime=sin(x)-y;

- You can watch a solution evolve if you call the solver
without any output variables. Use
`ode45`to solve the IVP in Equation (1) on the interval`[0,15]`with the following Matlab command.ode45('ex1_ode',[0,15],1)

You should see that the solution becomes very much like after a short period of adjustment to the initial condition that is different from . The circles on the plot indicate points at which Matlab has chosen to evaluate solution, based on its default error criteria. Please include this plot with your summary. - If you wish to examine the solution in more detail, or
manipulate it, you need the solution values, not a picture. You
can get the solution values with commands similar to the following:
[x,y]=ode45('ex1_ode',[0,15],1);

If you wish, you can plot the solution, compute its error, etc. Plot this solution (plot`y`vs`x`) to see that it looks the same as before. Please send me this plot with your summary also. - For this solution, what is the value of
at
`x=15`? (Please use at least 6 decimal places for this answer. You can use the command`format long`to get more than the default amount of printed accuracy.) How many steps did it take? (The length of`x`is one more than the number of steps because it includes the initial value.) - Suppose you want a very accurate solution. The default
tolerance is .001 relative accuracy in Matlab, but suppose you
want a relative accuracy of 1.e-8? There is an extra variable
to provide options to the solver. It works in the following manner:
myoptions=odeset('RelTol',1.e-8) [x,y]=ode45('ex1_ode',[0,15],1,myoptions);

(The semicolon is left off the odeset call so you can see some of the other options and their defaults. Use`help odeset`for more detail.) How many steps did it take this time? What is the value of at`x=15`(to 6 decimal places)?

- Copy the following code to a function m-file named

ODE systems

The vast majority of ODEs that are solved are actually
*systems* of differential
equations, that is, a set of equations for which the unknown function
is a vector, say of dimension
. In this case, we might
write an individual component of the vector as
or
. Assuming that the
right hand side function
returns a vector, then the system of
equations has the same form as the scalar equation, and we can apply
the same ODE solvers to this system.

Here's an implementation issue: Matlab uses the variables `x`
and `y` as column vectors. The value `y(k)` is the
approximate solution at `x(k)`.
Since this as a column vector,
it's possible to handle the case of a system by assuming the
individual solution values are *row* vectors. In this
setting, a numerical solution will be a *table* of values, an
initial condition specifies the first row of that table, and each
row specifies the behavior of one of the components of the solution.

Goodwin's economic model

A model for cyclic economic behavior due to R. M. Goodwin ("A Growth Cycle", 1967, in Feinstein, editor, Socialism, Capitalism and Economic Growth) will be described below. This model results in a system of ordinary differential equations that we will be able to solve using Euler's method. The primary role of this model in the study of economics is that it predicts cyclical behavior without a cyclical driving force-it just arises naturally from the model itself. It is not necessary for you to understand the discussion of the foundation of the model or its economic interpretation. We will be using it as a source of a system of differential equations to be solved.

The model involves two primary dependent variables, and . The variable is the share of the total economy represented by employee salaries

where

The employment rate, is the second primary dependent variable:

where

The first assumption is that there is a constant, , representing the average return on capital for all investments. Thus, the gross national product is a fixed percentage of total capital invested.

The second assumption is that the money left over after wages have been paid is assumed to be invested, increasing the total capital available. Thus the total capital increases in the following way.

The following expression, in terms of relative rates of change, is a consequence of (2) and (3).

The third assumption is that the relative rates of change of productivity and population are both constants. Productivity is the average value of goods produced by employees

and

A simple consequence of (5) is that

The fourth assumption is that the relative change in wages is a function of employment.

A consequence of the derivative of definition of employment rate, , (7), (4) and (6) is one of the two differential equations for Goodwin's model

A consequence of the derivative of the definition of wage share, , is the other equation for Goodwin's model

In summary, Goodwin's model can be written

where the Greek letters signify positive constants.

**Exercise 2**:- Copy the following code to a Matlab function m-file
called
`goodwin_ode.m`, that returns the derivative function as a column vector:function yprime = goodwin_ode ( x, y ) % comments % your name and the date ALPHA =0.4; BETA =0.9; DELTA =0.05; NU =0.01; RHO =5.0; yprime=zeros(2,1); yprime(1)=-(DELTA+ALPHA)*???+BETA*y(1)*y(2); yprime(2)=(1/RHO-DELTA-NU)*???-y(1)*y(2)/RHO;

- Add comments just below the signature to indicate the calling and return sequences and to provide a summary of what the function does.
- Replace the symbols
`???`with the proper quantities so`goodwin_ode`implements Goodwin's model equations. - What is the role of the line
yprime=zeros(2,1);

What would happen if it were left out? - Note that the variable
`x`, corresponding to time, is unused. Nonetheless, it must be present in the signature. - Now use
`ode45`to integrate from to starting from y(1)=0.5, y(2)=0.5;yInit = [ 0.5 0.5 ]; ode45 ( 'goodwin_ode', [ 0.0, 100.0 ], yInit );

- Why do you get two curves? Please include a copy of this plot with your summary.
- Open a blank plot window with the command
`figure(2)`. Solve the same system, but this time use`ode15s`. Your solution should look the same, but you should see a slightly different distribution of circles indicating a different set of step sizes.

- Copy the following code to a Matlab function m-file
called

Higher Order ODE's

The previous examples are for first order equations and systems. What about higher order differential equations, ones that include derivatives such as ? There is a standard trick for transforming higher order equations into first order systems. You may have seen this trick in an ODE course. It also shows up in control theory in the conversion to ``state-space form'' of a system model and, I am sure, many other places. Atkinson describes the trick in Section 6.1, page 340. The trick will be illustrated by example here.

Consider a fourth-order differential equation

with initial conditions

The first step is to define new variables

(

**Exercise 3**: A pendulum swings through a circular arc. At any time, the angle the pendulum makes with the downward vertical is described by Newton's law ( ):

0

- Rewrite this as a first order system. (Hint: set
and
.) Write a function m-file named
`pendulum_ode.m`to represent this system. Be sure to put comments after the signature line and wherever else they are needed. Be sure that you return a*column*vector. - Use
`ode45`to solve for up to . Please include a plot of your solution with your summary file.

- Rewrite this as a first order system. (Hint: set
and
.) Write a function m-file named

**Exercise 4**: Your team has a cannon that shoots water balloons with a velocity of 100 ft/sec. You can adjust its angle of inclination, , to change the range. The other team is in front of a wall 10 ft high, 300 ft away from you. If you hit the wall with a water balloon, you will soak the other team. Your plan is to use Matlab to solve the equations of motion and trial and error to choose the correct value of . The trick is to pick a value of that results in the balloon being between 0 and 10 ft above ground at 300 ft.The equations of motion of a cannonball of unit mass under the influence of gravity are:

0

where denotes distance along the ground, and denotes height. The independent variable, , is time. The initial conditions are

0 0

where is a constant you will discover by trial and error.This is a combination of two second-order systems. Define the variables:

- Write a function m-file called
`cannon_ode.m`to implement the system. Be sure to include comments and to have the function return a*column*vector. - Use
`ode45`to integrate the equations of motion and determine, by trial and error, a value of so that the cannon ball hits the wall. You can determine success by checking that when . Watch out: the condition is on , not . You can find the first step beyond using the Matlab`find`statement. The first index,`k`for which is`k=min(find(y(:,1)>=300))`. Alternatively, you can plot the balloon's trajectory and see where it hits the ground or the wall. You can plot the trajectory withplot ( y(:,1), y(:,3) )

You might find the`axis([xmin,xmax,ymin,ymax])`command in Matlab useful to zoom in on the target, or you can use the magnifying icon on the plot window. (Click on the magnifying glass and then use the mouse to outline a square by clicking on the upper left and dragging. Warning: this might cause Matlab to crash on the computers in GSCC.)

`cannon_ode.m`and the value of`alpha`you used when you turn in this exercise.- Write a function m-file called

Phase Plane Plots

In many problems, a system is defined by a second order ODE whose right hand side does not depend on time. Such a system is called ``autonomous.'' The behavior of the solution doesn't depend on when you start the problem, only on the initial condition. Most of the above examples are autonomous. Especially for these kind of problems, it can be useful to look at the solution on a special ``phase plane plot.'' This simply means that instead of plotting versus the solution , we plot versus . This is an especially good way to detect periodic solutions, decaying transients, and so on. Phase plane plots are most interesting when the differential equation is not linear.

Phase plane plots are valuable because they highlight periodic or almost periodic behavior. Think about it: the solution to the pendulum case is , and that makes . If you plot this in the phase plane ( on the vertical axis, on the horizontal axis) you get a circle. Of course, once you go around the circle once, the autonomous nature of the system guarantees that you will continue going around forever, so the periodic behavior is highlighted. If you add a little damping to the pendulum equation, the phase plane plot becomes a slow spiral into the center.

**Exercise 5**: An equation for motion of a pendulum including a frictional term proportional to velocity is- Copy your file
`pendulum_ode.m`to`decay_ode.m`and modify it to include the frictional term . - Use
`ode45`to solve the resulting ODE for between 0 and 20, starting from the same initial condition as in Exercise 3 above,`[1;0]`. - Plot it in phase space (
along the horizontal axis
and
on the vertical axis) using the command
plot(y(:,1),y(:,2));

Please send me the plot.

- Copy your file

**Exercise 6**: In the previous exercise, all three coefficients on the left in the ODE are positive. If the coefficient of velocity is negative instead of positive, the physical result is that the solution grows instead of decaying.- Copy your file
`decay_ode.m`to`grow_ode.m`and modify it by changing the coefficient of the velocity term from`.2`to`-.2`. - Use
`ode45`to solve the resulting ODE for between 0 and 20, starting from`[1;0]`. - Plot
`y`vs.`x`to see one view of what the solution looks like. You do not have to send me this plot. - Plot it in phase space (
along the horizontal axis
and
on the vertical axis) using the command
plot(y(:,1),y(:,2));

You should observe an outward spiral. Please send me this phase-plane plot.

- Copy your file

One simple nonlinear equation with quite complicated behavior is van der Pol's equation. This equation is written

For positive , it can be seen by analogy with

**Exercise 7**: This exercise illustrates phase plane plots using the van der Pol ODE.- Rewrite the van der Pol equation as a first-order system, and
place the code in a function m-file named
`vanderpol_ode.m`. Use the value . Be sure you place comments where they belong. - Use
`ode45`to solve the system of ODEs for`x`between`x=0`to`x=20`, starting from an initial condition of`yInitial = [0; 0.25]`(notice the semicolon, making it a column vector). Make a phase plane plot from your solution. - Now, on the same plot (use the Matlab command
`hold on`), plot the solution starting from the initial condition`[0;3]`. Send me the plot and your*guesses*about:- What would the plot look like if you started at other points close to the ones you already did?
- What would the plot look like if you started very near the origin?
- What would the plot look like if you started exactly on the origin?

- Rewrite the van der Pol equation as a first-order system, and
place the code in a function m-file named

Stiff systems

ODEs can be either stiff or non-stiff. We will look at stiff ODEs again later in the term, but it is instructive to look at one practical method for distinguishing stiff from non-stiff ODEs. Basically, you try both non-stiff and stiff solvers and note which is faster.

Take my word that the IVP

can be regarded as a stiff system.

**Exercise 8**:- Implement the right side of the above IVP, Equation
(14) in a function m-file named
`stiff_ode.m`. You may want to start from`ex1_ode.m`. - Integrate it from
`x=0`to`x=15`, starting from`y=1`using`ode45`. Measure the time it takes with the commandstic;[x,y]=ode45('stiff_ode',[0,15],1);toc

How long does it take? How many steps did it take (length(x)-1)? - Use the following command to plot your result, using a blue line.
plot(x,y,'b')

- Repeat the above steps using
`ode15s`. Again, how much time does it take (should be much smaller)? How many steps does it take (length(x)-1)? Plot the results in red on the same frame as above. You should see the lines almost on top of one another, so the solutions are essentially the same. Please include this plot with your summary.

- Implement the right side of the above IVP, Equation
(14) in a function m-file named

Roundoff errors

Last term you saw some effects of roundoff errors. Later in this term you will look at roundoff errors again. Right now, though, is a good time to look at how some roundoff errors come about.

In the exercise below we will have occasion to use a special matrix called the Frank matrix. Row of the Frank matrix has the formula:

The Frank matrix for looks like:

The determinant of the Frank matrix is 1, but is difficult to compute numerically. This matrix has a special form called

**Exercise 9**: Let's look carefully at the Frank matrix and its inverse. For convenience, define`A`to be the Frank matrix of order 6, and`Ainv`its inverse, computed using`frank`and`frank_inv`. Similarly, let`B`and`Binv`be the Frank matrix of order 24 and its inverse.*Do not use the Matlab*You know that both`inv`function for this exercise!`A*Ainv`and`B*Binv`should equal the identity matrices of order 6 and 24 respectively. Let's just look at the top left entry. ComputeA(1,:)*Ainv(:,1)= __________ B(1,:)*Binv(:,1)= __________

Both of these answers should equal 1. The first does and the second does not. To see what goes right, compute the terms:A(1,6)*Ainv(6,1)= __________ A(1,5)*Ainv(5,1)= __________ A(1,4)*Ainv(4,1)= __________ A(1,3)*Ainv(3,1)= __________ A(1,2)*Ainv(2,1)= __________ A(1,1)*Ainv(1,1)= __________ sum = __________

Note that the signs alternate, so that when you add them up, each term tends to cancel part of the preceeding term. Now, to see what goes wrong, compute the terms:B(1,24)*Binv(24,1)= __________ B(1,23)*Binv(23,1)= __________ B(1,22)*Binv(22,1)= __________ B(1,21)*Binv(21,1)= __________ B(1,20)*Binv(20,1)= __________ B(1,16)*Binv(16,1)= __________ B(1,11)*Binv(11,1)= __________ B(1,6) *Binv(6,1) = __________ B(1,1) *Binv(1,1) = __________

You can see what happens to the sum. The first few terms are huge compared with the final sum of 1. Matlab uses 64-bit floating point numbers, so you can only rely on the first thirteen or fourteen significant digits in numbers like`B(1,24)*Binv(24,1)`. Further, they are of opposing signs so that there is extensive cancellation. There simply are not enough bits in the calculation to get anything like the correct answer. (Remark: It would not have been productive to compute each of the products`B(1,k)*Binv(k,1)`for each`k`, so I had you do the five largest and then sampled the rest. I chose to sample the terms with an*odd*-sized interval between adjacent terms. Had I chosen an even interval-say every other term-the alternating sign pattern would have been obscured.*When you are sampling errors or residuals for any reason,*)**never**take every other term!

Back to MATH2071 page.

Mike Sussman 2007-12-24