Numerical Analysis by Dr. Anita Pal Assistant Professor Department of Mathematics National ...

15 Pages·2016·0.93 MB·English

Numerical Analysis
by
Dr. Anita Pal
Assistant Professor
Department of Mathematics
National Institute of Technology Durgapur
Durgapur-713209
email: [email protected]
1
.
Chapter 8
Numerical Solution of Ordinary Differential Equations
Module No. 1
Runge-Kutta Methods
......................................................................................
One of the most useful tools of applied mathematics is differential equation, it may
be ordinary differential equation (ODE) or partial differential equation (PDE). Both
ODEandPDEarewidelyusedtomodellotsofmathematicalandengineeringproblems.
Unfortunately,itisnotpossibletofindanalyticalsolutionofallODEsorPDEs. Finding
of analytic solution of an ODE or a PDE is a very difficult task. But, several numerical
techniques are available to solve ODEs and PDEs.
Let us consider the following first order differential equation
dy
= f(x,y) (1.1)
dx
with initial condition
y(x ) = y . (1.2)
0 0
If the analytic solution of an ODE is not available then we will go for numerical
solution. The numerical solution of a differential equation can be determined in one of
the following two forms:
(i) A power series solution for y in terms of x. Then the values of y can be obtained
by substituting x = x ,x ,...,x .
0 1 n
(ii) A set of tabulated values of y for x = x ,x ,...,x with spacing h.
0 1 n
The general solution of an ODE of order n contains n arbitrary constants. To find
such constants, n conditions are required. The problems in which all the conditions are
specified at the initial point only, are called initial value problems (IVPs). The
ODEs of order two or more and for which the conditions are specified at two or more
points are called boundary value problems (BVPs).
There may not exist a solution of an ordinary differential equation always. The
sufficient condition for existence of unique solution is provided by Lipschitz.
Themethodstofindapproximatesolutionofaninitialvalueproblemarecalledfinite
difference methods or discrete variable methods. The solutions are determined
at a set of discrete points called a grid or mesh of points. That is, a given differential
equation is converted to a discrete algebraic equation. By solving such equation, we get
an approximate solution of the given differential equation.
1
...........................................................Runge-Kutta Methods
A finite difference method is convergent if the solution of the finite difference equa-
tion approaches to a limit as the size of the grid spacing tends to zero. But, in general,
there is no guarantee, that this limit corresponds to the exact solution of the differential
equation.
Mainly, two types of numerical methods are used, viz. explicit method and implicit
method. If the value of y depends only on the values of y , h and f(x ,y ) then the
i+1 i i i
method is called explicit method, otherwise the method is called implicit method.
In this chapter, we discuss some useful methods to solve ODEs.
1.1 Euler’s method
Euler’s method is the most simple but crude method to solve the differential equation
of the form
dy
= f(x,y), y(x ) = y . (1.3)
0 0
dx
Let x = x +h, where h is small. Then by Taylor’s series
1 0
(cid:16)dy(cid:17) h2(cid:16)d2y(cid:17)
y = y(x +h)=y +h + ,
1 0 0 dx x0 2 dx2 c1
where c lies between x and x
1 0
h2
=y +hf(x ,y )+ y(cid:48)(cid:48)(c ) (1.4)
0 0 0 1
2
Thus, by neglecting second order term, we get
y = y +hf(x ,y ). (1.5)
1 0 0 0
In general,
y = y +hf(x ,y ), n = 0,1,2,... (1.6)
n+1 n n n
This is a very slow method. To find a reasonable accurate solution, the value of
h must be taken as small. The Euler’s method is less efficient in practical problems
because if h is not sufficiently small then this method gives inaccurate result.
This method is modified to get more better result as follows:
The value of y obtained by Euler’s method is repeatedly modified to get a more
1
accurate result. For this purpose, the term f(x ,y ) is replaced by the average of
0 0
(0) (0)
f(x ,y ) and f(x +h,y ), where y = y +hf(x ,y ).
0 0 0 1 1 0 0 0
2
1.1. Euler’s method..............................................................
That is, the modified Euler’s method is
h
(1) (0)
y =y + [f(x ,y )+f(x ,y )] (1.7)
1 0 2 0 0 1 1
(0)
where y = y +hf(x ,y ).
1 0 0 0
Note that the value of y depends on the value of y obtained in previous iteration.
1 1
So this method is known as implicit method, whereas Euler method is explicit method.
Other simple methods, viz. Taylor’s series, Picard, Runge-Kutta, etc. are also avail-
able to solve the differential equation of the form (1.1).
Now, we describe Runge-Kutta method. In this method, several function evaluations
are required at each step and it avoids the computation of higher order derivatives.
There are several types of Runge-Kutta methods, such as second, third, fourth, fifth,
etc. order. The fourth-order Runge-Kutta method is more popular. These are single-
step explicit methods.
1.2 Second-order Runge-Kutta method
First we deduce second order Runge-Kutta method from modified Euler’s method.
The modified Euler’s method is
h
(0)
y = y + [f(x ,y )+f(x ,y )] (1.8)
1 0 2 0 0 1 1
(0)
where y = y +hf(x ,y ).
1 0 0 0
(0)
Now, we substitute the value of y to the equation (1.8). Then
1
h
y = y + [f(x ,y )+f(x +h,y +hf(x ,y ))].
1 0 0 0 0 0 0 0
2
Let
k =hf(x ,y ) and
1 0 0
k =hf(x +h,y +hf(x ,y )) = hf(x +h,y +k ). (1.9)
2 0 0 0 0 0 0 1
Using this notation the equation (1.8) becomes
1
y = y + (k +k ). (1.10)
1 0 1 2
2
3
...........................................................Runge-Kutta Methods
This method is known as second-order Runge-Kutta method. Note that this is
an explicit formula, whereas modified Euler’s method is an implicit method.
The local truncation error of this method is O(h3).
This formula can also be derived independently.
Independent derivation of second-order Runge-Kutta method
Suppose the solution of the equation
dy
= f(x,y), y(x ) = y (1.11)
0 0
dx
is of the following form:
y = y +ak +bk , (1.12)
1 0 1 2
where
k =hf(x ,y ) and
1 0 0
k =hf(x +αh,y +βk ),a,b,α and β are constants.
2 0 0 1
By Taylor’s theorem
h2 h3
y =y(x +h) = y +hy(cid:48) + y(cid:48)(cid:48)+ y(cid:48)(cid:48)(cid:48)+···
1 0 0 0 2 0 6 0
h2(cid:104)(cid:16)∂f(cid:17) (cid:16)∂f(cid:17) (cid:105)
=y +hf(x ,y )+ +f(x ,y ) +O(h3)
0 0 0 0 0
2 ∂x (x0,y0) ∂y (x0,y0)
(cid:20) (cid:21)
df ∂f ∂f
since = +f(x,y)
dx ∂x ∂y
k =hf(x +αh,y +βk )
2 0 0 1
(cid:104) (cid:16)∂f(cid:17) (cid:16)∂f(cid:17) (cid:105)
=h f(x ,y )+αh +βk +O(h2)
0 0 1
∂x (x0,y0) ∂y (x0,y0)
(cid:16)∂f(cid:17) (cid:16)∂f(cid:17)
=hf(x ,y )+αh2 +βh2f(x ,y ) +O(h3).
0 0 0 0
∂x (x0,y0) ∂y (x0,y0)
Substituting these values to the equation (1.12), we get
h2
y +hf(x ,y )+ [f (x ,y )+f(x ,y )f (x ,y )]+O(h3)
0 0 0 x 0 0 0 0 y 0 0
2
= y +(a+b)hf(x ,y )+bh2[αf (x ,y )+βf(x ,y )f (x ,y )]+O(h3).
0 0 0 x 0 0 0 0 y 0 0
4
1.1. Euler’s method..............................................................
Equating the coefficients of f, f and f both sides we get the following equations.
x y
1 1
a+b = 1, bα = and bβ = . (1.13)
2 2
Obviously,α = β. Here,numberofequationsisthreeandvariablesisfour. Therefore,
thesystemofequationshasmanysolutions. However,usuallytheparametersarechosen
1
as α = β = 1, then a = b = .
2
For this set of parameters, the equation (1.12) becomes
1
y = y + (k +k )+O(h3), (1.14)
1 0 1 2
2
where k = hf(x ,y ) and k = hf(x +h,y +k ). (1.15)
1 0 0 2 0 0 1
1.3 Fourth-order Runge-Kutta method
The fourth-order Runge-Kutta formula to calculate y is
1
1
y = y + (k +2k +2k +k ) (1.16)
1 0 1 2 3 4
6
where
k =hf(x ,y )
1 0 0
k =hf(x +h/2,y +k /2)
2 0 0 1
k =hf(x +h/2,y +k /2)
3 0 0 2
k =hf(x +h,y +k ).
4 0 0 3
Using the value of y one can find the value of y as
1 2
1
y = y + (k +2k +2k +k ) (1.17)
2 1 1 2 3 4
6
where
k =hf(x ,y )
1 1 1
k =hf(x +h/2,y +k /2)
2 1 1 1
k =hf(x +h/2,y +k /2)
3 1 1 2
k =hf(x +h,y +k ).
4 1 1 3
5
...........................................................Runge-Kutta Methods
In general, 1 (i) (i) (i) (i)
y = y + (k +2k +2k +k ) (1.18)
i+1 i 6 1 2 3 4
where
(i)
k =hf(x ,y )
1 i i
(i) (i)
k =hf(x +h/2,y +k /2)
2 i i 1
(i) (i)
k =hf(x +h/2,y +k /2)
3 i i 2
(i) (i)
k =hf(x +h,y +k ),
4 i i 3
for i = 0,1,2,....
Example 1.1 Given y(cid:48) = x−y2 with x = 0,y = 1. Find y(0.2) by second and fourth-
order Runge-Kutta methods.
Solution. Here, x = 0,y = 1,f(x,y) = x−y2.
0 0
By Second order Runge-Kutta method
Let h = 0.1.
k =hf(x ,y ) = 0.1×(0−12) = −0.10000.
1 0 0
k =hf(x +h,y +k ) = 0.1×f(0.1,0.9)
2 0 0 1
=−0.07100.
Now, we calculate,
1
y =y(0.1) = y(0)+ (k +k )
1 1 2
2
=1+0.5×(−0.1−0.071) = 0.91450.
To determine y = y(0.2),x = 0.1 and y = 0.91450.
2 1 1
k =hf(x ,y ) = 0.1×{0.1−(0.91450)2} = −0.073633.
1 1 1
k =hf(x +h,y +k )
2 1 1 1
=0.1×[0.2−(0.9145−0.07363)2] = −0.05071.
Therefore,
1
y =y(0.2) = y(0.1)+ (k +k )
2 1 2
2
=0.91450+0.5×(−0.07363−0.05071) = 0.85233.
6
1.1. Euler’s method..............................................................
By fourth-order Runge-Kutta methods:
Here h = 0.1,x = 0,y = 1,f(x,y) = x−y2.
0 0
Then,
k =hf(x ,y ) = 0.1×(0−(1)2) = −0.1.
1 0 0
k =hf(x +h/2,y +k /2)
2 0 0 1
0.1 0.1
=0.1×{( )−(1− )2} = −0.08525.
2 2
k =hf(x +h/2,y +k /2)
3 0 0 2
0.1 0.08525
=0.1×{( )−(1− )2} = −0.08666.
2 2
0.1 0.08666
k =0.1×{( )−(1− )2} = −0.07342.
4
2 2
Therefore,
1
y =y(0.1) = y + (k +2k +2k +k )
1 0 1 2 3 4
6
1
=1+ [−0.1000+2×{(−0.08525)+(−0.08666)}−0.07342] = 0.91379.
6
To find, y = y(0.2),h = 0.1,x = 0.1,y = 0.91379.
2 1 1
Then,
k =hf(x ,y ) = 0.1×{0.1−(0.91379)2} = −0.07350.
1 1 1
k =hf(x +h/2,y +k /2)
2 1 1 1
0.1 0.07350
=0.1×{(0.1+ )−(0.91379− )2} = −0.06192.
2 2
k =hf(x +h/2,y +k /2)
3 1 1 2
0.1 0.06192
=0.1×{(0.1+ )−(0.91379− )2} = −0.06294.
2 2
k =hf(x +h,y +k )
4 1 1 3
=0.1×{(0.1+0.1)−(0.91379−0.06294)2} = −0.05239.
Therefore,
1
y =y(0.2) = y(0)+ (k +2k +2k +k )
2 1 2 3 4
6
1
=0.91379+ [−0.07350+2×{(−0.06192)+(−0.06294)}−0.05239]
6
=0.85119.
7
...........................................................Runge-Kutta Methods
Note 1.1 The fourth order Runge-Kutta method gives better result, though, this method
has some disadvantages. In this method, a lot of function calculations are required.
Thus, if the function f(x,y) is complicated, then the Runge-Kutta method is very labo-
rious.
Error in fourth-order Runge-Kutta method
The fourth-order Runge-Kutta formula can be written as
(h/2)(cid:104) 4(k +k ) (cid:105)
2 3
y = y + k + +k .
1 0 1 4
3 2
This form is similar to the Simpson’s 1/3 formula with step size h/2. Thus, the local
h5
truncation error of this formula is − yiv(c ), i.e. of O(h5). After n steps, the
1
2880
accumulated error is
−(cid:88)n h5 yiv(c ) = −xn−x0yiv(c)h4 = O(h4).
i
2880 5760
i=1
1.4 Runge-Kutta method for a pair of equations
The second and fourth-order Runge-Kutta methods can also used to solve a pair of
first order differential equations.
Let us consider a pair of first-order differential equations
dy
=f(x,y,z)
dx
(1.19)
dz
=g(x,y,z)
dx
with initial conditions
x = x , y(x ) = y , z(x ) = z . (1.20)
0 0 0 0 0
Here, x is the independent variable and y,z are the dependent variables. So, in this
problem we have to determine the values of y,z for different values of x.
8