|Matlab ODE solvers||Exercise 2|
|ODE systems||Exercise 3|
|Goodwin's economic model||Exercise 4|
|Higher Order ODE's||Exercise 5|
|Phase Plane Plots||Exercise 6|
|Stiff systems||Exercise 7|
|Roundoff errors||Exercise 8|
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:
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 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|
|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).
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;
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.
[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.
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)?
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.
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 wages refers to the average amount paid to each employee, numberEmployed refers to the number of employed people, and GNP refers to the gross national product, the total value of all production. Of course, the total of all wages paid to all workers is given by (wages)*(numberEmployed).
The employment rate, is the second primary dependent variable:
where population is the total number of people. The independent variable is time, .
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 population was used above in the definition of . Thus
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
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;
yprime=zeros(2,1);What would happen if it were left out?
yInit = [ 0.5 0.5 ]; ode45 ( 'goodwin_ode', [ 0.0, 100.0 ], yInit );
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
with initial conditions
The equations of motion of a cannonball of unit mass under the
influence of gravity are:
This is a combination
of two second-order systems. Define the variables:
plot ( 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.)
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.
plot(y(:,1),y(:,2));Please send me the plot.
plot(y(:,1),y(:,2));You should observe an outward spiral. Please send me this phase-plane plot.
One simple nonlinear equation with quite complicated behavior is van der Pol's equation. This equation is written
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
tic;[x,y]=ode45('stiff_ode',[0,15],1);tocHow long does it take? How many steps did it take (length(x)-1)?
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 Hessenberg form wherein all elements below the first subdiagonal are zero. Matlab provides the Frank matrix in its ``gallery'' of matrices, gallery('frank',n), but we will use an m-file frank.m. The inverse of the Frank matrix also consists of integer entries and an m-file for it can be downloaded as frank_inv.m. You can find more information about the Frank matrix from the Matrix Market, and the references therein.
A(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.