Transcript File

Numerical Solutions of
Ordinary Differential
Equations
4/8/2015
Dr.G.Suresh Kumar@KL
University
1
Ordinary Differential Equations
Equations which are composed of an unknown function
and its derivatives are called differential equations.
dv
c
g v
dt
m
v - dependent variable
t - independent variable
Differential equations play a fundamental role in
engineering because many physical phenomena are best
formulated mathematically in terms of their rate of
change.
4/8/2015
Dr.G.Suresh Kumar@KL University
2
Types of Differential Equations
• When a function involves one independent
variable, the equation is called an ordinary
differential equation (ODE).
•A partial differential equation (PDE)
involves two or more independent
variables.
4/8/2015
Dr.G.Suresh Kumar@KL University
3
Model Ordinary Differential Equations
4/8/2015
Dr.G.Suresh Kumar@KL University
4
Why Numerical solutions of ODE ?
Many differential equations cannot be
solved analytically, in which case we have
to satisfy ourselves with an approximation
to the solution.
4/8/2015
Dr.G.Suresh Kumar@KL University
5
Initial value problem
Initial value problem associated with first
order differential equations is
dy
 f(x,y)
dx
y(x0 )  y 0 .
4/8/2015
Dr.G.Suresh Kumar@KL University
6
Methods of solving IVP
Taylor series
 Picard successive approximation
 Euler
 Modified Euler
 Runge-Kutta (RK)

4/8/2015
Dr.G.Suresh Kumar@KL University
7
Taylor Series Method
The Taylor series method is a straight
forward adaptation of classic calculus to
develop the solution as an infinite series.
 The method is not strictly a numerical
method but it is used in conjunction with
numerical schemes.

4/8/2015
Dr.G.Suresh Kumar@KL University
8
Taylor Series Method
The Taylor series expansion of y(x) about
x = x0 is

x - x 0 2
yx   yx 0   x - x 0  yx 0  
yx 0 
2!

x - x0 

3
3!

x - x0 
yx 0  
4
4!
4/8/2015
y
(IV)
x 0   ...
Dr.G.Suresh Kumar@KL University
9
Problem#1
The differential equation y'(x) = x + y satisfying
y(0)=1. Determine the Taylor series solution.
Also compare with the exact solution.
{Analytical solution: y(x)=2ex - x -1(Cauchy
linear differential equation).}
4/8/2015
Dr.G.Suresh Kumar@KL University
10
Solution
First, we find derivatives of y' = x + y,
yx   x  yx 
yx   1  yx 
yx   yx 
y
(IV)
y0   x  y0   0  1  1
y0   1  y0   1  1  2
 y0  y0  2
x   yx 
y (IV) 0   y0   2.
4/8/2015
Dr.G.Suresh Kumar@KL University
11
Solution
The Taylor series expansion of y(x) about
x = x0 is

x - x 0 2
yx   yx 0   x - x 0  yx 0  
2!

x - x 0 3

3!

x - x 0 4
yx 0  
4!
yx 0 
y (IV) x 0   ...
x2
x3
x4
2  2  2  ...
yx   1  x 1 
2!
3!
4!
3
4
x
x
1 x  x2 

 ...
3
12
4/8/2015
Dr.G.Suresh Kumar@KL University
12
Comparison of solutions
The exact solution is


x 2 x3 x 4
y(x) 21  x 


 ...  (x  1)
2!
3!
4!


3
4
x
x
2
 1 x  x 

 ...
3
24
Therefore, Taylor series solution and exact
solution are equal.
4/8/2015
Dr.G.Suresh Kumar@KL University
13
Picard’s Method
The initial value Problem
dy
 f(x,y)
dx
y(x0 )  y 0
is equivalent to the following integral equation
x
y(x)  y 0 
 f(t, y(t))dt .
x0
4/8/2015
Dr.G.Suresh Kumar@KL University
14
Picard’s Method
The Picard’s nth approximate solution is
y n (x)  y0 
x
 f(t,yn 1(t))dt,
n  1, 2, 3, ...
x0
yn (x)  y(x)
as
4/8/2015
n  .
Dr.G.Suresh Kumar@KL University
15
Problem#1
Using Picard’s process of successive
approximation, obtain a solution up to the
fifth approximation of the equation
dy
xy
dx
such that y = 1 when x = 0. Check your
answer by finding the exact solution.
Ans.
3
4
5
6
x
x
x
x
y5 (x)  1  x  x 2    
3 12 60 720
Exact solution y(x) = 2ex – (x+1).
4/8/2015
Dr.G.Suresh Kumar@KL University
16
Application Problem
For the free-falling bungee jumper with
linear drag, compute the velocity at 10
seconds of a free-falling parachutist of
mass 80kg and drag coefficient 10 kg/s at
10 seconds, with initial velocity 20 m/s,
using
(i) Taylor series method
(ii) Picard’s method.
4/8/2015
Dr.G.Suresh Kumar@KL University
17
Solution
Let v(t) be the velocity and ‘a’ be the
acceleration of the parachutist at any time t.
Then the total force acting on the parachutist is
equal to the sum of the
 downward force FD due to gravity and
 upward force FU due to air resistance.
i.e. F = FD + FU
FD = mg
FU α v (linear drag)
FU = - Cd v(t).
4/8/2015
Dr.G.Suresh Kumar@KL University
18
Solution
By Newton second law of motion F = ma.
ma = mg – Cd v
Therefore,
Cd
dv
g
v
dt
m
Given that m = 80kg, Cd = 10 kg/s,
and initial velocity v(0)=20 m/s.
Now solve the above IVP using Taylor and
Picard’s method.
4/8/2015
Dr.G.Suresh Kumar@KL University
19
Euler’s Method
Problem:
Given the first order ODE :
with the initial condition:
Determine:
y ' ( x )  f ( x, y )
y0  y ( x0 )
yi  y ( x0  ih)
for i  1,2,...
Euler Method:
y0  y ( x0 )
yi 1  yi  h f ( xi , yi )
for i  1,2,...
4/8/2015
Dr.G.Suresh Kumar@KL University
20
Interpretation of Euler Method
y2
y1
y0
x0
x1
x2
4/8/2015
x
Dr.G.Suresh Kumar@KL University
21
Interpretation of Euler Method
y1
Slope=f(x0,y0)
y1=y0+hf(x0,y0)
hf(x0,y0)
y0
x0
h
x1
x2
4/8/2015
x
Dr.G.Suresh Kumar@KL University
22
Interpretation of Euler Method
y2
y2=y1+hf(x1,y1)
Slope=f(x1,y1)
hf(x1,y1)
Slope=f(x0,y0)
y1
y1=y0+hf(x0,y0)
hf(x0,y0)
y0
x0
h
x1
h
x2
4/8/2015
x
Dr.G.Suresh Kumar@KL University
23
Interpretation of Euler Method
dy
 f ( x, y )
dx
Solution :
yi 1  yi  f ( xi , yi )  h
4/8/2015
Dr.G.Suresh Kumar@KL University
24
Problem#1
Using Euler’s method, find an approximate
value of y corresponding to x=1, given
that yʹ = x + y such that y=1 when x=0.
Ans: Taking h = 0.1, then
y(0.1) = 1.1, y(0.2) = 1.22, y(0.3) = 1.36
y(0.4) = 1.53, y(0.5) = 1.72, y(0.6) = 1.94
y(0.7) = 2.19, y(0.8) = 2.48, y(0.9) = 2.81
y(1.0) = 3.18.
4/8/2015
Dr.G.Suresh Kumar@KL University
25
Problem#2
A cup of a coffee originally has temperature
of 68oC. The temperature of the ambient
is 21oC and the thermal constant is
0.017oC/min. Determine the temperature
of the coffee from t = 0 to 10 minutes
insteps of 2 minutes.
Sol. The differential equation is
Tʹ(t) = -k(T - Ta) = -0.017(T - 21)
T(0) = 68
4/8/2015
Dr.G.Suresh Kumar@KL University
26
Improvements of Euler’s method
A fundamental source of error in Euler’s
method is that the derivative at the
beginning of the interval is assumed to
apply across the entire interval.
Two simple modifications are available to
circumvent this shortcoming:
 Heun’s (Modified Euler) Method
 Midpoint (Improved Polygon) Method

4/8/2015
Dr.G.Suresh Kumar@KL University
27
Modified Euler’s Method
To improve the estimate of the slope, determine two
derivatives for the interval:
(1) At the initial point
(2) At the end point
The two derivatives are then averaged to obtain
an improved estimate of the slope for the entire
interval.
Predictor:
Corrector :
yi01  yi  f ( xi , yi )h
yi 1
f ( xi , yi )  f ( xi 1 , yi01 )
 yi 
h
2
4/8/2015
Dr.G.Suresh Kumar@KL University
28
Euler Modified method

To improve the estimate of the
slope, determine two derivatives
for the interval:
 At the initial point
 At the end point

Predictor:
The two derivatives are then
averaged to obtain an improved
estimate of the slope for the entire
interval.
yi01  yi  f ( xi , yi )h
f ( xi , yi )  f ( xi 1 , yi01 )
Corrector : yi 1  yi 
h
2
29
Heun’s method (improved)
Original Huen’s:
Predictor:
yi01  yi  f ( xi , yi )h
Corrector :
f ( xi , yi )  f ( xi 1, yi01 )
1
yi 1  yi 
h
2
Note that the corrector can be iterated to
improve the accuracy of yi+1 .
Predictor: y  yi  f ( xi , yi )h
0
i 1
j 1
f
(
x
,
y
)

f
(
x
,
y
j
i
i
i 1
i 1 )
Corrector : y i1  yi 
h
2
j  1,2,...
30
Problem#1
Using Modified Euler’s method, find an
approximate value of y when x= 0.3, given
that yʹ = x + y and y=1 when x=0.
Sol. Compare the IVP with the general IVP
yʹ = f(x, y), y(x0) = x0, we have
f(x, y) = x + y, x0 = 0 and y0 = 1.
Taking the step size h = 0.1, then x0 = 0, x1 =
0.1, x2 = 0.2, x3 = 0.3. Now we find the
corresponding value of y, using modified
Euler’s method.
31
Solution
Step-1. To find y1 = y(x1) = y(0.1)
Predictor formula (Euler’ formula)
y1(0) = y0 + h f(x0, y0)
= 1 + 0.1 f(0, 1) = 1+0.1(0+1) =1.1
Corrector formula
y1(k) = y0 + (h/2) [f(x0, y0) + f(x1, y1(k-1))]
y1(1) = y0 + (h/2) [f(x0, y0) + f(x1, y1(0))]
= 1 + (0.1/2) [f(0,1) + f(0.1, 1.1)]
= 1 + 0.05 [0+1 + 0.1+1.1] =1.11
32
Solution
y1(2) = y0 + (h/2) [f(x0, y0) + f(x1, y1(1))]
= 1 + (0.1/2) [f(0,1) + f(0.1, 1.11)]
= 1 + 0.05 [0+1 + 0.1+1.11]
=1.1105.
y1(3) = y0 + (h/2) [f(x0, y0) + f(x1, y1(2))]
= 1 + (0.1/2) [f(0,1) + f(0.1, 1.1105)]
= 1 + 0.05 [0+1 + 0.1+1.1105]
=1.1105.
Therefore, y1 = 1.1105.
33
Solution
Step-2. To find y2 = y(x2) = y(0.2)
Predictor formula (Euler’ formula)
y2(0) = y1 + h f(x1, y1)
= 1.1105 + 0.1 f(0.1, 1.1105)
= 1.1105+0.1(0.1 + 1.1105) =1.2316.
Corrector formula
y2(k) = y1 + (h/2) [f(x1, y1) + f(x2, y2(k-1))]
y2(1) = y1 + (h/2) [f(x1, y1) + f(x2, y2(0))]
= 1.2316 + (0.1/2) [f(0.1,1.1105) + f(0.2, 1.2316)]
= 1.2316 + 0.05 [0.1+1.1105 + 0.2+1.2316]
34
=1.2426.
Solution
y2(2) = y1 + (h/2) [f(x1, y1) + f(x2, y2(1))]
= 1.1105 + (0.1/2) [f(0.1,1.1105) + f(0.2, 1.2426)]
= 1.1105 + 0.05 [0.1+1.1105 + 0.2+1.2426]
=1.2432.
y2(3) = y1 + (h/2) [f(x1, y1) + f(x2, y2(2))]
= 1.1105 + (0.1/2) [f(0.1,1.1105) + f(0.2, 1.2432)]
= 1.1105 + 0.05 [0.1+1.1105 + 0.2+1.2432]
=1.2432.
Therefore, y2 = y(0.2)= 1.2432.
Similarly, y3 = y(0.3) = 1.4004.
35
Application Problem

A storage tank (Fig.) contains a liquid at
depth y where y = 0 when the tank is half
full. Liquid is withdrawn at a constant flow
rate Q=500m3/h to meet demands. The
contents are resupplied at a sinusoidal
rate 3Q sin2(t). The surface area of the
tank is 1200 m2. Determine the depth of
the tank from t= 0 to t=6 in steps of 2
hours, using modified Euler’s method.
4/8/2015
Dr.G.Suresh Kumar@KL University
36
Problem#1
3Q
Q =500 m3/h
sin2(t)
4/8/2015
Dr.G.Suresh Kumar@KL University
37
Solution
Let y(t) be the depth of the tank at any time t
and A be the surface area of the tank.
The depth of the tank is changed when
volume of the liquid in the tank changed.
Therefore,
 Rate of change 

  inflowrate– out flow rate
 in volume 
d
Ay(t)   3Q sin 2 (t)  Q
dt
4/8/2015
Dr.G.Suresh Kumar@KL University
38
Solution
Since A is constant, it implies that
dy 3Q 2
Q

sin (t) 
dt
A
A
Given that Q=500m3/h, A =1200 m2 and y(0) = 0.
Substituting these value, we have
dy 5 2
5
 sin (t)
dt 4
12
y(0) 0
Now solving the above initial value problem,
using modified Euler’s method.
4/8/2015
Dr.G.Suresh Kumar@KL University
39
Runge - Kutta Method
4/8/2015
Dr.G.Suresh Kumar@KL University
40
Carl Runge
(1856-1927)
Martin Wilhelm Kutta
(1867-1944)
4/8/2015
Dr.G.Suresh Kumar@KL University
41
Introduction
Runge-Kutta methods are very popular because of
their good efficiency; and are used in most
computer programs for differential equations.
 These methods agree with Taylor’s series solution
up to the term in hr, where r differs method to
method and is called the order of the method.
 Euler’s and modified Euler’s are called the first
and second order Runge-Kutta methods.
 The fourth order R-K method is most commonly
used and is often referred to as R-K method only.

4/8/2015
Dr.G.Suresh Kumar@KL University
42
Working Rule

To find increment ‘k’ of y corresponding to an
increment h of x by R-K method from the initial
value problem
y´ = f(x, y) with y(x0) = y0 is as follows
Calculate
k1 = h f(x0, y0)
k2 = h f(x0+0.5h, y0+0.5k1)
k3 = h f(x0+0.5h, y0+0.5k2)
k4 = h f(x0+h, y0+k3)
4/8/2015
Dr.G.Suresh Kumar@KL University
43
Working Rule
Finally compute
k = 1/6 (k1 + 2k2 + 2k3 + k4)
is the weighted mean of k1, k2, k3 and k4.
 Therefore the required approximate value
is
y1 = y(x1) = y(x0 + h) = y0 + k.

4/8/2015
Dr.G.Suresh Kumar@KL University
44
Graphical Representation
k2
k4
k3
k
k1
k
xi
xi + h/2
4/8/2015
1
k1  2k2  2k3  k4 
6
xi + h
Dr.G.Suresh Kumar@KL University
45
Problem#1
Using Runge-Kutta method of fourth
order, solve yʹ = (y2 – x2)/(y2 + x2) with
y(0)=1 at x = 0.2, 0.4.
Sol : - Compare the given problem with the
general first order initial value problem
yʹ=f(x, y) satisfying y(x0)=y0, we have
f(x, y) = (y2 – x2)/(y2 + x2)
x0 = 0, y0 =1.
Here we have to find value of y at x=0.2, 0.4.

4/8/2015
Dr.G.Suresh Kumar@KL University
46
Solution
Taking step size h = 0.2.
Step(1) To find y(0.2)
k1 = h f(x0, y0)
= 0.2 f(0, 1) = 0.2 (12-02)/(12+02)
= 0.2
k2 = h f(x0+0.5h, y0+0.5k1)
= 0.2 f(0+0.5(0.2), 1+ 0.5(0.2))
= 0.2 f(0.1, 1.1)
= 0.19672
4/8/2015
Dr.G.Suresh Kumar@KL University
47
Solution
k3 = h f(x0+0.5h, y0+0.5k2)
= 0.2 f(0+0.5(0.2), 1+0.5(0.19672))
= 0.2 f(0.1, 1.09836)
= 0.1967
k4 = h f(x0+h, y0+k3)
= 0.2 f(0+0.2, 1+0.1967)
= 0.2 f(0.2, 1.1967)
= 0.1891
4/8/2015
Dr.G.Suresh Kumar@KL University
48
Solution
k = 1/6 (k1 + 2k2 + 2k3 + k4)
= 1/6 (0.2 +2(0.19672)+2(0.1967)+0.1891)
= 0.19599.
Hence y1 = y0 + k
y(0.2) = 1 + 0.19599 = 1.19599.
Step(2) To find y(0.4)
k1 = h f(x1, y1)
= 0.2 f(0.2, 1.19599)
= 0.1891
4/8/2015
Dr.G.Suresh Kumar@KL University
49
solution
k2 = h f(x1+0.5h, y1+0.5k1)
= 0.2 f(0.2+0.5(0.2), 1.19599+ 0.5(0.1891))
= 0.2 f(0.3, 1.29054)
= 0.1795
k3 = h f(x1+0.5h, y1+0.5k2)
= 0.2 f(0.2+0.5(0.2), 1.19599+0.5(0.1795))
= 0.2 f(0.3, 1.28574)
= 0.1793
4/8/2015
Dr.G.Suresh Kumar@KL University
50
Solution
k4 = h f(x1+h, y1+k3)
= 0.2 f(0.2+0.2, 1.19599+0.1793)
= 0.2 f(0.4, 1.37529)
= 0.1688
k = 1/6 (k1 + 2k2 + 2k3 + k4)
= 1/6 (0.1891 +2(0.1795)+2(0.1793)+0.1688)
= 0.1792.
Hence y2 = y1 + k
y(0.4) = 1.19599 + 0.1792 = 1.37519.
4/8/2015
Dr.G.Suresh Kumar@KL University
51
Problem#2
Apply Runge-Kutta method to find an
approximate value of y for x=0.2 in steps
of 0.1, if yʹ = x + y2 given that y=1 where
x = 0.
Sol : - Compare the given problem with the
general first order initial value problem
yʹ=f(x, y) satisfying y(x0)=y0, we have
f(x, y) = x + y2 , x0 = 0, y0 =1.
Here we have to find value of y at x = 0.2.

4/8/2015
Dr.G.Suresh Kumar@KL University
52
Solution
Taking step size h = 0.1.
Step(1) To find y(0.1)
k1 = h f(x0, y0)
= 0.1 f(0, 1) = 0.1 (0 + 12)
= 0.1
k2 = h f(x0+0.5h, y0+0.5k1)
= 0.1 f(0+0.5(0.1), 1+ 0.5(0.1))
= 0.2 f(0.05, 1.05)
= 0.1152
4/8/2015
Dr.G.Suresh Kumar@KL University
53
Solution
k3 = h f(x0+0.5h, y0+0.5k2)
= 0.1 f(0+0.5(0.1), 1+0.5(0.1152))
= 0.1 f(0.05, 1.0576)
= 0.11685
k4 = h f(x0+h, y0+k3)
= 0.1 f(0+0.1, 1+0.11685)
= 0.1 f(0.1, 1.11685)
= 0.13473
4/8/2015
Dr.G.Suresh Kumar@KL University
54
Solution
k = 1/6 (k1 + 2k2 + 2k3 + k4)
= 1/6 (0.1 +2(0.1152)+2(0.11685)+0.13473)
= 0.11647 ≈ 0.1165.
Hence y1 = y0 + k
y(0.1) = 1 + 0.1165 = 1.1165.
Step(2) To find y(0.2)
k1 = h f(x1, y1)
= 0.1 f(0.1, 1.1165)
= 0.1347
4/8/2015
Dr.G.Suresh Kumar@KL University
55
solution
k2 = h f(x1+0.5h, y1+0.5k1)
= 0.1 f(0.1+0.5(0.1), 1.1165+ 0.5(0.1347))
= 0.2 f(0.15, 1.1838)
= 0.1551
k3 = h f(x1+0.5h, y1+0.5k2)
= 0.1 f(0.1+0.5(0.1), 1.1165+0.5(0.1551))
= 0.1 f(0.15, 1.194)
= 0.1575
4/8/2015
Dr.G.Suresh Kumar@KL University
56
Solution
k4 = h f(x1+h, y1+k3)
= 0.1 f(0.1+0.1, 1.1165+0.1576)
= 0.1 f(0.2, 1.2741)
= 0.1823
k = 1/6 (k1 + 2k2 + 2k3 + k4)
= 1/6 (0.1347 +2(0.1551)+2(0.1576)+0.1823)
= 0.1571.
Hence y2 = y1 + k
y(0.2) = 1.1165 + 0.1571 = 1.2736.
4/8/2015
Dr.G.Suresh Kumar@KL University
57