Download Ordinary Differential Equations - Lecture Notes | PHYS 321 and more Study notes Mechanics in PDF only on Docsity!
Chapter 7
Ordinary Differential
Equations
Matlab has several different functions for the numerical solution of ordinary dif- ferential equations. This chapter describes the simplest of these functions and then compares all of the functions for efficiency, accuracy, and special features. Stiffness is a subtle concept that plays an important role in these comparisons.
7.1 Integrating Differential Equations
The initial value problem for an ordinary differential equation involves finding a function y(t) that satisfies
dy(t) dt = f (t, y(t))
together with the initial condition
y(t 0 ) = y 0
A numerical solution to this problem generates a sequence of values for the indepen- dent variable, t 0 , t 1 ,... , and a corresponding sequence of values for the dependent variable, y 0 , y 1 ,... , so that each yn approximates the solution at tn
yn ≈ y(tn), n = 0, 1 ,...
Modern numerical methods automatically determine the step sizes
hn = tn+1 − tn
so that the estimated error in the numerical solution is controlled by a specified tolerance. The Fundamental Theorem of Calculus gives us an important connection be- tween differential equations and integrals.
y(t + h) = y(t) +
∫ (^) t+h
t
f (s, y(s))ds
2 Chapter 7. Ordinary Differential Equations
We cannot use numerical quadrature directly to approximate the integral because we do not know the function y(s) and so cannot evaluate the integrand. Nevertheless, the basic idea is to choose a sequence of values of h so that this formula allows us to generate our numerical solution. One special case to keep in mind is the situation where f (t, y) is a function of t alone. The numerical solution of such simple differential equations is then just a sequence of quadratures.
yn+1 = yn +
∫ (^) tn+
tn
f (s)ds
Through this chapter, we frequently use “dot” notation for derivatives,
y˙ =
dy(t) dt and ¨y =
d^2 y(t) dt^2
7.2 Systems of Equations
Many mathematical models involve more than one unknown function, and second and higher order derivatives. These models can be handled by making y(t) a vector- valued function of t. Each component is either one of the unknown functions or one of its derivatives. The Matlab vector notation is particularly convenient here. For example, the second-order differential equation describing a simple har- monic oscillator
x¨(t) = −x(t)
becomes two first-order equations. The vector y(t) has two components, x(t) and its first derivative ˙x(t),
y(t) =
[
x(t) x˙(t)
]
Using this vector, the differential equation is
y˙(t) =
[
x˙(t) −x(t)
]
[
y 2 (t) −y 1 (t)
]
The Matlab function defining the differential equation has t and y as input arguments and should return f (t, y) as a column vector. For the harmonic oscillator, the function is
function ydot = harmonic(t,y) ydot = [y(2); -y(1)]
A fancier version uses matrix multiplication in an inline function.
f = inline(’[0 1; -1 0]*y’,’t’,’y’);
4 Chapter 7. Ordinary Differential Equations
The most important term in this series is usually the one involving J, the Jacobian. For a system of differential equations with n components,
d dt
y 1 (t) y 2 (t) .. . yn(t)
f 1 (t, y 1 ,... , yn) f 2 (t, y 1 ,... , yn) .. . fn(t, y 1 ,... , yn)
the Jacobian is an n-by-n matrix of partial derivatives
J =
∂f 1 ∂y 1
∂f 1 ∂y 2...^
∂f 1 ∂yn ∂f 2 ∂y 1
∂f 2 ∂y 2...^
∂f 2 ∂yn .. .
∂fn ∂y 1
∂fn ∂y 2...^
∂fn ∂yn
The influence of the Jacobian on the local behavior is determined by the solution to the linear system of ordinary differential equations
y˙ = Jy
Let λk = μk + iνk be the eigenvalues of J and Λ = diag(λk) the diagonal eigenvalue matrix. If there is a linearly independent set of corresponding eigenvectors V , then
J = V ΛV −^1
The linear transformation
V x = y
transforms the local system of equations into a set of decoupled equations for the individual components of x,
x˙k = λkxk
The solutions are
xk(t) = eλk^ (t−tc)x(tc)
A single component xk(t) grows with t if μk is positive, decays if μk is negative, and oscillates if νk is nonzero. The components of the local solution y(t) are linear combinations of these behaviors. For example the harmonic oscillator
y˙ =
[
]
y
is a linear system. The Jacobian is simply the matrix
J =
[
]
7.4. Single-Step Methods 5
The eigenvalues of J are ±i and the solutions are purely oscillatory linear combi- nations of eit^ and e−it. A nonlinear example is the two-body problem
y˙(t) =
y 3 (t) y 4 (t) −y 1 (t)/r(t)^3 −y 2 (t)/r(t)^3
where
r(t) =
y 1 (t)^2 + y 2 (t)^2
In an exercise, we ask you to show that the Jacobian for this system is
J =
r^5
0 0 r^5 0 0 0 r^5 2 y 12 − y 22 3 y 1 y 2 0 0 3 y 1 y 2 2 y^22 − y^21 0
It turns out that the eigenvalues of J just depend on the radius r(t)
λ =
r^3 /^2
i −
−i
We see that one eigenvalue is real and positive, so the corresponding component of the solution is growing. One eigenvalue is real and negative, corresponding to a decaying component. Two eigenvalues are purely imaginary, corresponding to os- cillatory components. However, the overall global behavior of this nonlinear system is quite complicated, and is not described by this local linearized analysis.
7.4 Single-Step Methods
The simplest numerical method for the solution of initial value problems is Euler’s method. It uses a fixed step size h and generates the approximate solution by
yn+1 = yn + hf (tn, yn) tn+1 = tn + h
The Matlab code would use an initial point t0, a final point tfinal, an initial value y0, a step size h, and an inline function or function handle f. The primary loop would simply be
t = t0; y = y0; while t <= tfinal y = y + h*feval(f,t,y) t = t + h end
7.4. Single-Step Methods 7
s 2 = f (tn +
h 2 , yn +
h 2 s 1 )
s 3 = f (tn +
h 2
, yn +
h 2
s 2 ) s 4 = f (tn + h, yn + hs 3 )
yn+1 = yn + h 6
(s 1 + 2s 2 + 2s 3 + s 4 ) tn+1 = tn + h
If f (t, y) does not depend on y, then classical Runge-Kutta has s 2 = s 3 and the method reduces to Simpson’s quadrature rule. Classical Runge-Kutta does not provide an error estimate. The method is sometimes used with a step size h and again with step size h/2 to obtain an error estimate, but we now know more efficient methods. Several of the ODE solvers in Matlab, including the textbook solver we describe later in this chapter, are single-step or Runge-Kutta solvers. A general single-step method is characterized by a number of parameters, αi, βi,j , γi, and δi. There are k stages. Each stage computes a slope, si, by evaluating f (t, y) for a particular value of t and a value of y obtained by taking linear combinations of the previous slopes.
si = f (tn + αih, yn + h
∑^ i−^1
j=
βi,j sj ), i = 1,... , k
The proposed step is also a linear combination of the slopes.
yn+1 = yn + h
∑^ k
i=
γisi
An estimate of the error that would occur with this step is provided by yet another linear combination of the slopes.
en+1 = h
∑^ k
i=
δisi
If this error is less than the specified tolerance, then the step is successful and yn+ is accepted. If not, the step is a failure and yn+1 is rejected. In either case, the error estimate is used to compute the step size h for the next step. The parameters in these methods are determined by matching terms in Taylor series expansions of the slopes. These series involve powers of h and products of various partial derivatives of f (t, y). The order of a method is the exponent of the smallest power of h that cannot be matched. It turns out that one, two, three, or four stages yield methods of order one, two, three or four, respectively. But it takes six stages to obtain a fifth-order method. The classical Runge-Kutta method has four stages and is fourth-order. The names of the Matlab ODE solvers are all of the form odennxx with digits nn indicating the order of the underlying method and a possibly empty xx indicating
8 Chapter 7. Ordinary Differential Equations
some special characteristic of the method. If the error estimate is obtained by comparing formulas with different orders, the digits nn indicate these orders. For example, ode45 obtains its error estimate by comparing a fourth-order and a fifth- order formula.
7.5 The BS23 Algorithm
Our textbook function ode23tx is a simplified version of the function ode23 that is included with Matlab. The algorithm is due to Bogacki and Shampine [3, 6]. The “ 23 ” in the function names indicates that two simultaneous single step formulas, one of second order and one of third order, are involved. The method has three stages, but there are four slopes si because, after the first step, the s 1 for one step is the s 4 from the previous step. The essentials are
s 1 = f (tn, yn)
s 2 = f (tn + h 2
, yn + h 2
s 1 )
s 3 = f (tn +
h, yn +
hs 2 ) tn+1 = tn + h
yn+1 = yn +
h 9
(2s 1 + 3s 2 + 4s 3 ) s 4 = f (tn+1, yn+1)
en+1 =
h 72 (− 5 s 1 + 6s 2 + 8s 3 − 9 s 4 )
The simplified pictures in figure 7.1 show the starting situation and the three stages. We start at a point (tn, yn) with an initial slope s 1 = f (tn, yn) and an estimate of a good step size, h. Our goal is to compute an approximate solution yn+1 at tn+1 = tn + h that agrees with true solution y(tn+1) to within the specified tolerances. The first stage uses the initial slope s 1 to take an Euler step halfway across the interval. The function is evaluated there to get the second slope, s 2. This slope is used to take an Euler step three quarters of the way across the interval. The function is evaluated again to get the third slope, s 3. A weighted average of the three slopes
s =
(2s 1 + 3s 2 + 4s 3 )
is used for the final step all the way across the interval to get a tentative value for yn+1. The function is evaluated once more to get s 4. The error estimate then uses all four slopes.
en+1 = h 72
(− 5 s 1 + 6s 2 + 8s 3 − 9 s 4 )
If the error is within the specified tolerance, then the step is successful, the tentative value of yn+1 is accepted, and s 4 becomes the s 1 of the next step. If the error is
10 Chapter 7. Ordinary Differential Equations
A fourth input argument is optional, and can take two different forms. The simplest, and most common, form is a scalar numerical value, rtol, to be used as the relative error tolerance. The default value for rtol is 10−^3 , but you can provide a different value if you want more or less accuracy. The more complicated possibility for this optional argument is the structure generated by the Matlab function odeset. This function takes pairs of arguments that specify many different options for the Matlab ODE solvers. For ode23tx, you can change the default values of three quantities, the relative error tolerance, the absolute error tolerance, and the M-file that is called after each successful step. The statement
opts = odeset(’reltol’,1.e-5, ’abstol’,1.e-8, ... ’outputfcn’,@myodeplot)
creates a structure that specifies the relative error tolerance to be 10−^5 , the absolute error tolerance to be 10−^8 , and the output function to be myodeplot. The output produced by ode23tx can be either graphic or numeric. With no output arguments, the statement
ode23tx(F,tspan,y0);
produces a dynamic plot of all the components of the solution. With two output arguments, the statement
[tout,yout] = ode23tx(F,tspan,y0);
generates a table of values of the solution.
7.6 ode23tx
Let’s examine the code for ode23tx. Here is the preamble.
function [tout,yout] = ode23tx(F,tspan,y0,arg4,varargin) %ODE23TX Solve non-stiff differential equations. % Textbook version of ODE23. % % ODE23TX(F,TSPAN,Y0) with TSPAN = [T0 TFINAL] integrates % the system of differential equations y’ = f(t,y) from % t = T0 to t = TFINAL. The initial condition is % y(T0) = Y0. The input argument F is the name of an % M-file, or an inline function, or simply a character % string, defining f(t,y). This function must have two % input arguments, t and y, and must return a column vector % of the derivatives, y’. % % With two output arguments, [T,Y] = ODE23TX(...) returns % a column vector T and an array Y where Y(:,k) is the % solution at T(k). %
7.6. ode23tx 11
% With no output arguments, ODE23TX plots the solution. % % ODE23TX(F,TSPAN,Y0,RTOL) uses the relative error % tolerance RTOL instead of the default 1.e-3. % % ODE23TX(F,TSPAN,Y0,OPTS) where OPTS = ODESET(’reltol’, ... % RTOL,’abstol’,ATOL,’outputfcn’,@PLTFN) uses relative error % RTOL instead of 1.e-3, absolute error ATOL instead of % 1.e-6, and calls PLTFN instead of ODEPLOT after each step. % % More than four input arguments, ODE23TX(F,TSPAN,Y0,RTOL, % P1,P2,..), are passed on to F, F(T,Y,P1,P2,..). % % ODE23TX uses the Runge-Kutta (2,3) method of Bogacki and % Shampine. % % Example % tspan = [0 2pi]; % y0 = [1 0]’; % F = ’[0 1; -1 0]y’; % ode23tx(F,tspan,y0); % % See also ODE23.
Here is the code that parses the arguments and initializes the internal variables.
rtol = 1.e-3; atol = 1.e-6; plotfun = @odeplot; if nargin >= 4 & isnumeric(arg4) rtol = arg4; elseif nargin >= 4 & isstruct(arg4) if ~isempty(arg4.RelTol), rtol = arg4.RelTol; end if ~isempty(arg4.AbsTol), atol = arg4.AbsTol; end if ~isempty(arg4.OutputFcn), plotfun = arg4.OutputFcn; end end t0 = tspan(1); tfinal = tspan(2); tdir = sign(tfinal - t0); plotit = (nargout == 0); threshold = atol / rtol; hmax = abs(0.1*(tfinal-t0)); t = t0; y = y0(:);
% Make F callable by feval.
7.6. ode23tx 13
Here is the error estimate. The norm of the error vector is scaled by the ratio of the absolute tolerance to the relative tolerance. The use of the smallest floating-point number, realmin, prevents err from being exactly zero.
e = h(-5s1 + 6s2 + 8s3 - 9*s4)/72; err = norm(e./max(max(abs(y),abs(ynew)),threshold), ... inf) + realmin;
Here is the test to see if the step is successful. If it is, the result is plotted or appended to the output vector. If it is not, the result is simply forgotten.
if err <= rtol t = tnew; y = ynew; if plotit if feval(plotfun,t,y,’’); break end else tout(end+1,1) = t; yout(end+1,:) = y.’; end s1 = s4; % Reuse final function value to start new step. end
The error estimate is used to compute a new step size. The ratio rtol/err is greater than one if the current step is successful, or less than one if the current step fails. A cube root is involved because the BS23 is a third-order method. This means that changing tolerances by a factor of eight will change the typical step size, and hence the total number of steps, by a factor of two. The factors 0.8 and 5 prevent excessive changes in step size.
% Compute a new step size. h = hmin(5,0.8(rtol/err)^(1/3));
Here is the only place where a singularity would be detected.
if abs(h) <= hmin warning(sprintf( ... ’Step size %e too small at t = %e.\n’,h,t)); t = tfinal; end end
That ends the main loop. The plot function might need to finish its work.
if plotit feval(plotfun,[],[],’done’); end
14 Chapter 7. Ordinary Differential Equations
7.7 Examples
Please sit down in front of a computer running Matlab. Make sure ode23tx is in your current directory or on your Matlab path. Start your session by entering
ode23tx(’0’,[0 10],1)
This should produce a plot of the solution of the initial value problem
dy dt
y(0) = 1 0 ≤ t ≤ 10
The solution, of course, is a constant function, y(t) = 1. Now you can press the up arrow key, use the left arrow key to space over to the ’0’, and change it to something more interesting. Here are some examples. At first, we’ll change just the ’0’ and leave the [0 10] and 1 alone.
F Exact solution ’0’ 1 ’t’ 1+t^2/ ’y’ exp(t) ’-y’ exp(-t) ’1/(1-3t)’ 1-log(1-3t)/3 (Singular) ’2y-y^2’ 2/(1+exp(-2t))
Make up some of your own examples. Change the initial condition. Change the accuracy by including 1.e-6 as the fourth argument. Now let’s try the harmonic oscillator, a second-order differential equation writ- ten as a pair of two first-order equations. First, create an inline function to specify the equations. Use either
F = inline(’[y(2); -y(1)]’,’t’,’y’)
or
F = inline(’[0 1; -1 0]*y’,’t’,’y’)
Then, the statement
ode23tx(F,[0 2*pi],[1; 0])
plots two functions of t that you should recognize. If you want to produce a phase plane plot, you have two choices. One possibility is to capture the output and plot it after the computation is complete.
[t,y] = ode23tx(F,[0 2*pi],[1; 0]) plot(y(:,1),y(:,2),’-o’) axis([-1.2 1.2 -1.2 1.2]) axis square
16 Chapter 7. Ordinary Differential Equations
y0 = [1; 0; 0; change_this]; [t,y] = ode23tx(@twobody,[0 2*pi],y0); plot(y(:,1),y(:,2),’-’,0,0,’ro’) axis equal
You might also want to use something other than 2π for tfinal.
7.8 Lorenz Attractor
One of the world’s most extensively studied ordinary differential equations is the Lorenz chaotic attractor. It was first described in 1963 by Edward Lorenz, an MIT mathematician and meteorologist, who was interested in fluid flow models of the earth’s atmosphere. An excellent reference is a book by Colin Sparrow [8]. We have chosen to express the Lorenz equations in a somewhat unusual way, involving a matrix-vector product.
y˙ = Ay
The vector y has three components that are functions of t,
y(t) =
y 1 (t) y 2 (t) y 3 (t)
Despite the way we have written it, this is not a linear system of differential equa- tions. Seven of the nine elements in the 3-by-3 matrix A are constant, but the other two depend on y 2 (t).
A =
−β 0 y 2 0 −σ σ −y 2 ρ − 1
The first component of the solution, y 1 (t), is related to the convection in the atmo- spheric flow, while the other two components are related to horizontal and vertical temperature variation. The parameter σ is the Prandtl number, ρ is the normal- ized Rayleigh number, and β depends on the geometry of the domain. The most popular values of the parameters, σ = 10, ρ = 28, and β = 8/3, are outside the ranges associated with the earth’s atmosphere. The deceptively simple nonlinearity introduced by the presence of y 2 in the system matrix A changes everything. There are no random aspects to these equa- tions, so the solutions y(t) are completely determined by the parameters and the initial conditions, but their behavior is very difficult to predict. For some values of the parameters, the orbit of y(t) in three-dimensional space is known as a strange attractor. It is bounded, but not periodic and not convergent. It never intersects itself. It ranges chaotically back and forth around two different points, or attractors. For other values of the parameters, the solution might converge to a fixed point, diverge to infinity, or oscillate periodically.
7.8. Lorenz Attractor 17
0 5 10 15 20 25 30
y
y
y
t
Figure 7.2. Three components of Lorenz attractor
−25−25 −20 −15 −10 −5 0 5 10 15 20 25
−
−
−
−
0
5
10
15
20
25
y
y
Figure 7.3. Phase plane plot of Lorenz attractor
Let’s think of η = y 2 as a free parameter, restrict ρ to be greater than one, and study the matrix
A =
−β 0 η 0 −σ σ −η ρ − 1
It turns out that A is singular if and only if
η = ±
β(ρ − 1)
7.9. Stiffness 19
function ydot = lorenzeqn(t,y,A) A(1,3) = y(2); A(3,1) = -y(2); ydot = A*y;
Most of the complexity of lorenzgui is contained in the plotting subfunction, lorenzplot. It not only manages the user interface controls, it must also anticipate the possible range of the solution in order to provide appropriate axis scaling.
7.9 Stiffness
Stiffness is a subtle, difficult, and important concept in the numerical solution of ordinary differential equations. It depends on the differential equation, the initial conditions, and the numerical method. Dictionary definitions of the word “stiff” involve terms like “not easily bent,” “rigid,” and “stubborn.” We are concerned with a computational version of these properties.
A problem is stiff if the solution being sought is varying slowly, but there are nearby solutions that vary rapidly, so the numerical method must take small steps to obtain satisfactory results.
Stiffness is an efficiency issue. If we weren’t concerned with how much time a computation takes, we wouldn’t be concerned about stiffness. Nonstiff methods can solve stiff problems; they just take a long time to do it. A model of flame propagation provides an example. We learned about this example from Larry Shampine, one of the authors of the Matlab ODE suite. If you light a match, the ball of flame grows rapidly until it reaches a critical size. Then it remains at that size because the amount of oxygen being consumed by the combustion in the interior of the ball balances the amount available through the surface. The simple model is
y˙ = y^2 − y^3 y(0) = δ 0 ≤ t ≤ 2 /δ
The scalar variable y(t) represents the radius of the ball. The y^2 and y^3 terms come from the surface area and the volume. The critical parameter is the initial radius, δ, which is “small.” We seek the solution over a length of time that is inversely proportional to δ. At this point, we suggest that you start up Matlab and actually run our examples. It is worthwhile to see them in action. We will start with ode45, the workhorse of the Matlab ODE suite. If δ is not very small, the problem is not very stiff. Try δ = .01 and request a relative error of 10−^4.
delta = 0.01; F = inline(’y^2 - y^3’,’t’,’y’); opts = odeset(’RelTol’,1.e-4); ode45(F,[0 2/delta],delta,opts);
20 Chapter 7. Ordinary Differential Equations
With no output arguments, ode45 automatically plots the solution as it is computed. You should get a plot of a solution that starts at y = .01, grows at a modestly increasing rate until t approaches 100, which is 1/δ, then grows rapidly until it reaches a value close to 1, where it remains. Now let’s see stiffness in action. Decrease δ by a couple of orders of magnitude. (If you run only one example, run this one.)
delta = 0.0001; ode45(F,[0 2/delta],delta,opts);
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 x 10^4
0
0.
0.
0.
0.
1
ode
0.98 1 1.02 1.04 1.06 1.08 1.1 1. x 10^4
0.
1
1.
Figure 7.4. Stiff behavior of ode
You should see something like figure 7.4, although it will take a long time to complete the plot. If you get tired of watching the agonizing progress, click the stop button in the lower left corner of the window. Turn on zoom, and use the mouse to explore the solution near where it first approaches steady state. You should see something like the detail in the figure 7.4. Notice that ode45 is doing its job. It’s keeping the solution within 10−^4 of its nearly constant steady state value. But it certainly has to work hard to do it. If you want an even more dramatic demonstration of stiffness, decrease the tolerance to 10−^5 or 10−^6. This problem is not stiff initially. It only becomes stiff as the solution ap- proaches steady state. This is because the steady state solution is so “rigid.” Any solution near y(t) = 1 increases or decreases rapidly toward that solution. (We
