Calculo Numérico II

IF-392

Semana 03
Dr. Pierre Giovanny Ramos Apestegui

Cálculo Numérico II-IF392

 Runge-Kutta Method
They are a family of methods developed from the work of the Germans Carl David Tolmé
Runge, 1856-1927, and Martin Wilhelm Kutta, 1867-1944.

They achieve the precision of the Taylor method without the need to calculate
high-order derivatives.
Taylor methods have the property of a higher-order local truncation error, but the
disadvantage of requiring the calculation and evaluation of the derivatives of f(t,y).
This is somewhat slow and complicated in most problems, which is why in
practice they are almost not used.
Euler's method unfortunately requires a very small step for reasonable accuracy.
Cálculo Numérico II-IF392
 Higher-Order Runge-Kutta
Higher order Runge-Kutta methods are available.

Derived similar to second-order Runge-Kutta.

Higher order methods are more accurate but

require more calculations.

Cálculo Numérico II-IF392

3rd Order Runge-Kutta RK3

Know as RK3
k1 = f ( xi , yi )
h 1
k2 = f ( xi + , yi + k1h )
2 2
k3 = f ( xi + h, yi − k1h + 2k2h )

yi +1 = yi + (k1 + 4k2 + k3 )
h
6
Local error is O ( h 4 ) and Global error is O ( h 3 )

Cálculo Numérico II-IF392

4th Order Runge-Kutta
One of the most popular Runga-Kutta formulas, and most used due to its high
precision, is that of order 4, which is presented below.
The geometric figure of the formula scheme of order 4

Each of the ks represents a slope. The result is a weighted average of those slopes.
Also known as the classical Runge-Kutta method, abbreviated as RK4, with local
truncation error O(h5), and global error O(h4). The price to pay for the
improvement in the error is a greater number of evaluations of the function,
resulting in greater calculation time if the function is complicated. The advantage,
over the Taylor method of order 4 (whose global error is also O(h4), is that it does
not require the calculation of the derivatives of f.
Cálculo Numérico II-IF392
4th Order Runge-Kutta RK4

k1 = f ( xi , yi )
h 1
k2 = f ( xi + , yi + k1h )
2 2
h 1
k3 = f ( xi + , yi + k2h )
2 2
k4 = f ( xi + h, yi + k3h )

yi +1 = yi + (k1 + 2k2 + 2k3 + k4 )

h
6
5 4
Local error is O ( h ) and global error is O ( h )
Cálculo Numérico II-IF392
Higher-Order Runge-Kutta
Método RK de quinto orden de Butcher RK5
k1 = f ( xi , yi )
1 1
k2 = f ( xi + h, yi + k1h )
4 4
1 1 1
k3 = f ( xi + h, yi + k1h + k2h )
4 8 8
1 1
k4 = f ( xi + h, yi − k2h + k3h )
2 2
3 3 9
k5 = f ( xi + h, yi + k1h + k4h )
4 16 16
3 2 12 12 8
k6 = f ( xi + h, yi − k1h + k2h + k3h − k4h + k5h )
7 7 7 7 7
yi +1 = yi + (7k1 + 32k3 + 12k4 + 32k5 + 7k6 )
h
90
Cálculo Numérico II-IF392
Solving Systems of ODEs

Cálculo Numérico II-IF392

Learning Objectives:

• Convert a single (or a system of) high order

ODE to a system of first order ODEs.

• Use the methods discussed earlier in this

topic to solve systems of first order ODEs.

Cálculo Numérico II-IF392

Solving a System of First Order
ODEs
• Methods discussed earlier such as Euler, Runge-
Kutta,… are used to solve first order ordinary
differential equations.
• The same formulas will be used to solve a system
of first order ODEs.
• In this case, the differential equation is a vector
equation and the dependent variable is a vector
variable.

Cálculo Numérico II-IF392

Euler Method for Solving a System of First
Order ODEs
Recall Euler method for solving a first order ODE:

dy ( x)
Given = f ( y, x), y (a) = ya
dx

Euler Method :
y (a + h) = y (a) + h f ( y (a), a)
y (a + 2h) = y (a + h) + h f ( y (a + h), a + h)
y (a + 3h) = y (a + 2h) + h f ( y (a + 2h), a + 2h)

Cálculo Numérico II-IF392

Example - Euler Method
Euler method to solve a system of n first order ODEs.
 f1 (Y , x)   y1 (a) 
 f (Y , x)   y (a)
d Y ( x)
Given = F (Y , x) =  2 , Y (a ) =  2 
dx  ...   ... 
   
 f n (Y , x)  yn (a )
Euler Method :
Y (a + h) = Y (a ) + h F (Y (a ), a )
Y (a + 2h) = Y (a + h) + h F (Y (a + h), a + h)
Y (a + 3h) = Y (a + 2h) + h F (Y (a + 2h), a + 2h)

Cálculo Numérico II-IF392

Solving a System of n First Order ODEs
 y1 ( x) 
• Exactly the same formula  y ( x)
is used but the scalar Y ( x) =  2  Y is n 1 vector
 ... 
variables and functions are  
replaced by vector  n 
y ( x )

variables and vector values  d y1 

 dx   f1 (Y , x) 
functions. d y   
d Y ( x)  2   f 2 (Y , x) 
• Y is a vector of length n. = dx = = F (Y , x)
dx    ... 
 ...   
• F(Y,x) is a vector valued  n  n
d y f (Y , x ) 
function.  dx 

Y (a + h) = Y (a ) + h F (Y (a ), a )
Y (a + 2h) = Y (a + h) + h F (Y (a + h), a + h)
Y (a + 3h) = Y (a + 2h) + h F (Y (a + 2h), a + 2h)
Cálculo Numérico II-IF392
Example :
Euler method for solving a system of first order ODEs.
𝑦ሶ1 (𝑥) 𝑦2 𝑦1 (0) −1
= 1 − 𝑦 = 𝐹 𝑌, 𝑥 , 𝑌(0) = =
𝑦ሶ 2 (𝑥) 1 𝑦2 (0) 1

𝑇𝑤𝑜 𝑠𝑡𝑒𝑝𝑠 𝑜𝑓 𝐸𝑢𝑙𝑒𝑟 𝑀𝑒𝑡ℎ𝑜𝑑 𝑤𝑖𝑡ℎ ℎ = 0.1

𝑆𝑇𝐸𝑃 1:
𝑌 0 + ℎ = 𝑌 0 + ℎ𝐹 𝑌 0 , 0
𝑦1 (0.1) 𝑦1 (0) 𝑦2 (0) −1 + 0.1
= + 0.1 =
𝑦2 (0.1) 𝑦2 (0) 1 − 𝑦1 (0) 1 + 0.1(1 + 1)
−0.9
=
1.2
𝑆𝑇𝐸𝑃 2:
𝑌(0 + 2ℎ) = 𝑌(ℎ) + ℎ𝐹(𝑌(ℎ), ℎ)
𝑦1 (0.2) 𝑦1 (0.1) 𝑦2 (0.1) −0.9 + 0.12
= + 0.1 =
𝑦2 (0.2) 𝑦2 (0.1) 1 − 𝑦1 (0.1) 1.2 + 0.1(1 + 0.9)
−0.78
=
1.39
Cálculo Numérico II-IF392
Exercise
Implement a Matlab program to solve the following IVP,
using the Euler method:

Compare the results to the actual solutions:

Cálculo Numérico II-IF392

Example :
RK2 method for solving a system of first order ODEs
𝑦ሶ1 (𝑥) 𝑦2 𝑦1 (0) −1
= 1 − 𝑦 = 𝐹(𝑌, 𝑥), 𝑌(0) = =
𝑦ሶ 2 (𝑥) 1 𝑦2 (0) 1

𝑇𝑤𝑜 𝑠𝑡𝑒𝑝𝑠 𝑜𝑓second order Runge − Kutta Method 𝑤𝑖𝑡ℎ ℎ = 0.1

𝑆𝑇𝐸𝑃 1:
𝑦2 (0) 0.1
𝐾1 = ℎ𝐹(𝑌(0), 0) = 0.1 =
1 − 𝑦1 (0) 0.2

𝑦2 0 + 0.1(1 − 𝑦1 0 ) 0.12
𝐾2 = ℎ𝐹 𝑌 0 + 𝐾1,0 + ℎ = 0.1 =
1 − 𝑦1 0 + 0.1𝑦2 0 0.19

𝑌 0 + ℎ = 𝑌 0 + 0.5 𝐾1 + 𝐾2
𝑦1 (0.1) −1 1 0.1 0.12 −0.89
= + + =
𝑦2 (0.1) 1 2 0.2 0.19 1.195

Cálculo Numérico II-IF392

Example :
RK2 method for solving a system of first order ODEs

𝑦ሶ 1 (𝑥) 𝑦2 𝑦 (0) −1
= 1 − 𝑦 = 𝐹(𝑌, 𝑥), 𝑌(0) = 1 =
𝑦ሶ 2 (𝑥) 1 𝑦2 (0) 1

𝑆𝑇𝐸𝑃 2:
𝑦2 (0.1) 0.1195
𝐾1 = ℎ𝐹(𝑌(0.1), 0.1) = 0.1 =
1 − 𝑦1 (0.1) 0.1890

𝑦2 0.1 + 0.1(1 − 𝑦1 0.1 ) 0.1384

𝐾2 = ℎ𝐹 𝑌 0.1 + 𝐾1,0.1 + ℎ = 0.1 =
1 − 𝑦1 0.1 + 0.1𝑦2 0.1 0.1771

𝑌 0.1 + ℎ = 𝑌 0.1 + 0.5 𝐾1 + 𝐾2

𝑦1 (0.2) −0.89 1 0.1195 0.1384 −0.7611

= + + =
𝑦2 (0.2) 1.195 2 0.1890 0.1771 1.3780

Cálculo Numérico II-IF392

 Example :
RK4 method for solving a system of first order ODEs
𝑦ሶ 1 (𝑥) 𝑦2 𝑦1 (0) −1
= 1 − 𝑦 = 𝐹(𝑌, 𝑥), 𝑌(0) = =
𝑦ሶ 2 (𝑥) 1 𝑦2 (0) 1
One step of fourth order Runge-Kutta Method with
h=0.1
𝑘1 = 𝑓(𝑥𝑖 , 𝑦𝑖 )
Remember: ℎ 1
𝑘2 = 𝑓(𝑥𝑖 + , 𝑦𝑖 + 𝑘1 ℎ)
2 2
ℎ 1
𝑘3 = 𝑓(𝑥𝑖 + , 𝑦i + 𝑘2 ℎ)
2 2
𝑘4 = 𝑓(𝑥𝑖 + ℎ, 𝑦i + 𝑘3 ℎ)
ℎ
𝑦𝑖+1 = 𝑦𝑖 + 𝑘1 + 2𝑘2 + 2𝑘3 + 𝑘4
6

Cálculo Numérico II-IF392

 Example :
RK4 method for solving a system of first order ODEs
𝑦ሶ 1 (𝑥) 𝑦2 𝑦1 (0) −1
= 1 − 𝑦 = 𝐹(𝑌, 𝑥), 𝑌(0) = =
𝑦ሶ 2 (𝑥) 1 𝑦2 (0) 1
Step 1:
𝑦2 (0) 1 𝟏
𝑲𝟏 = 𝐹 𝑌 0 , 0 = = =
1 − 𝑦1 (0) 1 − (−1) 𝟐

ℎ 0 0.05 0.05
𝐾2 =? ? ? ? 𝑋 0 + = + =
2 0 0.05 0.05
𝐾1ℎ −1 1 −0.95
𝑌 0 + = + 0.05 =
2 1 2 1.1

𝑦2 (0.05) 1.1 𝟏. 𝟏
𝑲𝟐 = 𝐹 𝑌 0 + 𝑘1ℎ/2,0.05 = = =
1 − 𝑦1 (0.05) 1 − (−0.95) 𝟏. 𝟗𝟓

Cálculo Numérico II-IF392

 Example :
RK4 method for solving a system of first order ODEs
𝑦ሶ 1 (𝑥) 𝑦2 𝑦1 (0) −1
= 1 − 𝑦 = 𝐹(𝑌, 𝑥), 𝑌(0) = =
𝑦ሶ 2 (𝑥) 1 𝑦2 (0) 1
Step 1:
𝐾3 =? ? ? ? ℎ 0 0.05 0.05
𝑋 0 + = + =
2 0 0.05 0.05

𝐾2ℎ −1 1.1 −0.945

𝑌 0 + = + 0.05 =
2 1 1.95 1.0975

𝑦2 (0.05) 1.0975 𝟏. 𝟎𝟗𝟕𝟓

𝑲𝟑 = 𝐹 𝑌 0 + 𝑘2ℎ/2,0.05 = = =
1 − 𝑦1 (0.05) 1 − (−0.945) 𝟏. 𝟗𝟒𝟓

Cálculo Numérico II-IF392

 Example :
RK4 method for solving a system of first order ODEs
𝑦ሶ 1 (𝑥) 𝑦2 𝑦1 (0) −1
= 1 − 𝑦 = 𝐹(𝑌, 𝑥), 𝑌(0) = =
𝑦ሶ 2 (𝑥) 1 𝑦2 (0) 1
Step 1:
𝐾4 =? ? ? ? 0 0.1 0.1
𝑋 0 +ℎ = + =
0 0.1 0.1

−1 1.0975 −0.89025
𝑌 0 + ℎ𝐾3 = + 0.1 =
1 1.945 1.1945

𝑦2 (0.1) 1.1945 𝟏. 𝟏𝟗𝟒𝟓

𝑲𝟒 = 𝐹 𝑌 0 + 𝑘3ℎ, 0.1 = = =
1 − 𝑦1 (0.1) 1 − (−0.89025) 𝟏. 𝟖𝟗𝟎𝟐𝟓

Cálculo Numérico II-IF392

 Example :
RK4 method for solving a system of first order ODEs
𝑦ሶ 1 (𝑥) 𝑦2 𝑦1 (0) −1
= 1 − 𝑦 = 𝐹(𝑌, 𝑥), 𝑌(0) = =
𝑦ሶ 2 (𝑥) 1 𝑦2 (0) 1
Step 1:
ℎ
𝑌 0.1 = 𝑌 0 + ( )(𝑘1 + 2𝑘2 + 2𝑘3 + 𝑘4 )
6

−1 0.1 1 1.1 1.0975 1.1945

𝑌 0.1 = + ( +2 +2 + )
1 6 2 1.95 1.945 1.89025

−𝟎. 𝟖𝟗𝟎𝟏𝟕𝟓
𝒀 𝟎. 𝟏 =
𝟏. 𝟏𝟗𝟒𝟔𝟕𝟏

Cálculo Numérico II-IF392

Methods for Solving a System of First Order ODEs

• We have extended Euler, RK2 and RK4 methods to solve

systems of first order ODEs.

• Other methods used to solve first order ODE can be easily

extended to solve systems of first order ODEs.

Cálculo Numérico II-IF392

High Order ODEs
• How do solve a second order ODE?

x + 3 x + 6 x = 1
• How do solve high order ODEs?

Cálculo Numérico II-IF392

The General Approach to Solve ODEs

High order ODE Convert System of first order ODEs Solve

 z1   z2 
 z  = 1 − 3 z − 6 z ,
x + 3x + 6 x = 1 Convert
 2  2 1
Solve
x (0) = 1; x(0) = 4 4
Z ( 0) =  
1 
Second order ODE Two first order ODEs

Cálculo Numérico II-IF392

Conversion Procedure
High order ODE Convert System of first order ODEs Solve

1. Select the dependent variables

One way is to take the original dependent variable and its
derivatives up to one degree less than the highest order
derivative.
2. Write the Differential Equations in terms of the new
variables. The equations come from the way the new
variables are defined or from the original equation.
3. Express the equations in a matrix form.

Cálculo Numérico II-IF392

Example of Converting a High Order
ODE to First Order ODEs
Convert x + 3 x + 6 x = 1, x (0) = 1; x(0) = 4
to a system of first order ODEs

1. Select a new set of variables

(Second order ODE  We need two variables)
z1 = x
One degree less than the highest
z 2 = x order derivative

Cálculo Numérico II-IF392

Example of Converting a High Order
ODE to First Order ODEs
old new Initial Equation
name name cond.
x z1 4 z1 = z 2
x z2 1 z2 = 1 − 3 z 2 − 6 z1

 z1   z2   4
 z  = 1 − 3 z − 6 z , Z (0) = 1 
 2  2 1  
Cálculo Numérico II-IF392
Example of Converting a High Order
ODE to First Order ODEs
Convert
x + 2 x + 7 x + 8 x = 0
x(0) = 9, x (0) = 1; x(0) = 4

1. Select a new set of variables (3 of them)

z1 = x
z 2 = x One degree less than the
z3 = x highest order derivative
Cálculo Numérico II-IF392
Example of Converting a High Order
ODE to First Order ODEs
old new Initial Equation
name name cond.
x z1 4 z1 = z 2
x z2 1 z2 = z3
x z3 9 z3 = −2 z3 − 7 z 2 − 8 z1
 z1   z2   4
 z  =  z  , Z ( 0) =  1 
 2  3   
 z3  − 2 z3 − 7 z 2 − 8 z1  9 
Cálculo Numérico II-IF392
 Exercise: Consider the third-order IVP

Using Euler’s method for systems, with step size h = 0.1, find an estimate for y2 = y(0.2).

Solution:

Cálculo Numérico II-IF392

Exercise: Consider the third-order IVP

Cálculo Numérico II-IF392

Exercise: Consider the third-order IVP

Cálculo Numérico II-IF392

Conversion Procedure for Systems of High
Order ODEs
System of high order ODEs Convert System of first order ODE Solve

1. Select the dependent variables

Take the original dependent variables and their derivatives up
to one degree less than the highest order derivative for each
variable.
2. Write the Differential Equations in terms of the new
variables. The equations come from the way the new
variables are defined or from the original equation.
3. Express the equations in a matrix form.

Cálculo Numérico II-IF392

Example of Converting a High Order
ODE to First Order ODEs
Convert
x + 5x + 2 x + 8 y = 0
y + 2 xy + x = 2
x(0) = 4; x (0) = 2; x(0) = 9; y (0) = 1; y (0) = −3
1. Select a new set of variables ((3 + 2) variables)
z1 = x
z 2 = x One degree less than
the highest order
z3 = x
derivative
z4 = y
z5 = y
Cálculo Numérico II-IF392
 Example of Converting a High Order
ODE to First Order ODEs
old new Initial Equation
name name cond.
x z1 4 z1 = z 2
x z2 2 z2 = z3
x z3 9 z3 = −5 z3 − 2 z 2 − 8 z 4
y z4 1 z4 = z5
y z5 −3 z5 = 2 − z 2 − 2 z1 z 4
Cálculo Numérico II-IF392
 Solution of a Second Order ODE
• Solve the equation using Euler method. Use h=0.1

x + 2 x + 8 x = 2
x(0) = 1; x (0) = −2
Select a new set of variables : z1 = x, z 2 = x
The second order equation is expressed as :

  1  
z z  1
Z = F (Z ) =   =  2
 , Z ( 0) =  
 z2  2 − 2 z 2 − 8 z1   − 2
Cálculo Numérico II-IF392
Solution of a Second Order ODE
 z2  1
F (Z ) =   , Z (0) =  , h = 0.1
2 − 2 z 2 − 8 z1   − 2
Z (0 + 0.1) = Z (0) + hF ( Z (0))
1  −2   0.8 
=   + 0.1  = 
 − 2 2 − 2(−2) − 8(1) − 2.2
Z (0.2) = Z (0.1) + hF ( Z (0.1))
 0.8   − 2.2   0.58 
=  + 0.1  = 
 − 2 .2   2 − 2( −2 . 2) − 8( 0.8)   − 2 . 2 
Cálculo Numérico II-IF392
Summary
• Formulas used in solving a first order ODE are used
to solve systems of first order ODEs.
• Instead of scalar variables and functions, we have vector
variables and vector functions.

• High order ODEs are converted to a set of first

order ODEs.

Cálculo Numérico II-IF392

Example
Solve this exercise using the Runge-Kutta method of order 4 for this
system with ℎ = 0.1 in the interval 0 ≤ 𝑡 ≤ 1 𝑠

Kirchhoff's circuit laws

Cálculo Numérico II-IF392

Example
Solve this exercise using the Runge-Kutta method of order 4 for this
system with ℎ = 0.01 in the interval 0 ≤ 𝑡 ≤ 5 𝑠

Cálculo Numérico II-IF392