equations and state the types of problems we will solve using Matlab and. Maple in this chapter. ..... 8. Figure 4.2: Plot of a solution to a system of ODEs. 152 ...... For our model, we add a source block that represents the signal M(t). By editing its ...
Chapter 4
Solving Differential Equations 4.1
Introduction
In this chapter, we explore the use of both Matlab and Maple for solving what are, besides linear systems, among the most fundamental problems in applied mathematics: differential equations. We will see that in many cases, differential equations may be solved analytically, i.e. a formula for the solution can be obtained, in which case Maple is an ideal tool. More often, though, an analytical solution is not possible, in which case we must solve the equation numerically, and in this scenario both Matlab and Maple offer a suite of tools. Before we begin, we review some basic concepts regarding differential equations and state the types of problems we will solve using Matlab and Maple in this chapter.
4.1.1
Ordinary Differential Equations
An ordinary differential equation, or ODE, is an equation relating its solution y(t), which is unknown function of one variable, to its derivatives or other functions. The quantity y(t) is known as the dependent variable, as it depends on the quantity t, known as the independent variable. In most cases, the independent variable represents time, or spatial position in one dimension. The most common ODE is a first-order equation, which can be written 143
in the form
dy = f (t, y(t)) (4.1) dt where f is some function of both the independent and dependent variable. In general, the order of a differential equation is equal to the order of the highest derivative appearing in the equation. For example, the equation d2 y =y dt2
(4.2)
is a second-order equation. The solution to 4.1, if it exists, is not unique, but rather depends on an arbitrary constant. Example The ODE
dy = λy, (4.3) dt where λ is a constant, is y(t) = Ceλt , where C is an arbitrary constant. A unique solution may be obtained by prescribing an initial condition on the solution y(t), which typically takes the form y(t0 ) = y0
(4.4)
The problem described by the ODE (4.1) and initial condition (4.4) is known as the initial value problem. Example The initial value problem dy = λy, dt
y(0) = 1
(4.5)
has the unique solution y(t) = eλt . For a higher-order ODE, additional initial conditions must be specified to ensure the existence of a unique solution. In general, the number of initial conditions is equal to the order of the ODE. As we will see, however, ODE’s of higher order may be rewritten as a system of first-order ODE’s, where the solution is a vector-valued function of the independent variable.
4.1.2
Partial Differential Equations
A partial differential equation, or PDE, is an equation that relates its solution u, a function of several variables, to its partial derivatives with respect to these variables. u is the dependent variable, a function of several independent variables, which typically represent both spatial position and time. 144
4.2
Solving Differential Equations in Matlab
Matlab provides several methods for solving the general first-order differential equation y 0 (t) = f (t, y(t)), y(0) = y0 , (4.6) where f is a given function and y0 is a given vector of initial data. This differential equation is known as the initial value problem. Partial differential equations that involve only first derivatives in time can be solved numerically by solving an initial value problem.
4.2.1
Runge-Kutta methods
Runge-Kutta methods are used to compute y(t n+1 ) from y(tn ) as accurately as possible by evaluating f (t, y) at select points in the interval [t n , tn+1 ]. The more points that are used, the greater the accuracy, but also the greater the computational effort required. Matlab provides functions for Runge-Kutta methods of varying order of accuracy. Runge-Kutta Method of order 2 Let h be our time step tn+1 − tn , and let yn = y(tn ) = y(hn). We try to construct an approximation of the form yn+1 = yn + ak1 + bk2 ,
(4.7)
where k1 = hf (tn , yn ) k2 = hf (tn + αh, yn + βk1 )
(4.8) (4.9)
where a, b, α, and β are constants to be determined so that our approximation will be as accurate as possible. Using a Taylor series for y(t) centered at t = t n , we obtain, using the original differential equation, h3 h2 00 y (tn ) + y 000 (tn ) + · · · 2 6 h2 = y(tn ) + hf (tn , yn ) + (ft + f fy ) + O(h3 ) 2
yn+1 = yn + hy 0 (tn ) +
145
(4.10)
Furthermore, again using a Taylor series centered at t = t n , we have f (tn + αh, yn + βk1 ) = f (tn , yn ) + αhft + βk1 fy + O(h2 ).
(4.11)
Substituting the previous expression for f (t n + αh, yn + βk1 ) into the Taylor series for yn+1 yields yn+1 = yn + (a + b)hf + bh2 (αft + βf fy ) + O(h3 ), so by setting
1 a=b= , 2
α = β = 1,
(4.12) (4.13)
we obtain h f (tn , yn ) + 2 f (tn + h, yn + hf (tn , yn ))] + O(h3 )
yn+1 = yn +
(4.14)
Runge-Kutta method of order 4 Using a similar process, we can derive the 4th-order method, whose error is O(h5 ): 1 (4.15) yn+1 = yn + (k1 + 2k2 + 2k3 + k4 ), 6 where k1 = hf (tn , yn ) 1 h k2 = hf tn + , yn + k1 2 2 h 1 k3 = hf tn + , yn + k2 2 2 k4 = hf (tn + h, yn + k3 )
(4.16)
Matlab provides the function ode23 to use a combination of secondand third-order Runge-Kutta methods to solve y 0 (t) = f (t, y). Given the function f (t, y), an initial time t0 , initial data y0 = y(t0 ), and a final time tf , ode23 will compute the solution y(t) at various times between t 0 and tf and return a column vector of the values of t at which the solution is computed, as well as an array containing the values of y(t) at those times. In the following example, we wish to solve a scalar ODE, y 0 (t) = f (t, y(t)) = −y(t),
y(0) = 1
(4.17)
which has the exact solution y(t) = exp(−t). First, we must create an M-file, in this case called an ODE file, that defines our function f (t, y). 146
elaine12:~> cat > simpleode.m function f = odefile(t,y) f = -y; Now we can call ode23. Our arguments are the name of our ODE file, in single quotes and without the .m extension, a row vector containing our initial and final t-values, and our value of y 0 . >> [T,Y]=ode23(’simpleode’,[0 10],1); Since we know the exact solution in this case, we will plot it against our approximate solution. >> >> >> >> >> >> >>
Yexact=exp(-T); plot(T,Y,’ro’,T,Yexact,’g-’) xlabel(’t’) ylabel(’y’) legend(’Computed solution’,’Exact solution’) title(’Solutions of y’’=-y, y(0)=1’) norm(Yexact-Y)
ans = 0.0014 As the plot in Figure 4.1 shows, our approximation is quite accurate. Alternatively, we can specify specific times at which we wish y(t) to be computed. >> [T,Y]=ode23(’simpleode’,0:2:10,1) T = 0 2 4 6 8 10
147
Solutions of y’=−y, y(0)=1 1 Computed solution Exact solution
0.9 0.8 0.7
y
0.6 0.5 0.4 0.3 0.2 0.1 0
0
1
2
3
4
5 t
6
7
8
9
10
Figure 4.1: Plot of solution to 4.17 (red circles) vs. exact solution (green solid curve)
148
Y = 1.0000 0.1348 0.0182 0.0024 0.0003 0.0000 Options to ODE solvers The odeset function can be used to further customize the solution process. It accepts pairs of option names and option values, where option names are strings, in single quotes. It returns a structure containing the new options, which may be passed back to odeset if options in an existing structure are to be modified. The odeget function retrieves option values from an existing option structure. It accepts the option structure and an option name as arguments. Some commonly used options are: • RelTol is used to specify the relative error tolerance for each component yi of the solution y, i = 1, . . . , n where n = length(y). For each step of Runge-Kutta that is performed, the error e i in yi satisfies ei = max(RelTol|yi |, AbsTol(i)).
(4.18)
The default value is 0.001. • AbsTol is used to specify the absolute error. Set AbsTol to a scalar to impose an error bound that applies to all components of the solution, or a vector to impose different bounds on different components. The default value is 0.001. • Refine increases the number of output points to produce smoother output. The number of output points is multipled by Refine. Its default value is 1. • OutputFcn refers to a function that should be called after every time step. If an ODE solver is called with no output arguments, the default function is odeplot. Otherwise, no output function is called, by default. Output functions are called with three arguments at the start of integration: the tspan and y0 arguments passed to the solver, and 149
the string ’init’. After Then, it is called with two arguments, t and y, representing the time and solution value at that time, after each time step. • OutputSel can be set to a vector of indices, indicating which components of the solution should be passed to OutputFcn. By default, all components are passed. • Stats can either be on or off (default). If on, computational cost statistics are displayed. • Jacobian is set to on (default is off) if the ODE file is written so that F(t,y,’jacobian’) returns ∂F/∂y. • JConstant is used to indicate that the Jacobian matrix ∂F/∂y is constant. Its default value is off; if the Jacobian is constant, it should be set to on. • JPattern should be set to on (default is off) if the ODE file is coded so that F([],[],’jpattern’) returns a sparse matrix with 1’s corresponding to nonzero entries in ∂F/∂y. • Vectorized can be set to on (default is off) if the ODE file is coded so that a set of y-values may be passed to F, i.e. F(t,[y1 y2 ...]) returns [F(t,y1) F(t,y2) ... ]. • MaxStep is used to impose an upper bound on the step size. • InitialStep is used to suggest an initial step size. • NormControl can be turned on (default is off) to control the error in each time step by the relationship kek ≤ max(RelTolkyk, AbsTol),
(4.19)
rather than bounding the error component-wise. >> options=odeset(’RelTol’,1.e-6,’Refine’,2,’Stats’,’on’,... ’JConstant’,’on’) options = AbsTol: [] BDF: [] Events: [] 150
InitialStep: Jacobian: JConstant: JPattern: Mass: MassConstant: MassSingular: MaxOrder: MaxStep: NormControl: OutputFcn: OutputSel: Refine: RelTol: Stats: Vectorized:
[] [] ’on’ [] [] [] [] [] [] [] [] [] 2 1.0000e-06 ’on’ []
In the following example, we solve a system of ODE’s. A plot of the solution is given in Figure 4.2. function F = odefile(t,y) A=diag(-1:2); F=A*y; >> ode23(’odesystem’,[ 0 1 ],ones(4,1),options) 70 successful steps 0 failed attempts 211 function evaluations 0 partial derivatives 0 LU decompositions 0 solutions of linear systems The ode45 function can be used to obtain greater accuracy with fewer time steps, employing a combination of fourth- and fifth-order methods. The trade-off is that more computation effort is expended per time step. Its usage is identical to that of ode23.
4.2.2
Multistep Methods
The Runge-Kutta methods are known as one-step methods because they rely only the solution at one previous time to obtain the solution at the next time. 151
8 7 6 5 4 3 2 1 0
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
Figure 4.2: Plot of a solution to a system of ODEs
152
1
Multistep methods compute the solution at a given time using values of the solution at multiple times in the past. The Adams-Bashforth methods are a class of such methods, of varying orders of accuracy. The following method is 5th-order, i.e. the error per time step h is O(h 5 ). yn+1 = yn +
h (55fn − 59fn−1 + 37fn−2 − 9fn−3 ). 24
(4.20)
The Adams-Moulton methods are predictor-corrector methods, which use a multistep method to predict the value of the solution at a future time based on current knowledge, and then correct using another multistep scheme that takes the information from the predictor into account. The ode113 solver uses various Adams-Bashforth and Adams-Moulton methods of orders 1 through 12. Its usage is identical to that of ode23 and ode45.
4.2.3
Higher-order equations
While so far we have only discussed first-order ODE’s, these solvers can just as easily be applied to higher-order equations. A higher-order ODE can be written as a system of first-order ODE’s by introducing new variables to represent the higher-order derivatives. The number of equations in the system is equal to the order of the highest derivative. For example, for the second-order equation dy d2 y − 2 + y = 0, dt dt
(4.21)
we introduce two dependent variables: y1 = y,
y2 = y 0 .
(4.22)
Then we can rewrite this equation as y10 = y2 y20 = 2y2 − y1
(4.23)
Some ODE’s are more difficult to solve because their solutions contain components with widely varying time scales. Such equations are called stiff. A simple example is the equation dy d2 y + 1001 + 1000y = 0, 2 dt dt 153
(4.24)
which has the general solution y(t) = Ae−t + Be−1000t .
(4.25)
The second component, Be−1000t , decays much more rapidly than the first component, Ae−t . As a result, an extremely small time step must be used in order to obtain an accurate solution, or even any solution at all. Systems of ODEs derived from PDEs tend to be stiff, so such equations are extremely common. Matlab provides a family of ODE solvers for use with stiff equations. All of these are used in the same way as other ODE solvers. • ode15s uses backward-difference formulas of orders 1-5 to solve stiff equations. • ode23s uses second- and third-order methods. • ode23tb uses implicit Runge-Kutta formulas of second- and thirdorder. Partial differential equations are often rephrased as systems of ODE’s for numerical solution. To make this transition, we view a function f (x, t) as a vector-valued function f (t), where the elements of f (t) are the values of f (x, t) at selected x-values called gridpoints. The function f (t) is called a gridfunction. Using gridfunctions, An initial-boundary value problem or IBVP of the form ∂ ∂u = P x, , t, u , u(x, 0) = f (x), (4.26) ∂t ∂x where u(x, t) must satisfy given boundary conditions, can be viewed as an initial value problem du = A(t)u(t), dt
u(0) = f ,
(4.27)
where A is a matrix that approximates the differential operator P and also enforces the boundary conditions. We now use this technique to solve the first-order wave equation with periodic boundary conditions: ∂u ∂u =a , ∂t ∂x
0 < x < 2π,
u(x, 0) = f (x), 154
t > 0,
0 < x < 2π,
(4.28) (4.29)
u(0, t) = u(2π, t),
t > 0.
(4.30)
We represent functions on the interval [0, 2π] using N equally spaced gridpoints, with spacing h = 2π/N . We then approximate the differentiation operator as follows: f (x + h) − f (x − h) ∂f = . ∂x 2h Since we are working with periodic functions, this circulant matrix A: 0 1 0 ··· 0 −1 0 1 0 ··· .. . a 0 −1 0 1 A= .. . . . . 2h .. .. .. .. . 0 · · · 0 −1 0 1 0 · · · 0 −1
(4.31) approximation yields a −1 0 .. .
. 0 1 0
(4.32)
We would like to be able to create this matrix only once and then use it for every time step. To that end, we pass A as an additional argument the function defined in to our ODE file. Any options to the function odefile beyond the third argument are assumed to be such additional parameters passed to the solver. We construct our ODE fils as follows: elaine21:~> cat > waveeqn.m function F = odefile(t,y,flag,A) F=A*y; The following Matlab code then sets up and solves the PDE, and plots the solution: >> >> >> >> >> >> >> >> >> >>
a=1; N=256; h=2*pi/(N+1); x=h*(0:N)’; f=sin(x); A=(a/(2*h))*(diag(ones(N,1),1)-diag(ones(N,1),-1)); A(N+1,1)=1/(2*h); A(1,N+1)=-1/(2*h); [T,Y]=ode23s(’waveeqn’,[0 1],f,options,A); plot(x,Y’) 155
Solution to du/dt = du/dx, u(x,0) = sin(x) 1 0.8 0.6 0.4
u(x,t)
0.2 0 −0.2 −0.4 −0.6 −0.8 −1
0
1
2
3
4
5
6
x
Figure 4.3: Solution to the first-order wave equation >> xlabel(’x’) >> ylabel(’u(x,t)’) >> title(’Solution to du/dt = du/dx, u(x,0) = sin(x)’) Each row of the solution matrix Y represents the solution u(x, t) for a given time t. The solution is plotted in Figure 4.3. Visualizing Solutions in Higher Dimensions Matlab provides two functions, odephas2 and odephas3, that allow visualization of solutions in two and three dimensions, respectively. To use one of these functions, the OutputFcn option should be set to the name of the function you wish to use. >> options=odeset(options,’OutputFcn’,’odephas2’); 156
7
0.35 0.3 0.25 0.2 0.15 0.1 0.05 0 −0.05 −0.2
0
0.2
0.4
0.6
0.8
1
1.2
Figure 4.4: Orbit of the solution of a system of two ODEs The above code will allow us to solve a system of two ODE’s and view the orbit of the solution [x(t), y(t)] in the plane. We define our ODE file as follows: elaine21:~> cat > orbit.m function f = odefile(t,y) f = [ -exp(-t)*sin(t)-exp(-t)*cos(t); ... -exp(-t)*sin(t)+exp(-t)*cos(t) ]; and then solve: >> ode23(’orbit’,[0 10],[1 0],options); The orbit of the solution is plotted in Figure 4.4 157
4.3
Solving ODEs in Maple
Maple offers a package called DEtools that provides several functions for analyzing and solving ordinary differential equations (ODEs). To use the functions in this package, it is recommended to use the with statement in your Maple session. > with(DEtools):
4.3.1
Specifying ODEs in Maple
An ODE is described to Maple in the same way as any other equation: in an expression of the form expr = expr, where either expression includes the unknown solution and some of its derivatives. Derivatives are specified using Maple’s diff function, or its D notation. The following two equations are equivalent, and illustrate the use of the diff and D forms. > ODE1 := diff(y(x),x$2) - diff(y(x),x) + y(x) = sin(x); / 2 \ |d | /d \ ODE1 := |--- y(x)| - |-- y(x)| + y(x) = sin(x) | 2 | \dx / \dx / > ODE2 := (D@@2)(y)(x) - D(y)(x) + y(x) = sin(x); (2) ODE2 := (D )(y)(x) - D(y)(x) + y(x) = sin(x) The dsolve function uses a variety of methods to attempt to solve an ODE analytically. If it is unable to construct a closed-form or explicit solution, it will instead try to return an implicit solution, defined by an equation. Alternatively, dsolve can be asked to solve an ODE numerically, as in Matlab. > ode1 := diff(y(x),x)-y(x)^2+y(x)*sin(x)-cos(x); /d \ 2 ode1 := |-- y(x)| - y(x) + y(x) sin(x) - cos(x) \dx /
158
> ans1 := dsolve(ode1); _C1 exp(-cos(x)) ans1 := y(x) = sin(x) - --------------------------/ | _C1 | exp(-cos(x)) dx + 1 | /
In the previous example, the solution contained an arbitrary constant C1. This is because the solution of the given ODE is only unique up to that constant. Alternatively, initial conditions can be supplied to dsolve, similar to those given to Matlab’s ODE solvers, in order to obtain a unique solution. Typically, initial conditions are specified by passing a set as the first argument to dsolve, consisting of all ODEs and any initial conditions. Sets in Maple are delimited by curly braces. We use the same ODE as in the previous example, but with an initial condition. > dsolve({ode1,y(0)=1}); exp(-cos(x)) y(x) = sin(x) - ----------------------------------------x / | -cosh(1) + sinh(1) + | exp(-cos(u)) du | / 0
> dsolve({diff(x(t),t)=-x(t),x(0)=2}); x(t) = 2 exp(-t) dsolve can use a number of methods to solve ODEs numerically as well as analytically. To use numerical methods, pass the argument type=numeric to dsolve. 159
Optionally, one can specify a numerical method to be used, using the method argument, which can have one of these values: • rkf45: Fehlberg fourth-fifth order Runge-Kutta method • dverk78: seventh-eighth order continuous Runge-Kutta method • classical: by default, the forward Euler method, which is simple but has poor convergence properties. One can specify method=classical[choice] where choice can be: – foreuler is forward Euler – heunform is the trapezoidal rule – impoly is the modified Euler method – rk2, rk3, and rk4 are second, third, and fourth-order Runge Kutta methods – adambash is the Adams-Bashforth method – abmoulton is the Adams-Moulton method • gear is the Gear single-step extrapolation method • mgear is the Gear multi-step method • lsode uses the Livermore Stiff ODE solver • taylorseries uses a Taylor series method We will solve this simple ODE numerically, using dsolve. > ode := diff(x(t),t)=-x(t); d ode := -- x(t) = -x(t) dt > ic := x(0)=2; ic := x(0) = 2 By default, dsolve returns a procedure, whose arguments are values of the independent variables of the ODE, and whose output is the value of the solution at those values. 160
> g := dsolve({ode,ic},x(t),type=numeric, method=classical[rkf45]); g := proc(rkf45_x) ... end proc > g(0); [t = 0, x(t) = 2.] > g(1); [t = 1, x(t) = .735758877735078864] Alternatively, one can use the value argument to specify a list of values at which the solution should be evalutaed. In this case, Maple returns a 2 × 1 matrix identifying all relevant variables and values for the solution. > dsolve({ode,ic},x(t),type=numeric, value=array([0,0.1,0.2])); [ [t, x(t)] ] [ ] [[0 2. ]] [[ ]] [[.1 1.809674836]] [[ ]] [[.2 1.637461506]] dsolve can also return a solution in series form. > ode := diff(y(t),t,t)+diff(y(t),t)^2=0; / 2 \ |d | /d \2 ode := |--- y(t)| + |-- y(t)| = 0 | 2 | \dt / \dt / When the initial conditions are not given, the answer is expressed in terms of the solution and its derivatives evaluated at the origin. > ans := dsolve( {ode}, y(t), type=series); 2 ans := y(t) = y(0) + D(y)(0) t - 1/2 D(y)(0)
161
2 t
3 + 1/3 D(y)(0)
3 t
-
4 4 5 5 6 1/4 D(y)(0) t + 1/5 D(y)(0) t + O(t ) If initial conditions are given, the series is centered at that the initial point. > ans := dsolve({ode,y(a)=Y_a,D(y)(a)=DY_a}, y(t),type=series); 2 ans := y(t) = Y_a + DY_a (t - a) - 1/2 DY_a
2 (t - a)
3 + 1/3 DY_a
4 4 5 5 6 1/4 DY_a (t - a) + 1/5 DY_a (t - a) + O((t - a) ) By default, the first six terms of the series are returned. The global variable Order can be used to obtain a different number of terms. > Order := 3; Order := 3 The odetest function can be used to verify that a solution returned by dsolve does in fact solve the ODE. > ODE := diff(y(x),x)=F((y(x)-x*ln(x))/x) + ln(x); d y(x) - x ln(x) ODE := -- y(x) = F(--------------) + ln(x) dx x > sol := dsolve(ODE,implicit); y(x) ---- - ln(x) x / | 1 sol := ln(x) + | -------------- d_a - _C1 = 0 | _a + 1 - F(_a) / > odetest(sol,ODE); 0 162
3 (t - a)
-
The odeadvisor function can be used to obtain information about an ODE that may be helpful in solving it or simply understanding it. Internally, dsolve uses the analysis techniques of odeadvisor to determine which method it should use to solve an ODE. > ODE := x*diff(y(x),x)+a*y(x)^2-y(x)+b*x^2; /d \ 2 2 ODE := x |-- y(x)| + a y(x) - y(x) + b x \dx / > odeadvisor(ODE); [[_homogeneous, class D], _rational, _Riccati]
4.3.2
Other solution techniques
Using the information obtained by odeadvisor, one can use a number of other Maple functions to solve an ODE. These methods are used by dsolve, but are still desirable if one wishes to have more control over which solution techniques are used. Exact Equations An ODE is exact if it can be written in the form u(x)
dv du + v(x) = 0, dx dx
(4.33)
since then we can apply the product rule for differentiation in reverse and obtain d(uv) = 0. (4.34) dx Maple’s exactsol determines whether an ODE is exact, and if so, solves it. > ode := 3*t^3*z(t)^2*diff(z(t),t)+3*t^2*z(t)^3 = 0: > exactsol( ode, z(t) ); 1/2 _C1 - 1/2 _C1 - 1/2 I 3 _C1 {z(t) = ---, z(t) = --------------------------, t t
163
1/2 - 1/2 _C1 + 1/2 I 3 _C1 z(t) = --------------------------} t Homogeneous Equations Many ODEs are homogeneous, meaning they can be written as functions of new variables, often making them easier to solve. For example, an ODE is homogeneous of class A if it can be written as y(x) dy =f , (4.35) dx x so that we can make the change of variable z = y/x. Maple’s genhomosol function can solve various types of homogeneous equations. > ode := -z(t) - sqrt(t^2 - z(t)^2) + t*D(z)(t)= 0: > genhomosol( ode, z(t) ); z(t) {-arctan(---------------) + ln(t) - _C1 = 0} 2 2 1/2 (t - z(t) ) Linear Equations An differential equation is linear if the coefficients of the solution, and of its derivatives, do not depend on the solution itself. For example, the equation dy + ty = 0 dt
(4.36)
is linear, whereas dy + y2 = 0 (4.37) dt is not. Maple’s linearsol function can solve first-order linear equations. > ode := diff(z(t),t) + p(t)*z(t) = q(t); /d \ ode := |-- z(t)| + p(t) z(t) = q(t) 164
\dt
/
> linearsol( ode, z(t) ); / / / \ / | | | | | {z(t) = | | q(t) exp( | p(t) dt) dt + _C1| exp( | -p(t) dt)} | | | | | \/ / / / Separable Equations An ODE is separable if the independent and dependent variables can be isolated on opposite sides of the equation. For example, dy = xy dx
(4.38)
is separable, as we can rewrite it as 1 dy = x dx, y
(4.39)
and solve by integrating both sides. Maple’s separablesol detects whether an equation is separable and if so, solves it. > ode := t^2*(z(t)+1) +z(t)^2*(t-1)*diff(z(t),t) = 0; 2 ode := t
2 (z(t) + 1) + z(t)
/d \ (t - 1) |-- z(t)| = 0 \dt /
> separablesol( ode, z(t) ); 2 {1/2 t
2 + t + ln(t - 1) + 1/2 z(t) - z(t) + ln(z(t) + 1) + _C1 = 0}
Equations with Constant Coefficients A differential equation has constant coefficients if the coefficents of the solution and its derivatives do not depend on either the dependent or inde165
pendent variables. For example, the equation d2 y dy −2 +y =0 dt2 dt
(4.40)
has constant coefficients, but dy dy −t =0 2 dt dt
(4.41)
does not. Maple’s constcoeffsols function can solve equations that have constant coefficients. > ode := 3*diff(z(t),t$3) + diff(z(t),t$2) diff(z(t),t) - 3*z(t) = 0: > constcoeffsols(ode, z(t) ); 1/2 [exp(t), exp(- 2/3 t) sin(1/3 5
1/2 t), exp(- 2/3 t) cos(1/3 5
The solution of such equations arises from the roots of a polynomial derived from the coefficients. If the roots of that polynomial cannot be found, Maple returns an implicitly defined solution. > ode := diff(z(t),t$5)+diff(z(t),t$2)-diff(z(t),t)-3*z(t): > constcoeffsols(ode, z(t) ); 2 5 [exp(RootOf(-3 - _Z + _Z + _Z , index = 5) t), 2
5 + _Z , index = 4) t),
2
5 + _Z , index = 3) t),
2
5 + _Z , index = 2) t),
2
5 + _Z , index = 1) t)]
exp(RootOf(-3 - _Z + _Z
exp(RootOf(-3 - _Z + _Z
exp(RootOf(-3 - _Z + _Z
exp(RootOf(-3 - _Z + _Z
166
t)]
Equations with Rational Coefficients A differential equation has rational coefficients if the coefficients can be written as the quotient of two polynomials of the independent variable. For example, the equation 3 d2 z 3 dz − + 2 z(t) = 0 dt2 t dt t
(4.42)
has rational coefficients. Maple’s ratsols function can solve equations with rational coefficients. > ode := (D@@2)(z)(t) - 3/t*D(z)(t) + 3/t^2*z(t): > ratsols(ode, z(t)); 3 [t, t ] Changing Variables In many cases, a change of variable may simplify an equation so that it can be solved. Maple’s dchange function, from its PDEtools package, can be used to make such a change. Its arguments are a set of equations defining the change of variable, the differential equation on which to perform the change, and a list of the new variables. The following ODE is invariant under rotation by an angle alpha. > ODE := (-y(x)+x*diff(y(x),x))/(x+y(x)*diff(y(x),x))=H(sqrt(x^2+y(x)^2)); /d \ -y(x) + x |-- y(x)| \dx / 2 2 1/2 ODE := ------------------- = H((x + y(x) ) ) /d \ x + y(x) |-- y(x)| \dx / > tr := { x = ys(xs)*sin(alpha) + xs*cos(alpha), > y(x) = cos(alpha)*ys(xs) - sin(alpha)*xs}: > newODE := simplify(dchange(tr,ODE,[xs,ys(xs)])); / d \ xs |--- ys(xs)| - ys(xs) 167
\dxs / 2 2 1/2 newODE := ------------------------ = H((ys(xs) + xs ) ) / d \ ys(xs) |--- ys(xs)| + xs \dxs /
4.3.3
Solving Systems of ODEs
In many applications, a system of ODEs must be solved. Frequently, such a system has the form Y 0 (t) = A(t)Y (t) + B(t), (4.43) where Y is a vector of solutions y1 (t), . . . , yn (t), A(t) is an n × n matrix whose entries are functions of t, and B(t) is a vector of functions of t. Such an equation can be solved using Maple’s matrixDE function. > A := matrix(2,2,[1,t^2,t,1]); [ 2] A := [1 t ] [ ] [t 1 ] > sol := matrixDE(A,t); sol := [ [ 3/2 5/2 3/2 5/2 ] [t BesselK(3/5, 2/5 t ) exp(t) , t BesselI(3/5, 2/5 t ) exp(t)] [ 5/2 5/2 ] [-exp(t) t BesselK(2/5, 2/5 t ) , exp(t) t BesselI(-2/5, 2/5 t )], [0, 0] ] Higher-order ODEs can be converted to a system of first-order ODEs. The convertsys function can be used to aid in such a conversion. > > >
deq2 := (D@@3)(y)(x) = y(x)*x: init2 := y(0) = 3, D(y)(0) = 2, (D@@2)(y)(0) = 1: convertsys(deq2, [init2], y(x), x);
[[YP[1] = Y[2], YP[2] = Y[3], YP[3] = Y[1] x], 168
2 d d [Y[1] = y(x), Y[2] = -- y(x), Y[3] = --- y(x)], 0, [3, 2, 1]] dx 2 dx
4.3.4
Displaying Solutions
Maple offers a set of functions that can be used to visualize the solution of an ODE. In particular, the DEplot function can be used to solve an ODE numerically and plot solution curves, given ranges for the dependent and independent variables to guide in the solution and plotting processes. > > > > >
DEplot(cos(x)*diff(y(x),x$3)-diff(y(x),x$2)+ Pi*diff(y(x),x)=y(x)-x,y(x), x=-2.5..1.4, [[y(0)=1,D(y)(0)=2,(D@@2)(y)(0)=1]], y=-4..5,stepsize=.05);
The solution to the above ODE is plotted in Figure 4.5. If the equation or system of equations is autonomous, i.e. the coefficients of the dependent variable and its derivatives do not depend on the independent variable, then a direction field is also plotted. The following Maple statement produces a direction field for the ODE dx/dt = −2x, which is shown in Figure 4.6. DEplot(diff(x(t),t)=-2*x(t),x(t),t=0..1,x=0..2);
4.4
Simulink
Simulink is one of Matlab’s toolboxes. It allows a user to describe to Matlab a dynamical system, which is a sytem that can be modeled by a differential equation and whose solution is dependent on time. Simulink provides an intuitive means of describing a model to Matlab that can be simulated. A mathematical model of a physical phenomenon can be viewed as a combination of several simpler processes. Using Simulink, you can link several “black boxes” together in order to create the entire model, and then have Matlab run simulate the model. 169
4
y(x) 2
–2
x
0
–1
–2
–4 Figure 4.5: Solution to ODE plotted using DEplot 170
1
2
1.5
x(t)
1
0.5
0
0.2
0.6
0.4
0.8
1
t Figure 4.6: Solution and direction field for autonomous ODE system
171
Figure 4.7: Main Simulink window
4.4.1
Invoking Simulink
Simulink can be used either through its set of commands, or through its GUI, which is recommended. To start the Simulink GUI, simply type simulink at the prompt. Matlab will then display a figure window with several choices of categories of components that can be used to construct a model. This window is shown in Figure 4.7. Using the File menu, you may create a new, empty model and begin construction.
4.4.2
Building Models in Simulink
In the main Simulink window, you may double-click on any of the icons representing categories of model components and obtain a new figure window that contains several “black boxes” from that category. You can include any of these black boxes, or blocks, in your model by dragging it into the figure window created for your own model. You may then link blocks within your model by clicking on one of the outputs of one block and dragging the mouse pointer to an input of another block.
4.4.3
Block Properties
Once a block is included in your model, you may double-click on the block to edit its properties. An example of these property dialogs is shown in Figure 4.8. After all of your blocks have been linked and their properties have been set appropriately, you may go to the Simulation menu and choose Start to run the simulation. 172
Figure 4.8: Properties of the From Workspace block
173
4.4.4
A Simulink Model
We will now use Simulink to simulate the motion of a pendulum. If there is viscous friction in the pivot and there is an applied moment M (t) around the pivot, then a nonlinear differential equation describes the angle θ(t) that the pendulum makes with the vertical at time t: l
d2 θ dθ + c + mgL sin θ = M (t), 2 dt dt
(4.44)
where l is the mass moment of inertia about the pivot, m is the mass of the pendulum, L is the length of the pendulum’s rod, and g is the gravity constant. For this simulation, we will assume that l = 4, mgL = 10, c = 0.8, and M (t) is a square wave with an amplitude of 3 and a frequency of 0.5Hz. Furthermore, since the differential equation is of second order, we will need to specify initial conditions on θ(t) and θ 0 (t). We prescribe θ(0) = π/4 and θ 0 (0) = 0. That is, the pendulum starts from a resting position, held at a 45-degree angle. Before constructing the model in Simulink, we rewrite the second-order equation as a system of two first order equations, and then rewrite them as integral equations. Introducting a new dependent variable ω(t) = θ 0 (t), we obtain Z θ(t) = ω(s) ds, Z ω(t) = −0.8ω(s) − 10 sin θ(s) + M (s) ds, (4.45) with initial conditions θ(0) = π/4,
4.4.5
ω(0) = 0.
(4.46)
Integrator blocks
An integrator block is used in a model to compute the integral of some quantity. In this case, we see that we will need two integrator blocks, one to obtain ω(t), and then a second block to integrate ω(t) to obtain θ(t). The values of ω(t) and θ(t) are then combined with the value of M (t) to produce the input back into the first integrator block, to compute ω(t) again. We will add two integrator blocks to our model, and link them together, to reflect that the output of one integration, ω(t), is the input to the other, that computes θ(t). We also edit the properties of both blocks to account for the initial conditions. Finally, we edit the labels of these blocks to indicate 174
Figure 4.9: Model with two integrator blocks that one produces ω as output, and the other produces θ. The updated model is displayed in Figure 4.9.
4.4.6
Gain blocks
Gain blocks are used to multiply some quantity by a scalar. In this model, we see that ω is multiplied by 0.8, and also that the entire integrand for ω is multiplied by 0.25 before integrating. We therefore add two gain blocks. The first has a value of 0.8, and receives the output from the first integrator, which represents ω. The second has a value of 0.25, and its output serves as the input to the integrator that yields ω. The updated model is displayed in Figure 4.10.
4.4.7
Function blocks
We observe that the quantity 10 sin θ is computed. Since θ is the output of the second integrator block, we add a function block that computes the quantity 10 sin u, where u is the input of the function block. We then set its input to be the output of the second integrator block, which represents θ. The model, with this function block added, is shown in Figure 4.11. 175
Figure 4.10: Model with two gain blocks added
Figure 4.11: Function block added to pendulum model
176
4.4.8
Sources and Sinks
A source block provides external input to the model, and a sink block receives the model’s output and stores it externally. A commonly used source is “from workspace”, which uses the value of an object in the Matlab workspace as input to the model. Similarly, the sink “to workspace” takes its input and stores it in the Matlab workspace. The source and sink windows, displaying the blocks available for each category, are shown in Figures 4.12 and 4.13, respectively. For our model, we add a source block that represents the signal M (t). By editing its properties, we can specify the shape of the wave, its amplitude, and its frequency. The properties dialog for this block is shown in Figure 4.14.
4.4.9
Sum blocks
Now, we have all of the terms in the expression that is integrated, and then multiplied by 0.25, to produce ω. To add them together, we add a sum block, available from the Math category. By editing its properties, we can dictate how many inputs the blocks should have, and whether the inputs should be added or subtracted. In this case, we will be adding one input, M (t), and subtracting two, 0.8ω and 10 sin θ. The properties dialog for the sum block is shown in Figure 4.15, and the updated model, with the sum block added, is shown in Figure 4.16.
4.4.10
The Scope Sink
Finally, we add a sink, a scope block. This block is used to plot its input, in this case θ(t), with respect to time. It provides toolbar buttons to zoom in or out, among other functionality. The completed model is shown in Figure 4.17. Once we have linked all of the blocks, we can run the simulation and observe the output. As expected, the pendulum swings back and forth, but to a lesser extent as time moves forward. However, because of the applied moment around the pivot, it does not come to rest. The plot of θ is shown in Figure 4.18.
4.5
Exercises
1. In this problem, we examine a slight variation on banded matrices: circulant matrices. These are used to solve partial differential equa177
Figure 4.12: Simulink source blocks
178
Figure 4.13: Simulink sink blocks
179
Figure 4.14: Properties of Signal Generator block
180
Figure 4.15: Properties of sum block
181
Figure 4.16: Model with sum block included
Figure 4.17: Completed Simulink model
182
Figure 4.18: Plot of simulation output θ(t) tions with periodic boundary conditions, such as the first-order wave equation ∂u ∂u =a , ∂t ∂x
t > 0, 0 < x < 1,
u(x, 0) = f (x),
0 < x < 1,
u(0, t) = u(1, t),
t > 0.
(4.47)
(4.48)
(4.49)
We are going to solve this problem using a central difference to approximate ∂u/∂x and a forward difference to approximate ∂u/∂t: ∂u(x, t) u(x, t + k) − u(x, t) ≈ , ∂t k
∂u(x, t) u(x + h, t) − u(x − h, t) ≈ ∂x 2h (4.50) In trying to compute u(x, t + k) from u(x, t), if we evaluate ∂u/∂x at time t + k instead of at time t, we obtain an implicit method for 183
computing u(x, t + k): u(x, t + k) = u(x, t) + ak
u(x + h, t + k) − u(x − h, t + k) . 2h
(4.51)
This method is called implicit because in order to obtain the values of u at time t + k from those at time t, it is necesssary to solve a system of linear equations. In this case, we have ak ak u(x − h, t + k) + u(x, t + k) − u(x + h, t + k) = u(x, t). (4.52) 2h 2h It follows that we must sovle a system Ax = b where the right-hand side b represents the values of the solution u at time t, and the coefficient matrix A has diagonal elements equal to 1, sub-diagonal elements equal to ak/2h, and super-diagonal elements equal to −ak/2h. Thus we have a banded matrix. What about the case where x is 0? Then we need the values of u at points outside of the interval [0, 1]. Since u is periodic, though, we can “wrap around” and use points at the right end of the interval to obtain the equations for u(0, t). Similarly, when computing u(1 − h, t) we must use u(0) in place of u(1). This requires us to add elements at the corners of our matrix A. In particular, if A is an n × n matrix, we must set A1n = ak/2h and An1 = −ak/2h. The matrix is no longer banded, but is called circulant because the super-diagonal and sub-diagonal both “wrap around” and continue through the opposite corners of the matrix. This is a special case of what is known as an anti-symmetric or Toeplitz‘ matrix. In this case, our matrix looks like this: 1 −ak/2h 0 ··· 0 ak/2h ak/2h 1 −ak/2h 0 ··· 0 . . . . . . . . . . . . . . 0 . . A= .. . . . . . . . . . . . . 0 . 0 ··· 0 ak/2h 1 −ak/2h −ak/2h 0 ··· 0 ak/2h 1 (4.53) We will now step through the process of constructing and then solving such a system. We will let a = 1, our time step k = 0.1, and choose N = 20 gridpoints within the interval (0, 1). They will be equally spaced, with spacing h = 1/N . 184
(a) Construct a banded matrix A of size N × N , with the diagonal equal to all 1’s, the super-diagonal equal to −ak/2h, and the sub-diagonal equal to ak/2h. (b) Fill in the corner elements A1N and AN 1 as described above to make the matrix circulant. (c) Compute the LU decomposition of A. What can you say about the structure of L and U ? Does your discussion of using LU or QR for banded matrices apply here as well? Discuss. (d) Construct a vector b representing the initial data f (x). In this case we will use f (x) = sin(2πx). Hints: first, construct a column vector x containing the gridpoints x j = jh, j = 0, . . . , N −1. The sin function in Matlab applies sin to all elements of its argument. (e) Solve the system Au1 = b. Then use u1 as the right-hand side and solve Au2 = u1 . Repeat this to obtain u1 , . . . , u5 . (f) Plot b versus x (you should have computed both vectors in part (d)). Then use Matlab’s hold command to hold that plot in place, and plot ui versus x for i = 1, . . . , 5. (g) Based on your plots, can you describe the behavior of the solution u(x, t) as time moves forward? 2. In this problem, we will use sparse matrices and iterative methods for solving a convection-diffusion problem: ∂u ∂2u ∂u =a + b 2, ∂t ∂x ∂x u(x, 0) = f (x),
0 < x < 2π,
t>0
0 < x < 2π,
u(0, t) = u(2π, t),
t>0
(4.54) (4.55) (4.56)
We can solve this problem using a method similar to the one discussed in Exercise 2 for the first-order wave equation. As in that case, we assume a time step of k and a grid of N equally spaced points, with spacing h = 2π/N . In this case, we must build a circulant matrix A, of size N × N , with the following entries: • The main diagonal has entries equal to 1 + 2b/h 2 . • The superdiagonal entries, and the lower left entry A N,1 , are all equal to ak/2h − bk/h2 185
• The subdiagonal entries, and the upper right entry A 1,N , are all equal to −ak/2h − bk/h2 In the case where a = 0, A is symmetric positive definite, but otherwise it is not. Once we have constructed A, we can compute the solution u at times t = k, 2k, 3k, . . . as follows: • Let u0 = f , where f is a vector of values of the initial data f (x) above, evaluated at the gridpoints x j = jh, for j = 0, 1, . . . , N . • Solve Auj+1 = uj , for j = 0, 1, . . .. Then uj approximates the solution u(x, jk) evaluated at each of the gridpoints. We can continue this process until we have evaluated the solution at some final time, say t = 1. (a) Write a Matlab function convecsolve that accepts three arguments, a, b, and f. The first two arguments refer to a and b above, while f is a column vector of length N + 1, with entries defined as follows: f(j) = f (jh),
h=
2π , N
j = 0, . . . , N.
(4.57)
Your function must solve the convection-diffusion problem with the given parameters and return a matrix U whose columns are defined as follows: The jth column of U, U(:,j), must hold the computed values of the solution u(x, jk), for j = 1, 2, . . . until the time jk = 1. For this problem, use k = 1/100. Thus U should be a N × 100 matrix. If the argument a is 0, you must solve the system Au k+1 = uk using the Conjugate Gradient method. Otherwise, you must use the Quasi-Minimal Residual (QMR) method. (b) Use your function to solve the convection-diffusion problem for the case where a = 0, b = 1, and f (x) = sin x, where f is to be evaluated at N = 256 equally spaced points. Plot the matrix U that is returned. How does the solution behave as time moves forward? (c) Repeat part (b), with a = −1, b = 2, f (x) = cos(4x) and N = 256. Again, how does the solution behave? (d) Modify your function as follows to produce a new function convecsolve2: 186
• Use preconditioning. If using Conjugate Gradients, use cholinc and use R’ and R as left and right preconditioners, where R is the Cholesky factor returned by cholinc. If using QMR, use luinc to produce factors L and U and use them as left and right preconditioners. In both cases, use a drop tolerance of 0.5 as the second argument to cholinc or luinc. • When using pcg or qmr, get the flag output and if its value is nonzero (i.e. the iteration failed), display a warning message and continue as before. 3. Let u(x, t) represent the temperature at any point along a thin wire at a specific time. The variable x denotes a position on the wire, with x = 0 denoting the left endpoint and x = 1 denoting the right endpoint. The variable t denotes time elapsed since the initial time t = 0. The temperature at point x at the initial time is given by the function g(x), 0 < x < 1. The temperature of the endpoints of the wire are held at a fixed value α. (a) Describe a mathematical problem that can be used to model the temperature of the wire for t > 0. (b) Solve the problem given the data g(x) = sin 2πx and α = 0 to obtain a prediction of what the temperature in the wire will be at t = 0.1. Try to obtain an answer that is accurate to within 1% relative error; i.e. try to ensure that your computed solution u ˜(t, 1) satisfies |˜ u(x, 1) − u(x, 1)| ≤ 10−2 , |u(x, 1)|
0 < x < 1.
(4.58)
Describe your approach to achieving this accuracy and explain why you believe that it is effective. 4. Consider the initial-boundary value problem with the periodic boundary conditions ut = uxx , 0 < x < 1, t > 0 (4.59) u(x, 0) = sin 4πx, u(0, t) = u(1, t),
0≤x≤1 t≥0
(4.60) (4.61)
We will solve this problem using a grid of n + 1 equally spaced points xi = i∆x, i = 0, . . . , n, where ∆x = 1/(n + 1). 187
(a) Using the method of lines, discretize this problem in space only and write down an initial value problem of the form u0 (t) = Au(t), where
u(t) =
t > 0,
u(x0 , t) u(x1 , t) .. . u(xn , t)
,
(4.62)
(4.63)
where you are to describe the entries of the matrix A. To account for these boundary conditions, we will prescribe that for any t, u(x0 , t) = u(xn+1 , t) and u(xn , t) = u(x−1 , t). (b) Solve the initial value problem from part (a) to obtain u(0.01) using the trapezoidal rule u(tn+1 ) = u(tn ) +
∆t (Au(tn ) + Au(tn+1 )). 2
(4.64)
To solve the system of linear equations at each time step, use ten iterations of the Gauss-Seidel method. Use n = 10k − 1, k = 1, . . . , 4, and ∆t = (∆x2 )/4. Submit your code and all computed solutions. (c) Explain why the Gauss-Seidel method will converge in this case, and why Gaussian elimination is not a practical method for this problem. (d) Using the trapezoidal rule on the initial value problem described in part (a) is equivalent to using what numerical method to solve the original initial-boundary value problem?
188
