MA1002 – Computational Mathematics
Module – 3: Numerical Methods - III
MA1002 – COMPUTATIONAL MATHEMATICS
Module 1: Numerical Methods III
Module 3
Numerical Methods III: Solutions of ordinary differential equations
(Picard’s Method, Euler’s and Modified Euler’s Method, Classical Fourh
Order Runge Kutta Method).
MA1002 - COMPUTATIONAL MATHEMATICS
Module 3: Numerical Methods III
Numerical Solution of First Order ODE
We have several analytical methods to solve the ordinary Differential
equations (ODE). Under this topic we will discuss the Numerical methods for
solving the ODE of first order –
𝑑𝑦/𝑑𝑥 = 𝑓(𝑥, 𝑦)
For these methods it is desired that the equation having numerical coefficients
and the Initial condition is given. i.e. by this way we will solve the Initial value
problems (IVP).
MA1002 - COMPUTATIONAL MATHEMATICS
Module 3: Numerical Methods III
• These methods are categorized in two classes-
1. Multistep Method
2. Single Step Method
• Single Step Methods are those methods by which in a single step we can
find the solution of ODE at the specific value of independent variable. i.e.
there is no need to use the previous values.
• In Multistep methods the solution of the problem be found in several steps
by using the values obtained in the last step.
MA1002 - COMPUTATIONAL MATHEMATICS
Module 3: Numerical Methods III
Picard’s Method
Let us consider the IVP as –
𝑑𝑦/𝑑𝑥 = 𝑓(𝑥, 𝑦) … (1)
with the initial condition-
𝑦 𝑥0 = 𝑦0 … (2)
The Picard’s Method for finding nth approximation of the value of ‘𝑦’ at some
specific value of ‘𝑥’ is given by-
𝑥
𝑦𝑛 = 𝑦0 + න 𝑓 𝑥, 𝑦𝑛−1 𝑑𝑥 ; 𝑛 = 1,2,3, …
𝑥0
The process of approximation is stopped when two consecutive values of 𝑦 (i.e.
𝑦𝑖−1 and 𝑦𝑖 ) are same as per desired degree of accuracy.
MA1002 - COMPUTATIONAL MATHEMATICS
Module 3: Numerical Methods III
Picard’s method can be applied to find precise accurate value, only to those
equations whose successive integrations can be performed.
Ex. Employ Picard’s method to obtain the solution of following differential
equation at 𝑥 = 0.4, correct to four places of decimals.
𝑑𝑦
= 𝑥2 + 𝑦2 ; 𝑦 0 = 0
𝑑𝑥
Sol. The given equation is
𝑑𝑦
= 𝑥 2 + 𝑦 2 = 𝑓 𝑥, 𝑦 ; 𝑦0 = 𝑦 0 = 0, 𝑥0 = 0
𝑑𝑥
By Picard’s method 𝑥
𝑦𝑛 = 𝑦0 + න 𝑓 𝑥, 𝑦𝑛−1 𝑑𝑥 ; 𝑛 = 1,2,3, …
𝑥0
MA1002 - COMPUTATIONAL MATHEMATICS
Module 3: Numerical Methods III
𝑥
𝑦𝑛 = 0 + න 𝑥 2 + 𝑦𝑛−1 2 𝑑𝑥 ; 𝑛 = 1,2,3, …
𝑥0
𝑥
⇒ 𝑦𝑛 = න 𝑥 2 + 𝑦𝑛−1 2 𝑑𝑥 ; 𝑛 = 1,2,3, … (1)
0
First Approximation (𝑛 = 1)
𝑥 𝑥
𝑥 3
𝑦1 = න 𝑥 2 + 𝑦0 2 𝑑𝑥 = න 𝑥 2 + 0 𝑑𝑥 = ;
3
0 0
(0.4)3
𝑦1 0.4 = = 0.021333333
3
MA1002 - COMPUTATIONAL MATHEMATICS
Module 3: Numerical Methods III
Second Approximation (𝑛 = 2)
𝑥 𝑥 2
𝑥 3 𝑥3 𝑥7
𝑦2 = න 𝑥 2 + 𝑦1 2 𝑑𝑥 = න 𝑥 2 + 𝑑𝑥 = + ;
3 3 63
0 0
(0.4)3 (0.4)7
𝑦2 0.4 = + = 0.021359339
3 63
∵ 𝑦2 − 𝑦1 = 0.021359339 − 0.021333333 = 0.000026006 < 𝜀
So, the accuracy is achieved. The solution of the problem correct to four
decimal places is
𝑦2 0.4 = 𝑦 0.4 = 0.02135
MA1002 - COMPUTATIONAL MATHEMATICS
Module 3: Numerical Methods III
Ex. Use Picard method to approximate 𝑦 when 𝑥 = 0.2, correct to three decimal
places, for the following equation
𝑑𝑦
+ 𝑦 = 𝑥 ;𝑦 0 = 1
𝑑𝑥
Sol. The given equation is
𝑑𝑦
= 𝑥 − 𝑦 = 𝑓 𝑥, 𝑦 ; 𝑦0 = 𝑦 0 = 1, 𝑥0 = 0
𝑑𝑥
By Picard’s method
𝑥
𝑦𝑛 = 𝑦0 + න 𝑓 𝑥, 𝑦𝑛−1 𝑑𝑥 ; 𝑛 = 1,2,3, …
𝑥0
MA1002 - COMPUTATIONAL MATHEMATICS
Module 3: Numerical Methods III
𝑥
⇒ 𝑦𝑛 = 1 + න 𝑥 − 𝑦𝑛−1 𝑑𝑥 ; 𝑛 = 1,2,3, …
0
First Approximation (𝑛 = 1)
𝑥 𝑥
𝑥2
𝑦1 = 1 + න 𝑥 − 𝑦0 𝑑𝑥 = 1 + න 𝑥 − 1 𝑑𝑥 = 1 − 𝑥 + ;
2
0 0
2
0.2
𝑦1 0.2 = 1 − 0.2 + = 0.82
2
Second Approximation (𝑛 = 2)
𝑥 𝑥
𝑥2 𝑥 3
𝑦2 = 1 + න 𝑥 − 𝑦1 𝑑𝑥 = 1 + න 𝑥 − 1 + 𝑥 − 𝑑𝑥 = 1 − 𝑥 + 𝑥 2 − ;
2 6
0 0
MA1002 - COMPUTATIONAL MATHEMATICS
Module 3: Numerical Methods III
3
0.2
𝑦2 0.2 = 1 − 0.2 + 0.2 2 − = 0.83867
6
∵ 𝑦2 − 𝑦1 = 0.83867 − 0.82 = 0.01867 ≰ 𝜀
So, the accuracy is not achieved.
Third Approximation (𝑛 = 3)
𝑥 𝑥
𝑥 3
𝑦3 = 1 + න 𝑥 − 𝑦2 𝑑𝑥 = න 𝑥 − 1 − 𝑥 + 𝑥 2 − 𝑑𝑥
6
0 0
𝑥 3 𝑥 4
= 1 − 𝑥 + 𝑥2 − +
3 24
3 4
0.2 0.2
𝑦3 0.2 = 1 − 0.2 − 0.2 2 − + = 0.8374
3 24
MA1002 - COMPUTATIONAL MATHEMATICS
Module 3: Numerical Methods III
∵ 𝑦3 − 𝑦2 = 0.8374 − 0.83867 = 0.00127 ≰ 𝜀
So, the accuracy is not achieved.
Fourth Approximation (𝑛 = 4)
𝑥 𝑥
𝑥 3 𝑥 4
𝑦4 = 1 + න 𝑥 − 𝑦3 𝑑𝑥 = න 𝑥 − 1 − 𝑥 + 𝑥 2 − + 𝑑𝑥
3 24
0 0
𝑥 3 𝑥 4 𝑥 5
= 1 − 𝑥 + 𝑥2 − + −
3 12 120
0.2 3 0.2 4 0.2 5
𝑦4 0.2 = 1 − 0.2 − 0.2 2 − + − = 0.837464
3 12 120
∵ 𝑦4 − 𝑦3 = 0.837464 − 0.8374 = 0.000064 < 𝜀
Hence the solution correct to three decimal places is y = 0.837.
MA1002 - COMPUTATIONAL MATHEMATICS
Module 3: Numerical Methods III
Euler’s Method
The differential equation is –
𝑑𝑦
= f x, y ; y 𝑥0 = 𝑦0 … (1)
𝑑𝑥
Let 𝑥 = 𝑎 be the point where the solution of the differential equation is
required.
Let us first divide the interval 𝑥0 , 𝑎 into finite number of subintervals (say 𝑛)
by taking the step size
𝑎 − 𝑥0
ℎ=
𝑛
The approximate value of 𝑦 at 𝑥 = 𝑥𝑛 by Euler’s method is given by
𝒚𝒏 = 𝒚𝒏−𝟏 + 𝒉𝒇 𝒙𝒏−𝟏 , 𝒚𝒏−𝟏 ; 𝒏 = 𝟏, 𝟐, 𝟑, …
MA1002 - COMPUTATIONAL MATHEMATICS
Module 3: Numerical Methods III
Ex. Find the solution of the following differential equation at 𝑥 = 0.1.
𝑑𝑦 𝑦 − 𝑥
= ; 𝑦 = 1 𝑎𝑡 𝑥 = 0
𝑑𝑥 𝑦 + 𝑥
𝑦−𝑥
Sol. Here 𝑥0 = 0, 𝑦0 = 1, 𝑓 𝑥, 𝑦 = , 𝑎 = 0.1
𝑦+𝑥
By the Euler’s method
𝑦𝑛 = 𝑦𝑛−1 + ℎ𝑓 𝑥𝑛−1 , 𝑦𝑛−1 ; 𝑛 = 1,2,3, …
𝑦𝑛−1 − 𝑥𝑛−1
⇒ 𝑦𝑛 = 𝑦𝑛−1 + ℎ ; 𝑛 = 1,2,3, … (1)
𝑦𝑛−1 + 𝑥𝑛−1
Let ℎ = 0.02 then we have
𝑥1 = 0.02, 𝑥2 = 0.04, 𝑥3 = 0.06, 𝑥4 = 0.08, 𝑥5 = 0.1
So, we need 𝑦5 = 𝑦 𝑥5 = 𝑦(0.1).
MA1002 - COMPUTATIONAL MATHEMATICS
Module 3: Numerical Methods III
Calculation of solution by Euler’s Formula stated in (1)
𝑦0 − 𝑥0 1−0
𝑦1 = 𝑦0 + ℎ = 1 + 0.02 = 1.02 = 𝑦(0.02)
𝑦0 + 𝑥0 1+0
𝑦1 − 𝑥1 1.02 − 0.02
𝑦2 = 𝑦1 + ℎ = 1.02 + 0.02 = 1.03923 = 𝑦(0.04)
𝑦1 + 𝑥1 1.02 + 0.02
𝑦2 − 𝑥2 1.03923 − 0.04
𝑦3 = 𝑦2 + ℎ = 1.03923 + 0.02 = 1.057747
𝑦2 + 𝑥2 1.03923 + 0.04
= 𝑦(0.06)
𝑦3 − 𝑥3 1.057747 − 0.06
𝑦4 = 𝑦3 + ℎ = 1.057747 + 0.02 = 1.0755998
𝑦3 + 𝑥3 1.057747 + 0.06
= 𝑦(0.08)
𝑦4 − 𝑥4 1.0755998 − 0.08
𝑦5 = 𝑦4 ℎ = 1.0755998 + 0.02 = 1.092830675
𝑦4 + 𝑥4 1.0755998 + 0.08
= 𝑦(0.1)
MA1002 - COMPUTATIONAL MATHEMATICS
Module 3: Numerical Methods III
Ex: Use Euler’s method (taking ℎ = 0.2) to find approximate value of 𝑦 when 𝑥 =
0.6 for the following differential equation
𝑑𝑦
= 1 − 2𝑥𝑦 ; 𝑦 0 = 0
𝑑𝑥
Sol. Here 𝑥0 = 0, 𝑦0 = 0, 𝑓 𝑥, 𝑦 = 1 − 2𝑥𝑦
As, ℎ = 0.2 so, 𝑥1 = 0.2, 𝑥2 = 0.4, 𝑥3 = 0.6
By Euler’s method
𝑦𝑛 = 𝑦𝑛−1 + ℎ𝑓 𝑥𝑛−1 , 𝑦𝑛−1 ; 𝑛 = 1,2,3, …
⇒ 𝑦𝑛 = 𝑦𝑛−1 + ℎ 1 − 2𝑥𝑛−1 𝑦𝑛−1 ; 𝑛 = 1,2,3, … … (1)
𝑦1 = 𝑦0 + ℎ 1 − 2𝑥0 𝑦0 = 0 + 0.2 1 − 2(0)(0) = 0.2 = 𝑦(0.2)
𝑦2 = 𝑦1 + ℎ 1 − 2𝑥1 𝑦1 = 0.2 + 0.2 1 − 2(0.2)(0.2) = 0.384 = 𝑦 0.4
𝑦3 = 𝑦2 + ℎ 1 − 2𝑥2 𝑦2 = 0.384 + 0.2 1 − 2(0.4)(0.384) = 0.52256 = 𝑦 0.6
MA1002 - COMPUTATIONAL MATHEMATICS
Module 3: Numerical Methods III
Modified Euler’s Method
The given differential equation is
𝑑𝑦
= f x, y ; y 𝑥0 = 𝑦0
𝑑𝑥
The First approximation of 𝑦 at 𝑥 = 𝑥𝑛 is given by
1
𝑦𝑛 = 𝑦𝑛−1 + ℎ𝑓 𝑥𝑛−1 , 𝑦𝑛−1 ; 𝑛 = 1,2,3, …
and the 𝑚𝑡ℎ modification (or approximation) of 𝑦 at 𝑥 = 𝑥𝑛 is given by
𝑚 ℎ 𝑚−1
𝑦𝑛 = 𝑦𝑛−1 + 𝑓 𝑥𝑛−1 , 𝑦𝑛−1 + 𝑓 𝑥𝑛 , 𝑦𝑛
2
𝑚 = 2,3,4, … . ; 𝑛 = 1,2,3, …
MA1002 - COMPUTATIONAL MATHEMATICS
Module 3: Numerical Methods III
Ex. Find the solution of the following differential by Modified Euler’s method at
𝑥 = 0.6, taking step size of 0.2.
𝑑𝑦
= 𝑥 + 𝑦 ; 𝑦 = 1 𝑎𝑡 𝑥 = 0
𝑑𝑥
Sol. Here
𝑥0 = 0, 𝑦0 = 1, 𝑓 𝑥, 𝑦 = 𝑥 + 𝑦 , 𝑥 = 0.6, ℎ = 0.2
∴ 𝑥1 = 0.2, 𝑥2 = 0.4, 𝑥3 = 0.6
By Modified Euler’s method, the First approximation of 𝑦 at 𝑥 = 𝑥𝑛 is given by
1
𝑦𝑛 = 𝑦𝑛−1 + ℎ 𝑥𝑛−1 + 𝑦𝑛−1 ; 𝑛 = 1,2,3, …
and the 𝑚𝑡ℎ approximation of y at 𝑥 = 𝑥𝑛 is given by
𝑚 ℎ 𝑚−1
𝑦𝑛 = 𝑦𝑛−1 + 𝑥𝑛−1 + 𝑦𝑛−1 + 𝑥𝑛 + 𝑦𝑛 ;
2
𝑚 = 2,3,4, … , 𝑛 = 1,2,3, …
MA1002 - COMPUTATIONAL MATHEMATICS
Module 3: Numerical Methods III
1 1
𝑦1 = 𝑦0 + ℎ 𝑥0 + 𝑦0 = 1.2 = 𝑦 (0.2)
2 ℎ 1
𝑦1 = 𝑦0 + 𝑥0 + 𝑦0 + 𝑥1 + 𝑦1 = 1.22954
2
3 ℎ 2
𝑦1 = 𝑦0 + 𝑥0 + 𝑦0 + 𝑥1 + 𝑦1 = 1.23088
2
4 ℎ 3
𝑦1 = 𝑦0 + 𝑥0 + 𝑦0 + 𝑥1 + 𝑦1 = 1.230945
2
5 ℎ 4
𝑦1 = 𝑦0 + 𝑥0 + 𝑦0 + 𝑥1 + 𝑦1 = 1.230947
2
4 5
As 𝑦1 = 𝑦1 = 1.23094 (up to five places of decimal)
So 𝑦1 = 𝑦 0.2 = 1.23094 (up to five places of decimal)
MA1002 - COMPUTATIONAL MATHEMATICS
Module 3: Numerical Methods III
1
𝑦2 = 𝑦1 + ℎ 𝑥1 + 𝑦1 = 1.492835 = 𝑦 1 (0.4)
2 ℎ 1
𝑦2 = 𝑦1 + 𝑥1 + 𝑦1 + 𝑥2 + 𝑦2 = 1.524069
2
3 ℎ 2
𝑦2 = 𝑦1 + 𝑥1 + 𝑦1 + 𝑥2 + 𝑦2 = 1.52534
2
4 ℎ 3
𝑦2 = 𝑦1 + 𝑥1 + 𝑦1 + 𝑥2 + 𝑦2 = 1.525392
2
5 ℎ 4
𝑦2 = 𝑦1 + 𝑥1 + 𝑦1 + 𝑥2 + 𝑦2 = 1.525392
2
4 5
As 𝑦2 = 𝑦2 = 1.525392(up to six places of decimal)
So 𝑦2 = 𝑦 0.4 = 1.525392(up to six places of decimal)
MA1002 - COMPUTATIONAL MATHEMATICS
Module 3: Numerical Methods III
For the value of y3
1
𝑦3 = 𝑦2 + ℎ 𝑥2 + 𝑦2 = 1.852405 = 𝑦 1 (0.6)
2 ℎ 1
𝑦3 = 𝑦2 + 𝑥2 + 𝑦2 + 𝑥3 + 𝑦3 = 1.88500
2
3 ℎ 2
𝑦3 = 𝑦2 + 𝑥2 + 𝑦2 + 𝑥3 + 𝑦3 = 1.886194
2
4 ℎ 3
𝑦3 = 𝑦2 + 𝑥2 + 𝑦2 + 𝑥3 + 𝑦3 = 1.886237
2
5 ℎ 4
𝑦3 = 𝑦2 + 𝑥2 + 𝑦2 + 𝑥3 + 𝑦3 = 1.886239
2
4 5
As 𝑦3 = 𝑦3 = 1.88623 (up to five places of decimal)
So 𝑦3 = 𝑦 0.6 = 1.88623 (up to five places of decimal)
Which is the required result.
MA1002 - COMPUTATIONAL MATHEMATICS
Module 3: Numerical Methods III
Classical Fourth Order Runge –Kutta Method:
𝑑𝑦
= f x, y ; y 𝑥0 = 𝑦0 … (1)
𝑑𝑥
Let h denote the length of the interval between equi-distant values of x.
If 𝑥1 = 𝑥0 + ℎ, then the corresponding value of 𝑦1 is
1
𝑦1 = 𝑦0 + [𝑘1 + 2 𝑘2 + 𝑘3 + 𝑘4 ]
6
𝑘1 = ℎ𝑓(𝑥0 , 𝑦0 )
𝑘2 = ℎ𝑓(𝑥0 + ℎ2, 𝑦0 + 𝑘21 )
𝑘3 = ℎ𝑓(𝑥0 + ℎ2, 𝑦0 + 𝑘22 )
𝑘4 = ℎ𝑓(𝑥0 + ℎ, 𝑦0 + 𝑘3 )
MA1002 - COMPUTATIONAL MATHEMATICS
Module 3: Numerical Methods III
Similarly, the value of 𝑦2 is given by -
1
𝑦2 = 𝑦1 + [𝑘1 + 2 𝑘2 + 𝑘3 + 𝑘4 ]
6
𝑘1 = ℎ𝑓 𝑥1 , 𝑦1 ; 𝑘2 = ℎ𝑓(𝑥1 + ℎ2, 𝑦1 + 𝑘21)
𝑘3 = ℎ𝑓 𝑥1 + ℎ2, 𝑦1 + 𝑘22 ; 𝑘4 = ℎ𝑓(𝑥1 + ℎ, 𝑦1 + 𝑘3 )
and in general, the formula for 𝑦𝑛 is –
1
𝑦𝑛 = 𝑦𝑛−1 + 𝑘1 + 2 𝑘2 + 𝑘3 + 𝑘4 ; 𝑛 = 1,2,3, …
6
𝑘1 = ℎ𝑓(𝑥𝑛−1 , 𝑦𝑛−1 )
𝑘2 = ℎ𝑓(𝑥𝑛−1 + ℎ2, 𝑦𝑛−1 + 𝑘21)
𝑘3 = ℎ𝑓(𝑥𝑛−1 + ℎ2, 𝑦𝑛−1 + 𝑘22)
𝑘4 = ℎ𝑓(𝑥𝑛−1 + ℎ, 𝑦𝑛−1 + 𝑘3 )
MA1002 - COMPUTATIONAL MATHEMATICS
Module 3: Numerical Methods III
Ex. Use Runge – Kutta Method of fourth order to solve the following equation
𝑎𝑡 𝑥 = 0.2.
𝑦 ′ = 𝑥 + 𝑦 ; 𝑦(0) = 1
Sol. Here it is given that 𝑥0 = 0, 𝑦0 = 1, 𝑓(𝑥, 𝑦) = 𝑥 + 𝑦
Let us take ℎ = 0.1, so 𝑥1 = 0.1 𝑎𝑛𝑑 𝑥2 = 0.2
By Runge-Kutta Method of fourth order-
1
𝑦𝑛 = 𝑦𝑛−1 + 𝑘1 + 2 𝑘2 + 𝑘3 + 𝑘4 ; 𝑛 = 1,2,3, …
6
Where, 𝑘1 = ℎ𝑓(𝑥𝑛−1 , 𝑦𝑛−1 )
𝑘2 = ℎ𝑓(𝑥𝑛−1 + ℎ2, 𝑦𝑛−1 + 𝑘21)
𝑘3 = ℎ𝑓(𝑥𝑛−1 + ℎ2, 𝑦𝑛−1 + 𝑘22)
𝑘4 = ℎ𝑓(𝑥𝑛−1 + ℎ, 𝑦𝑛−1 + 𝑘3 )
MA1002 - COMPUTATIONAL MATHEMATICS
Module 3: Numerical Methods III
So, for n=1
1
𝑦1 = 𝑦0 + [𝑘1 + 2 𝑘2 + 𝑘3 + 𝑘4 ]
6
and where, 𝑘1 = ℎ𝑓 𝑥0 , 𝑦0 =h[𝑥0 + 𝑦0 ] = 0.1 (0 + 1) = 0.1
𝑘2 = ℎ𝑓 𝑥0 + ℎ2, 𝑦0 + 𝑘21 = h 𝑥0 + ℎ2 + 𝑦0 + 𝑘21 = 0.1 ( 0.05 + 1.05)
= 0.110
𝑘3 = ℎ𝑓 𝑥0 + ℎ2, 𝑦0 + 𝑘22 = ℎ 𝑥0 + ℎ2 + 𝑦0 + 𝑘22
=0.1 [ 0.05 + 1.055] = 0.1105
𝑘4 = ℎ𝑓 𝑥0 + ℎ, 𝑦0 + 𝑘3 = ℎ 𝑥0 + ℎ + 𝑦0 + 𝑘3
= 0.1 [0.1 + 1.105] = 0.12105
1
So 𝑦1 = 1 + [0.1 + 2 0.110 + 0.1105 + 0.12105]
6
⇒ 𝑦1 = 𝑦 0.1 = 1.11034
MA1002 - COMPUTATIONAL MATHEMATICS
Module 3: Numerical Methods III
for n = 2
1
𝑦2 = 𝑦1 + [𝑘1 + 2 𝑘2 + 𝑘3 + 𝑘4 ]
6
where, 𝑘1 = ℎ𝑓 𝑥1 , 𝑦1 =h 𝑥1 + 𝑦1 = 0.1 (0.1 + 1.11034) = 0.12103
𝑘2 = ℎ𝑓 𝑥1 + ℎ2, 𝑦1 + 𝑘21 = h 𝑥1 + ℎ2 + 𝑦1 + 𝑘21
= 0.1 0.15 + 1.17081 = 0.13208
𝑘3 = ℎ𝑓 𝑥1 + ℎ2, 𝑦1 + 𝑘22 = ℎ 𝑥1 + ℎ2 + 𝑦1 + 𝑘22
= 0.1 (0.15 + 1.1763) = 0.13263
𝑘4 = ℎ𝑓 𝑥1 + ℎ, 𝑦1 + 𝑘3 = ℎ 𝑥1 + ℎ + 𝑦1 + 𝑘3
= 0.1 (0.2 + 1.24293) = 0.144293
1
So 𝑦2 = 1.11034 + [0.12103 + 20.13208 + 0.13263 + 0.144293]
6
⇒ 𝑦2 = 𝑦 0.2 = 1.242803
MA1002 - COMPUTATIONAL MATHEMATICS
Module 3: Numerical Methods III
Ex. Use Classical fourth order Runge – Kutta method to find the solution of
following ODE at 𝑥 = 0.4 by taking ℎ = 0.2.
𝑑𝑦
= −2𝑥𝑦 2 ; 𝑦 0 = 1
𝑑𝑥
Sol. Here, 𝑥0 = 0, 𝑦0 = 1, 𝑓 𝑥, 𝑦 = −2𝑥𝑦 2 , ℎ = 0.2
So, 𝑥1 = 0.2, 𝑥2 = 0.4 and we have to find 𝑦2 = 𝑦 𝑥2 = 𝑦(0.4)
By Runge-Kutta Method of fourth order-
1
𝑦𝑛 = 𝑦𝑛−1 + 𝑘1 + 2 𝑘2 + 𝑘3 + 𝑘4 ; 𝑛 = 1,2,3, …
6
Where, 𝑘1 = ℎ𝑓(𝑥𝑛−1 , 𝑦𝑛−1 )
𝑘2 = ℎ𝑓(𝑥𝑛−1 + ℎ2, 𝑦𝑛−1 + 𝑘21 )
𝑘3 = ℎ𝑓(𝑥𝑛−1 + ℎ2, 𝑦𝑛−1 + 𝑘22 )
𝑘4 = ℎ𝑓(𝑥𝑛−1 + ℎ, 𝑦𝑛−1 + 𝑘3 )
MA1002 - COMPUTATIONAL MATHEMATICS
Module 3: Numerical Methods III
So, for n=1
1
𝑦1 = 𝑦0 + [𝑘1 + 2 𝑘2 + 𝑘3 + 𝑘4 ]
6
and where, 𝑘1 = ℎ𝑓 𝑥0 , 𝑦0 =ℎ[−2𝑥0 𝑦0 2 ] =0
ℎ 𝑘1 ℎ 𝑘1 2
𝑘2 = ℎ𝑓 𝑥0 + +,𝑦
2 0
= ℎ −2 𝑥0 +
2 2
𝑦0 + 2
= 0.2 −2 0 + 0.1 1 + 0 2 = −0. 04
ℎ 𝑘2 ℎ 𝑘2 2
𝑘3 = ℎ𝑓 𝑥0 + 2 , 𝑦0
+ = ℎ −2 𝑥0 +
2 𝑦0 +
2 2
=0.2 −2 0 + 0.1 1 − 0.02 2 = −0.038416
𝑘4 = ℎ𝑓 𝑥0 + ℎ, 𝑦0 + 𝑘3 = ℎ −2 𝑥0 + ℎ 𝑦0 + 𝑘3 2
= 0.2 −2 0 + 0.2 1 − 0.038416 2 = −0.073797
1
So, 𝑦1 = 𝑦(0.2) = 1 + [0 + 2 −0.04 − 0.0398416 − 0.073797]
6
⇒ 𝑦1 = 𝑦 0.1 = 0.96108663
MA1002 - COMPUTATIONAL MATHEMATICS
Module 3: Numerical Methods III
1
for n = 2, 𝑦2 = 𝑦1 + [𝑘1 + 2 𝑘2 + 𝑘3 + 𝑘4 ]
6
where, 𝑘1 = ℎ𝑓 𝑥1 , 𝑦1 =ℎ[−2𝑥1 𝑦1 2 ]
= 0.2 −2(0.2)(0.96108663)2 = −0.073895
ℎ 𝑘1 2
𝑘2 = ℎ𝑓 𝑥1 + ,𝑦
2 1
+ 𝑘21
= ℎ −2 𝑥1 + ℎ
2
𝑦1 + 2
2 = −0. 102484
= 0.2 −2 0.2 + 0.1 0.96108663 − 0.0369475
𝑘2 2
𝑘3 = ℎ𝑓 𝑥1 + ℎ2, 𝑦1
+ 𝑘2
= ℎ −2
2
𝑥1 + ℎ2
𝑦0 + 2
= 0.2 −2 0.2 + 0.1 0.96108663 − 0.051242 2 = −0.099338
𝑘4 = ℎ𝑓 𝑥1 + ℎ, 𝑦1 + 𝑘3 = ℎ −2 𝑥1 + ℎ 𝑦0 + 𝑘3 2
= 0.2 −2 0.2 + 0.2 0.96108663 − 0.099338 2 = −0.1188177
1
So, 𝑦2 = 0.96108663 + −0.073895 + 2 −0. 102484 − 0.099338 − 0.1188177
6
⇒ 𝑦2 = 𝑦 0.2 = 0.861694
MA1002 – COMPUTATIONAL MATHEMATICS
Module 1: Numerical Methods III
Summary
In this topic, you learnt:
Concept of Numerical solution of ordinary differential
equation of first order.
Applications of numerical methods (Picard method,
Euler’s and Modified Euler’s method, and Classical Fourth
order Runge Kutta Method) to solve ordinary differential
equations of order one.
THANK YOU