6 ORDINARY DIFFERENTIAL EQUATIONS.docx

1、6 ORDINARY DIFFERENTIAL EQUATIONS6 ORDINARY DIFFERENTIAL EQUATIONSDifferential equations are mathematical descriptions of how the variables and their derivatives (rates of change) with respect to one or more independent variable affect each other in a dynamical way. Their solutions show us how the d

2、ependent variable(s) will change with the independent variable(s). Many problems in natural sciences and engineering fields are formulated into a scalar differential equation or a vector differential equationthat is, a system of differential equations. In this chapter, we look into several methods o

3、f obtaining the numerical solutions to ordinary differential equations (ODEs) in which all dependent variables (x) depend on a single independent variable (t). First, the initial value problems (IVPs) will be handled with several methods including RungeKutta method and predictorcorrector methods in

4、Sections 6.1 to 6.5. The final section (Section 6.6) will introduce the shooting method and the finite difference and collocation method for solving the two-point boundary value problem (BVP). ODEs are called an IVP if the values x(t0) of dependent variables are given at the initial point t0 of the

5、independent variable, while they are called a BVP if the values x(t0) and x(tf) are given at the initial/final points t0 and tf.6.1 EULERS METHODWhen talking about the numerical solutions to ODEs, everyone starts with the Eulers method, since it is easy to understand and simple to program. Even thou

6、gh its low accuracy keeps it from being widely used for solving ODEs, it gives us a clue to the basic concept of numerical solution for a differential equation simply and clearly. Lets consider a first-order differential equation:It has the following form of analytical solution:which can be obtained

7、 by using a conventional method or the Laplace transform technique. However, such a nice analytical solution does not exist for every differential equation; even if it exists, it is not easy to find even by using a computer equipped with the capability of symbolic computation. That is why we should

8、study the numerical solutions to differential equations.Then, how do we translate the differential equation into a form that can easily be handled by computer? First of all, we have to replace the derivative y(t) = dy/dt in the differential equation by a numerical derivative (introduced before), whe

9、re the step-size h is determined based on the accuracy requirements and the computation time constraints. Eulers method approximates the derivative in Eq. (6.1.1) with Eq. (5.1.2) asand solves this difference equation step-by-step with increasing t by h each time from t = 0.This is a numeric sequenc

10、e y(kh), which we call a numerical solution of Eq. (6.1.1).To be specific, let the parameters and the initial value of Eq. (6.1.1) be a = 1, r = 1, and y0 = 0. Then, the analytical solution (6.1.2) becomesy(t) = 1 eat (6.1.5)and the numerical solution (6.1.4) with the step-size h = 0.5 and h = 0.25

11、are as listed in Table 6.1 and depicted in Fig. 6.1. We make a MATLAB program “nm610.m”, which uses Eulers method for the differential equation (6.1.1), actually solving the difference equation (6.1.3) and plots the graphs of the numerical solutions in Fig. 6.1. %nm610: Euler method to solve a 1st-o

12、rder differential equationclear, clfa = 1; r = 1; y0 = 0; tf = 2;t = 0:0.01:tf; yt = 1 - exp(-a*t); %Eq.(6.1.5): true analytical solutionplot(t,yt,k), hold onklasts = 8 4 2; hs = tf./klasts;y(1) = y0;for itr = 1:3 %with various step size h = 1/8,1/4,1/2klast = klasts(itr); h = hs(itr); y(1)=y0;for k

13、 = 1:klasty(k + 1) = (1 - a*h)*y(k) +h*r; %Eq.(6.1.3):plot(k - 1 k*h,y(k) y(k+1),b, k*h,y(k+1),ro)if k 4, pause; endendendFigure 6.1 Examples of numerical solution obtained by using the Eulers method.The graphs seem to tell us that a small step-size helps reduce the error so as to make the numerical

14、 solution closer to the (true) analytical solution. But, as will be investigated thoroughly in Section 6.2, it is only partially true. In fact, a too small step-size not only makes the computation time longer (proportional as 1/h), but also results in rather larger errors due to the accumulated roun

15、d-off effect. This is why we should look for other methods to decrease the errors rather than simply reduce the step-size.Eulers method can also be applied for solving a first-order vector differential equationwhich is equivalent to a high-order scalar differential equation. The algorithm can be des

16、cribed byand is cast into the MATLAB routine “ode_Euler()”.function t,y = ode_Euler(f,tspan,y0,N)%Eulers method to solve vector differential equation y(t) = f(t,y(t)% for tspan = t0,tf and with the initial value y0 and N time stepsif nargin4 | N = 0, N = 100; endif nargin3, y0 = 0; endh = (tspan(2)

17、- tspan(1)/N; %stepsizet = tspan(1)+0:N*h; %time vectory(1,:) = y0(:); %always make the initial value a row vectorfor k = 1:Ny(k + 1,:) = y(k,:) +h*feval(f,t(k),y(k,:); %Eq.(6.1.7)end6.2 HEUNS METHOD: TRAPEZOIDAL METHODAnother method of solving a first-order vector differential equation like Eq. (6.

18、1.6) comes from integrating both sides of the equation.If we assume that the value of the (derivative) function f(t ,y) is constant as f(tk,y(tk) within one time step tk,tk+1), this becomes Eq. (6.1.7) (with h = t k+1 t k), amounting to Eulers method. If we use the trapezoidal rule (5.5.3), it becom

19、esfunction t,y = ode_Heun(f,tspan,y0,N)%Heun method to solve vector differential equation y(t) = f(t,y(t)% for tspan = t0,tf and with the initial value y0 and N time stepsif nargin4 | N = 0, N = 100; endif nargin3, y0 = 0; endh = (tspan(2) - tspan(1)/N; %stepsizet = tspan(1)+0:N*h; %time vectory(1,:

20、) = y0(:); %always make the initial value a row vectorfor k = 1:Nfk = feval(f,t(k),y(k,:); y(k+1,:) = y(k,:)+h*fk; %Eq.(6.2.3)y(k+1,:) = y(k,:) +h/2*(fk +feval(f,t(k+1),y(k+1,:); %Eq.(6.2.4)endBut, the right-hand side (RHS) of this equation has yk+1, which is unknown at tk. To resolve this problem,

21、we replace the yk+1 on the RHS by the following Eulers approximation:so that it becomesThis is Heuns method, which is implemented in the MATLAB routine “ode_Heun()”. It is a kind of predictor-and-corrector method in that it predicts the value of yk+1 by Eq. (6.2.3) at tk and then corrects the predic

22、ted value by Eq. (6.2.4) at tk+1. The truncation error of Heuns method is O(h2) (proportional to h2) as shown in Eq. (5.6.1), while the error of Eulers method is O(h).6.3 RUNGEKUTTA METHODAlthough Heuns method is a little better than the Eulers method, it is still not accurate enough for most real-w

23、orld problems. The fourth-order RungeKutta (RK4) method having a truncation error of O(h4) is one of the most widely used methods for solving differential equations. Its algorithm is described below. where Equation (6.3.1) is the core of RK4 method, which may be obtained by substituting Simpsons rul

24、e (5.5.4)into the integral form (6.2.1) of differential equation and replacing fk+1/2 with the average of the successive function values (fk2 + fk3)/2. Accordingly, the RK4 method has a truncation error of O(h4) as Eq. (5.6.2) and thus is expected to work better than the previous two methods.The fou

25、rth-order RungeKutta (RK4) method is cast into the MATLAB routine “ode_RK4()”. The program “nm630.m” uses this routine to solve Eq. (6.1.1) with the step size h = (tf t0)/N = 2/4 = 0.5 and plots the numerical result together with the (true) analytical solution. function t,y = ode_RK4(f,tspan,y0,N,va

26、rargin)%Runge-Kutta method to solve vector differential eqn y(t) = f(t,y(t)% for tspan = t0,tf and with the initial value y0 and N time stepsif nargin 4 | N = 0, N = 100; endif nargin 3, y0 = 0; endy(1,:) = y0(:); %make it a row vectorh = (tspan(2) - tspan(1)/N; t = tspan(1)+0:N*h;for k = 1:Nf1 = h*

27、feval(f,t(k),y(k,:),varargin:); f1 = f1(:); %(6.3.2a)f2 = h*feval(f,t(k) + h/2,y(k,:) + f1/2,varargin:); f2 = f2(:);%(6.3.2b)f3 = h*feval(f,t(k) + h/2,y(k,:) + f2/2,varargin:); f3 = f3(:);%(6.3.2c)f4 = h*feval(f,t(k) + h,y(k,:) + f3,varargin:); f4 = f4(:); %(6.3.2d)y(k + 1,:) = y(k,:) + (f1 + 2*(f2

28、+ f3) + f4)/6; %Eq.(6.3.1)end%nm630: Heun/Euer/RK4 method to solve a differential equation (d.e.)clear, clftspan = 0 2;t = tspan(1)+0:100*(tspan(2) - tspan(1)/100;a = 1; yt = 1 - exp(-a*t); %Eq.(6.1.5): true analytical solutionplot(t,yt,k), hold ondf61 = inline(-y + 1,t,y); %Eq.(6.1.1): d.e. to be s

29、olvedy0 = 0; N = 4;t1,ye = oed_Euler(df61,tspan,y0,N);t1,yh = ode_Heun(df61,tspan,y0,N);t1,yr = ode_RK4(df61,tspan,y0,N);plot(t,yt,k, t1,ye,b:, t1,yh,b:, t1,yr,r:)plot(t1,ye,bo, t1,yh,b+, t1,yr,r*)N = 1e3; %to estimate the time for N iterationstic, t1,ye = ode_Euler(df61,tspan,y0,N); time_Euler = to

30、ctic, t1,yh = ode_Heun(df61,tspan,y0,N); time_Heun = toctic, t1,yr = ode_RK4(df61,tspan,y0,N); time_RK4 = tocComparison of this result with those of Eulers method (“ode_Euler()”) and Heuns method (“ode_Heun()”) is given in Fig. 6.2, which shows that the RK4 method is better than Heuns method, while

31、Eulers method is the worst in terms of accuracy with the same step-size. But,Figure 6.2 Numerical solutions for a first-order differential equation.in terms of computational load, the order is reversed, because Eulers method, Heuns method, and the RK4 method need 1, 2, and 4 function evaluations (calls) per iteration, respectively.(cf) Note that a function call takes much more time than a multiplication and thus the number of function calls should be a criterion in estimating and comparing computational time.The MATLAB built-in routines “ode23()” and “ode45()” implement the RungeKutt