Solving differential equations
Jiří Chmelík, Marek Trtík
PA199
2
 Initial value problem for ordinary differential equations.
 Forward Euler’s method.
 Backward Euler’s method.
 Midpoint method.
 Runge-Kutta methods.
Outline
3
Initial value problem
4
 Initial value problem (IVP) for the 1st order ordinary differential
equations (ODE)s:
ሶ𝒚 = 𝑭 𝒚, 𝑡 , 𝒚 𝑡0 = 𝒚0
 𝒚 𝑡 = 𝑦1 𝑡 , … , 𝑦𝑛 𝑡
⊤
is a vector of unknown functions 𝑦𝑖: ℝ → ℝ.
 𝑭 𝒚, 𝑡 = 𝑓1 𝒚 𝑡 , 𝑡 , … , 𝑓𝑛 𝒚 𝑡 , 𝑡
⊤
is a vector of known fns 𝑓𝑖: ℝ 𝑛+1 → ℝ.
 The initial value condition:
 𝑡0 given time point.
 𝒚0 = 𝑦1 𝑡0 , … , 𝑦𝑛 𝑡0
⊤
is a vector of known values of functions 𝑦𝑖 at 𝑡0.
 Solution: Any vector of functions ෝ𝒚(𝑡) = (ො𝑦1 𝑡 , … , ො𝑦𝑛 𝑡 ) s.t.:
ሶෝ𝒚 = 𝑭 ෝ𝒚, 𝑡 , ෝ𝒚 𝑡0 = 𝒚0
 NOTE: We can extend to higher orders, e.g., ሷ𝒚 = 𝑭(𝒚, ሶ𝒚, 𝑡), 𝒚 𝑡0 = 𝒚0.
 We can also have initial condition for derivatives, e.g., ሶ𝒚 𝑡0 = ሶ𝒚0.
Initial value problem
5
 Example: Check that 𝑦 𝑡 =
3
4
+
𝑐
𝑡2 , 𝑐 ∈ ℝ is a general solution to
ሶ𝑦 =
3−4𝑦
2𝑡
. Find 𝑐 for which initial condition 𝑦 1 = −4 is satisfied.
 Solution:

𝑑
𝑑𝑡
3
4
+
𝑐
𝑡2 = −
2𝑐
𝑡3

3−4(
3
4
+
𝑐
𝑡2)
2𝑡
= −
4𝑐
𝑡2 ⋅
1
2𝑡
= −
2𝑐
𝑡3

3
4
+
𝑐
12 = −4 => 𝑐 = −
19
4
 In physics simulations:
 Initial conditions define current state of the system.
Initial value problem
6
 Plot of a function 𝑭(𝒚, 𝑡) for some values of 𝒚, 𝑡.
 Goal: Get visual impression about derivatives 𝒚.
 Example: Show direction field for ሶ𝑦 = 𝑦 − 𝑡. [axes: 𝑡 horizontal, 𝑦 vertical]
Direction field
Step 1
Find contours:
ሶ𝑦 = 𝑦 − 𝑡 = 𝑐, 𝑐 ∈ R.
Step 2
Draw 𝑦’s slopes
using arrows.
Step 3
Draw more arrows.
Step 4
Predict solutions.
Picturessource:[3]
7
 The goal is find 𝒚 𝑡1 , where 𝑡1 > 𝑡0, for a given IVP ሶ𝒚 = 𝑭 𝒚, 𝑡 , 𝒚 𝑡0 = 𝒚0:
 Start at the initial time 𝑡0 and the initial value 𝒚0.
 Compute a sequence of values 𝒚 𝑡0 + Δ𝑡 , 𝒚 𝑡0 + 2Δ𝑡 , … , 𝒚 𝑡0 + 𝑛Δ𝑡 , where
𝑡1 = 𝑡0 + 𝑛Δ𝑡.
 There are two kinds of methods:
 Explicit methods:
 Compute 𝒚 𝑡0 + Δ𝑡 by a function 𝓕(𝑭, 𝒚0, 𝑡0, Δ𝑡) of current state of the system.
 Implicit methods:
 Compute 𝒚 𝑡0 + Δ𝑡 by a solution of an equation 𝓕 𝑭, 𝒚0, 𝑡0, Δ𝑡, 𝒚 𝑡0 + Δ𝑡 = 0
over the current and future state of the system.
Numerical solution
8
 For a 𝑘-times differentiable function 𝑦: ℝ → ℝ at a point 𝑡0 ∈ 𝐷(𝑦) there
exists a polynomial 𝑃𝑘: ℝ → ℝ and a functions 𝑅: ℝ → ℝ s.t.
 𝑦 𝑡0 + Δ𝑡 = 𝑃𝑘 𝑡0 + Δ𝑡 + Δ𝑡 𝑘
𝑅(𝑡0 + Δ𝑡)
 𝑃𝑘 𝑡0 + Δ𝑡 = 𝑦 𝑡0 + Δ𝑡 ሶ𝑦 𝑡0 +
Δ𝑡2
2!
ሷ𝑦 𝑡0 + ⋯ +
Δ𝑡 𝑘
𝑘!
𝑦 𝑘
𝑡0
 lim
Δ𝑡→0
𝑅 𝑡0 + Δ𝑡 = 0
Taylor theorem
9
 Examples (Taylor approximation):
Numerical solution
𝑦 𝑡 = 𝑒 𝑡
, 𝑃1 𝑡 = 1 + 𝑡, 𝑡0 = 0. 𝑦 𝑡 = 𝑒 𝑡
, 𝑃2 𝑡 = 1 + 𝑡 +
𝑡2
2
, 𝑡0 = 0.
10
 What is the error from the approximation using 𝑃𝑘:
𝑦 𝑡0 + Δ𝑡 ≈ 𝑃𝑘 𝑡0 + Δ𝑡
 It is a distance from the exact value 𝑃𝑘+1 𝑡0 + Δ𝑡 + Δ𝑡 𝑘+1 𝑅(𝑡0 + Δ𝑡):
error = 𝑃𝑘+1 𝑡0 + Δ𝑡 + Δ𝑡 𝑘+1 𝑅 𝑡0 + Δ𝑡 − 𝑃𝑘 𝑡0 + Δ𝑡
=
Δ𝑡 𝑘+1
(𝑘 + 1)!
𝑦 𝑘+1 𝑡0 + Δ𝑡 𝑘+1 𝑅(𝑡0 + Δ𝑡)
 For small Δ𝑡 the error is proportional to the term Δ𝑡 𝑘+1. Therefore,
𝑦 𝑡0 + Δ𝑡 = 𝑃𝑘 𝑡0 + Δ𝑡 + 𝒪(Δ𝑡 𝑘+1)
“𝒪” error notation
11
 We get forward Euler’s method, when we approximate 𝑦 by 𝑃1 at 𝑡0:
𝑦 𝑡0 + Δ𝑡 ≈ 𝑦 𝑡0 + Δ𝑡 ሶ𝑦 𝑡0
= 𝑦 𝑡0 + Δ𝑡𝐹 𝑦(𝑡0), 𝑡0
 We see that forward Euler’s method is an explicit method.
Forward Euler’s method
12
 Example: Let IPV be ሶ𝑦 =
3−4𝑦
2𝑡
, 𝑦 1 = −4. Compute 𝑦 2 by forward
Euler’s method. [Note: Exact solution is 𝑦 𝑡 =
3
4
−
19
4𝑡2]
 Solution:
 Let’s choose a time step Δ𝑡 =
1
2
=> We must apply the method 2 times.
 𝑦 1 = −4  from the initial condition.
 𝑦
3
2
= 𝑦 1 +
1
2
≈ 𝑦 1 +
1
2
𝐹 𝑦 1 , 1 = −4 +
1
2
⋅
3−4 −4
2⋅1
=
3
4
 1st iteration
 𝑦 2 =
3
4
+
1
2
⋅
3−4⋅
3
4
2⋅
3
2
=
3
4
 2nd iteration
 We see the method is simple and fast.
Forward Euler’s method
13
 Low accuracy issue:
 𝒪 Δ𝑡2 error in each iteration.
 Example:
Forward Euler’s method
IVP: ሶ𝑦 =
3−4𝑦
2𝑡
, 𝑦 1 = −4.
Euler’s method: 𝛥𝑡 =
1
2
, 𝛥𝑡 =
1
4
.
Exact solution: 𝑦 𝑡 =
3
4
−
19
4𝑡2
14
 Instability issue:
 The iteration process may diverge.
 Example:
Forward Euler’s method
IVP: ሶ𝑦 = −2.3𝑦, 𝑦 0 = 1.
Euler’s method: 𝛥𝑡 = 1, 𝛥𝑡 =
1
2
.
Exact solution: 𝑦 𝑡 = 𝑒−2.3𝑡
.
15
 What can we do with the issues?
 Use smaller time step Δ𝑡 to reduce the error and/or avoid the instability.
 But we then need more iterations => slower simulation.
 Choose more accurate/stable solver.
 Suggestion for seminar: Implement method “ODE_Euler_forward”.
void ODE _Euler_forward(
std::vector<float> const& y0, // 𝒙, 𝒗 of particle(s)
std::vector<F_y_t> const& Fyt, // ሶ𝒙, ሶ𝒗 of particle(s), i.e. 𝒗, 𝑭/𝑚
float& t, // current time (to be updated)
float const dt, // time step
std::vector<float>& y) // integrated 𝒙, 𝒗 of particle(s)
{ TODO }
Forward Euler’s method
16
 Example: Let’s consider a particle 𝒫(𝑡) = 𝒙, 𝒗, 𝑭, 𝑚 , where 𝑚 = 0.1kg, in a
homogenous gravity field with 𝒈 = 0,0, −10 ⊤
m ⋅ s−2
. At time 𝑡 = 1 we have
𝒙 = 1, −1,5 ⊤m, 𝒗 = 1,0,0 ⊤m ⋅ s−1. Using forward Euler’s method with Δ𝑡 =
0.5s compute 𝒫 2 .
 Solution: Particle moves by Newton’s equations of motion:
ሶ𝒙 𝑡 = 𝒗 𝑡 , ሶ𝒗(𝑡) =
𝑭
𝑚
Therefore: 𝒙 1.5 =
1 + 0.5 ⋅ 1
−1 + 0.5 ⋅ 0
5 + 0.5 ⋅ 0
=
1.5
−1
5
, 𝒗 1.5 =
1 + 0.5
0.1⋅0
0.1
0 + 0.5
0.1⋅0
0.1
0 + 0.5
0.1⋅−10
0.1
=
1
0
−5
.
𝒙 2 = 2, −1,0 ⊤
, 𝒗 2 = 1,0, −10 ⊤
.
Forward Euler’s method
17
Backward Euler’s method
18
 From the fundamental theorem of the calculus:
න
𝑡0
𝑡0+Δ𝑡
ሶ𝑦 𝑡 𝑑𝑡 = 𝑦 𝑡0 + Δ𝑡 − 𝑦 𝑡0 .
 Therefore,
𝑦 𝑡0 + Δ𝑡 = 𝑦 𝑡0 + න
𝑡0
𝑡0+Δ𝑡
𝐹 𝑦 𝑡 , 𝑡 𝑑𝑡 .
 We can approximate the integral by
“right-hand” rectangle:
න
𝑡0
𝑡0+Δ𝑡
𝐹 𝑦 𝑡 , 𝑡 𝑑𝑡 ≈ Δ𝑡𝐹 𝑦 𝑡0 + Δ𝑡 , 𝑡0 + Δ𝑡 .
Backward Euler’s method
𝑡0 + Δ𝑡𝑡0
Δ𝑡
𝑡
𝐹(𝑦(𝑡), 𝑡)
𝐹(𝑦(𝑡0 + Δ𝑡), 𝑡0 + Δ𝑡)
න
𝑡0
𝑡0+Δ𝑡
𝐹 𝑦 𝑡 , 𝑡 𝑑𝑡
𝐹(𝑦(𝑡0), 𝑡0)
Δ𝑡𝐹 𝑦 𝑡0 + Δ𝑡 , 𝑡0 + Δ𝑡
19
 Backward Euler’s method leads to this equation:
𝑦 𝑡0 + Δ𝑡 ≈ 𝑦 𝑡0 + Δ𝑡𝐹 𝑦 𝑡0 + Δ𝑡 , 𝑡0 + Δ𝑡
 Backward Euler’s method is an implicit method.
 We must solve this equation to obtain the unknown 𝑦 𝑡0 + Δ𝑡 :
𝑦 𝑡0 + Δ𝑡 − 𝑦 𝑡0 − Δ𝑡𝐹 𝑦 𝑡0 + Δ𝑡 , 𝑡0 + Δ𝑡 = 0
 Use any available method for solving the equation, e.g., Newton’s
method (𝑦[𝑘+1]
= 𝑦[𝑘]
− 𝓕(𝑦[𝑘]
)/ ሶ𝓕(𝑦[𝑘]
)).
 Note: If we have system of ODEs, then we get system of equations.
Backward Euler’s method
20
 Example: Let IPV be ሶ𝑦 =
3−4𝑦
2𝑡
, 𝑦 1 = −4. Compute 𝑦 2 by backward
Euler’s method. [Note: Exact solution is 𝑦 𝑡 =
3
4
−
19
4𝑡2]
 Solution:
 Let’s choose a time step Δ𝑡 =
1
2
=> We must apply the method 2 times.
 𝑦 1 = −4  from the initial condition.
 𝑦
3
2
= 𝑦 1 +
1
2
= 𝑦 1 +
1
2
𝐹 𝑦
3
2
,
3
2
= −4 +
1
2
⋅
3−4𝑦
3
2
2⋅
3
2
= −
7
2
−
2
3
𝑦
3
2
 We solve: 𝑦
3
2
= −
7
2
−
2
3
𝑦
3
2
=> 𝑦
3
2
= −
7
2 1+
2
3
= −
21
10
 1st iteration
 𝑦 2 = −
21
10
+
1
2
⋅
3−4𝑦 2
2⋅2
= −
21
10
+
3
8
−
1
2
𝑦(2) => 𝑦 2 = −
23
20
 2nd iteration
Backward Euler’s method
21
 We can plot out result and compare it with forward Euler’s method:
Backward Euler’s method
IVP: ሶ𝑦 =
3−4𝑦
2𝑡
, 𝑦 1 = −4.
Backward Euler: 𝛥𝑡 =
1
2
.
Forward Euler: 𝛥𝑡 =
1
2
.
Exact solution: 𝑦 𝑡 =
3
4
−
19
4𝑡2
22
 Example 2: Let IPV be ሶ𝑦 = −2.3𝑦, 𝑦 0 = 1. Compute 𝑦 4 by backward
Euler’s method.
 Solution:
 Let’s choose a time step Δ𝑡 = 1 => We must apply the method 4 times.
 𝑦 0 = 1  from the initial condition.
 𝑦 1 = 1 + 1 ⋅ −2.3𝑦(1) => 𝑦 1 =
1
1+2.3
=
10
33
 1st iteration
 𝑦 2 =
10
33
+ 1 ⋅ −2.3𝑦(2) => 𝑦 2 =
10
33
⋅
10
33
=
10
33
2
 2nd iteration
 𝑦 3 =
10
33
2
− 2.3𝑦(3) => 𝑦 3 =
10
33
2
⋅
10
33
=
10
33
3
 3rd iteration
 𝑦 4 =
10
33
3
− 2.3𝑦(4) => 𝑦 4 =
10
33
4
 4th iteration
 Note: Observe the geometric progression 𝑦 𝑘+1 = 𝑞𝑦 𝑘, 𝑞 =
10
33
=> 𝑦 𝑘 = 𝑞 𝑘 𝑦0.
Backward Euler’s method
23
 We can plot out result and compare it with forward Euler’s method:
Backward Euler’s method
IVP: ሶ𝑦 = −2.3𝑦, 𝑦 0 = 1.
Backward Euler: 𝛥𝑡 = 1.
Forward Euler: 𝛥𝑡 = 1.
Exact solution: 𝑦 𝑡 = 𝑒−2.3𝑡
24
 Properties of backward Euler’s method
 Hard to implement.
 Requires solving an equation or a system of equations.
 𝒪 Δ𝑡2 error in each iteration.
 Stable for large time step Δ𝑡.
 Choice between forward/backward Euler’s method depends on a
problem. “Rule of thumb”:
 Prefer forward method for “stable” problems.
 Prefer backward method for “stiff” problems.
Backward Euler’s method
25
Midpoint method
26
 Let’s try to approximate 𝑦 by 𝑃2 at 𝑡0:
𝑦 𝑡0 + Δ𝑡 = 𝑦 𝑡0 + Δ𝑡 ሶ𝑦 𝑡0 +
Δ𝑡2
2!
ሷ𝑦 𝑡0 + 𝒪(Δ𝑡3)
= 𝑦 𝑡0 + Δ𝑡𝐹 𝑦(𝑡0), 𝑡0 +
Δ𝑡2
2
ሶ𝐹 𝑦(𝑡0), 𝑡0 + 𝒪(Δ𝑡3)
 How to compute ሶ𝐹?
 Using the chain rule, we get: ሶ𝐹 =
𝜕𝐹
𝜕𝑡
+
𝜕𝐹
𝜕𝑦
ሶ𝑦 =
𝜕𝐹
𝜕𝑡
+
𝜕𝐹
𝜕𝑦
𝐹
 Not much better, because we still do not know
𝜕𝐹
𝜕𝑡
,
𝜕𝐹
𝜕𝑦
.
 So, let’s try to approximate 𝐹 using 𝑃1 …
 Note: We must use a 2-variables version of Taylor’s theorem.
Midpoint method
(*)
27
𝐹 𝑦(𝑡0) + Δ𝑦, 𝑡0 + Δ𝑡 = 𝐹 𝑦 𝑡0 , 𝑡0 + Δ𝑦
𝜕𝐹
𝜕𝑦
𝑦 𝑡0 , 𝑡0 + Δ𝑡
𝜕𝐹
𝜕𝑡
𝑦 𝑡0 , 𝑡0 + 𝒪 Δ𝑦2
+ Δ𝑡2
 Let’s substitute: Δ𝑦 →
Δ𝑡
2
𝐹 𝑦(𝑡0), 𝑡0 , Δ𝑡 →
Δ𝑡
2
𝐹 𝑦(𝑡0) +
Δ𝑡
2
𝐹 𝑦 𝑡0 , 𝑡0 , 𝑡0 +
Δ𝑡
2
= 𝐹 𝑦 𝑡0 , 𝑡0 +
Δ𝑡
2
𝐹 𝑦 𝑡0 , 𝑡0
𝜕𝐹
𝜕𝑦
𝑦 𝑡0 , 𝑡0 +
Δ𝑡
2
𝜕𝐹
𝜕𝑡
𝑦(𝑡0), 𝑡0
+ 𝒪(Δ𝑡2
)
𝐹 𝑦(𝑡0) +
Δ𝑡
2
𝐹 𝑦 𝑡0 , 𝑡0 , 𝑡0 +
Δ𝑡
2
= 𝐹 𝑦 𝑡0 , 𝑡0 +
Δ𝑡
2
ሶ𝐹 𝑦 𝑡0 , 𝑡0 + 𝒪(Δ𝑡2
)
Δ𝑡
2
ሶ𝐹 𝑦 𝑡0 , 𝑡0 + 𝒪 Δ𝑡2
= 𝐹 𝑦 𝑡0 +
Δ𝑡
2
𝐹 𝑦 𝑡0 , 𝑡0 , 𝑡0 +
Δ𝑡
2
− 𝐹 𝑦 𝑡0 , 𝑡0
Δ𝑡2
2
ሶ𝐹 𝑦 𝑡0 , 𝑡0 + 𝒪(Δ𝑡3
) = Δ𝑡𝐹 𝑦(𝑡0) +
Δ𝑡
2
𝐹 𝑦 𝑡0 , 𝑡0 , 𝑡0 +
Δ𝑡
2
− Δ𝑡𝐹 𝑦 𝑡0 , 𝑡0
Midpoint method
Δ𝑡
2
ሶ𝐹 𝑦 𝑡0 , 𝑡0
(*)
28
𝑦 𝑡0 + Δ𝑡 = 𝑦 𝑡0 + Δ𝑡𝐹 𝑦(𝑡0), 𝑡0 + Δ𝑡𝐹 𝑦(𝑡0) +
Δ𝑡
2
𝐹 𝑦 𝑡0 , 𝑡0 , 𝑡0 +
Δ𝑡
2
− Δ𝑡𝐹 𝑦 𝑡0 , 𝑡0
𝑦 𝑡0 + Δ𝑡 = 𝑦 𝑡0 + Δ𝑡𝐹 𝑦(𝑡0) +
Δ𝑡
2
𝐹 𝑦 𝑡0 , 𝑡0 , 𝑡0 +
Δ𝑡
2
 This is an explicit method.
 This method is more accurate than Euler’s method:
 Euler: 𝒪(Δ𝑡2)
 Midpoint: 𝒪(Δ𝑡3)
Midpoint method
29
𝑦 𝑡0 + Δ𝑡 = 𝑦 𝑡0 + Δ𝑡𝐹 𝑦(𝑡0) +
Δ𝑡
2
𝐹 𝑦 𝑡0 , 𝑡0 , 𝑡0 +
Δ𝑡
2
Midpoint method
𝑡0 + Δ𝑡𝑡0
𝑡
𝑦
𝑦(𝑡0)
𝑡0 +
Δ𝑡
2
Slope: 𝐹(𝑦(𝑡0), 𝑡0)
𝑦 𝑡0 + Δ𝑡/2
𝑦(𝑡0 + Δ𝑡)
Slope: 𝐹(𝑦(𝑡0 +
Δ𝑡
2
), 𝑡0 +
Δ𝑡
2
)
2nd Euler step by
Δ𝑡
2
.
30
Midpoint method
IVP: ሶ𝑦 =
3−4𝑦
2𝑡
, 𝑦 1 = −4.
Midpoint method: 𝛥𝑡 =
1
2
.
Forward Euler: 𝛥𝑡 =
1
2
.
Exact solution: 𝑦 𝑡 =
3
4
−
19
4𝑡2
31
Midpoint method
IVP: ሶ𝑦 = −2.3𝑦, 𝑦 0 = 1.
Midpoint method: 𝛥𝑡 = 1.
Forward Euler: 𝛥𝑡 = 1.
Exact solution: 𝑦 𝑡 = 𝑒−2.3𝑡
32
 There is also implicit version of the midpoint method.
 From the fundamental theorem of the calculus:
𝑦 𝑡0 + Δ𝑡 = 𝑦 𝑡0 + න
𝑡0
𝑡0+Δ𝑡
𝐹 𝑦 𝑡 , 𝑡 𝑑𝑡 .
 We can approximate the integral by “midpoint” rectangle:
න
𝑡0
𝑡0+Δ𝑡
𝐹 𝑦 𝑡 , 𝑡 𝑑𝑡 ≈ Δ𝑡𝐹
𝑦 𝑡0 + 𝑦 𝑡0 + Δ𝑡
2
,
𝑡0 + (𝑡0 + Δ𝑡)
2
 Therefore, we get
𝑦 𝑡0 + Δ𝑡 = 𝑦 𝑡0 + Δ𝑡𝐹
𝑦 𝑡0 + 𝑦 𝑡0 + Δ𝑡
2
, 𝑡0 +
Δ𝑡
2
Midpoint method
33
Runge-Kutta methods
34
 In general, we can approximate the integral as follows:
න
𝑡0
𝑡0+Δ𝑡
𝐹 𝑦 𝑡 , 𝑡 𝑑𝑡 ≈ Δ𝑡 ෍
𝑖=1
𝑛
𝑏𝑖 𝐹 𝑦 𝑡0 + 𝑐𝑖Δ𝑡 , 𝑡0 + 𝑐𝑖Δ𝑡
 The problem is that values 𝑦 𝑡0 + 𝑐𝑖Δ𝑡 are unknown!
 Runge-Kutta methods solve the issue by this substitution:
𝑘1 = 𝐹 𝑦 𝑡0 , 𝑡0
𝑘𝑖 = 𝐹 𝑦 𝑡0 + Δ𝑡 ෍
𝑗=1
𝑖−1
𝑎𝑖,𝑗 𝑘𝑗 , 𝑡0 + 𝑐𝑖Δ𝑡 , s. t. ෍
𝑗=1
𝑖−1
𝑎𝑖,𝑗 = 𝑐𝑖
න
𝑡0
𝑡0+Δ𝑡
𝐹 𝑦 𝑡 , 𝑡 𝑑𝑡 ≈ Δ𝑡 ෍
𝑖=1
𝑛
𝑏𝑖 𝑘𝑖
Runge-Kutta methods
35
 Therefore, Runge-Kutta of order 𝑛 is defined as:
𝑦 𝑡0 + Δ𝑡 = 𝑦 𝑡0 + Δ𝑡 ෍
𝑖=1
𝑛
𝑏𝑖 𝑘𝑖 ,
where terms 𝑘𝑖 were defined on the previous slide.
 However, we must compute the numbers 𝑎𝑖,𝑗, 𝑏𝑖, 𝑐𝑖 so that resulting
expression yields an approximation by Taylor’s polynomial 𝑃𝑛.
Runge-Kutta methods
36
 Example: Runge-Kutta method of order 1 (i.e. 𝑛 = 1):
𝑘1 = 𝐹 𝑦 𝑡0 , 𝑡0
𝑦 𝑡0 + Δ𝑡 = 𝑦 𝑡0 + Δ𝑡𝑏1 𝑘1 = 𝑦 𝑡0 + Δ𝑡𝑏1 𝐹 𝑦 𝑡0 , 𝑡0
What value we should choose for 𝑏1? We compare 𝑦 𝑡0 + Δ𝑡 with 𝑃1.
𝑃1 𝑡0 + Δ𝑡 = 𝑦 𝑡0 + Δ𝑡 ሶ𝑦 𝑡0 = 𝑦 𝑡0 + Δ𝑡𝐹 𝑦 𝑡0 , 𝑡0
Therefore, 𝑏1 must be 1.
 Observation: Euler’s method is Runge-Kutta method of order 1.
Runge-Kutta methods
37
 Example: Runge-Kutta method of order 2 (i.e. 𝑛 = 2):
𝑘1 = 𝐹 𝑦 𝑡0 , 𝑡0
𝑘2 = 𝐹 𝑦 𝑡0 + Δ𝑡𝑎2,1 𝑘1, 𝑡0 + 𝑐2Δ𝑡
𝑦 𝑡0 + Δ𝑡 = 𝑦 𝑡0 + Δ𝑡𝑏1 𝑘1 + Δ𝑡𝑏2 𝑘2
We compute 𝑎2,1, 𝑏1, 𝑏2 by comparison of 𝑦 𝑡0 + Δ𝑡 with 𝑃2.
𝑃2 𝑡0 + Δ𝑡 = 𝑦 𝑡0 + Δ𝑡 ሶ𝑦 𝑡0 +
Δ𝑡2
2!
ሷ𝑦 𝑡0
= 𝑦 𝑡0 + Δ𝑡𝐹 𝑦(𝑡0), 𝑡0 +
Δ𝑡2
2
ሶ𝐹 𝑦(𝑡0), 𝑡0
= 𝑦 𝑡0 + Δ𝑡𝐹 𝑦(𝑡0), 𝑡0 +
Δ𝑡2
2
𝜕𝐹
𝜕𝑡
+
𝜕𝐹
𝜕𝑦
𝐹 𝑦(𝑡0), 𝑡0
Runge-Kutta methods
38
For the comparison let’s approximate 𝑘2 by 𝑃1:
𝑘2 = 𝐹 𝑦 𝑡0 + Δ𝑡𝑎2,1 𝑘1, 𝑡0 + 𝑐2Δ𝑡
≈ 𝐹 𝑦 𝑡0 , 𝑡0 + Δ𝑡(𝑐2
𝜕𝐹
𝜕𝑡
+ 𝑎2,1 𝑘1
𝜕𝐹
𝜕𝑦
) 𝑦 𝑡0 , 𝑡0
= 𝐹 𝑦 𝑡0 , 𝑡0 + Δ𝑡(𝑐2
𝜕𝐹
𝜕𝑡
+ 𝑎2,1 𝐹
𝜕𝐹
𝜕𝑦
) 𝑦 𝑡0 , 𝑡0
When we substitute the approximated 𝑘2 we get:
𝑦 𝑡0 + Δ𝑡 = 𝑦 𝑡0 + Δ𝑡𝑏1 𝐹 𝑦 𝑡0 , 𝑡0
+ Δ𝑡𝑏2(𝐹 𝑦 𝑡0 , 𝑡0 + Δ𝑡(𝑐2
𝜕𝐹
𝜕𝑡
+ 𝑎2,1 𝐹
𝜕𝐹
𝜕𝑦
) 𝑦 𝑡0 , 𝑡0 )
= 𝑦 𝑡0 + Δ𝑡 𝑏1 + 𝑏2 𝐹 𝑦 𝑡0 , 𝑡0
+ Δ𝑡2 𝑏2 𝑐2
𝜕𝐹
𝜕𝑡
+ 𝑎2,1
𝜕𝐹
𝜕𝑦
𝐹 𝑦 𝑡0 , 𝑡0
Runge-Kutta methods
39
So, we must solve this system of equations:
𝑏1 + 𝑏2 = 1, 𝑏2 𝑐2 =
1
2
, 𝑏2 𝑎2,1 =
1
2
.
One possible solution is: 𝑏1 = 0, 𝑏2 = 1, 𝑐2 =
1
2
, 𝑎2,1 =
1
2
.
(Note: Another solution is: 𝑏1 =
1
2
, 𝑏2 =
1
2
, 𝑐2 = 1, 𝑎2,1 = 1)
We get the result:
𝑘1 = 𝐹 𝑦 𝑡0 , 𝑡0
𝑘2 = 𝐹 𝑦 𝑡0 +
Δ𝑡
2
𝑘1, 𝑡0 +
Δ𝑡
2
𝑦 𝑡0 + Δ𝑡 = 𝑦 𝑡0 + Δ𝑡𝑘2 = 𝑦 𝑡0 + Δ𝑡𝐹 𝑦 𝑡0 +
Δ𝑡
2
𝐹 𝑦 𝑡0 , 𝑡0 , 𝑡0 +
Δ𝑡
2
.
 Observation: Midpoint method is Runge-Kutta method of order 2.
Runge-Kutta methods
40
 Example: Runge-Kutta method of order 4:
𝑘1 = 𝐹 𝑦(𝑡0), 𝑡0
𝑘2 = 𝐹 𝑦 𝑡0 +
𝑘1
2
, 𝑡0 +
Δ𝑡
2
𝑘3 = 𝐹 𝑦 𝑡0 +
𝑘2
2
, 𝑡0 +
Δ𝑡
2
𝑘4 = 𝐹 𝑦 𝑡0 + 𝑘3, 𝑡0 + Δ𝑡
𝑦 𝑡0 + Δ𝑡 = 𝑦 𝑡0 + Δ𝑡
1
6
𝑘1 +
1
3
𝑘2 +
1
3
𝑘3 +
1
6
𝑘4
Runge-Kutta methods
41
Schema of numerical methods
Picture source: [2]
42
[1] A. Witkin, D. Baraff; Differential Equation Basics; Physically Based
Modeling: Principles and Practice, 1997
[2] J.C.Butcher; Numerical methods for ordinary differential equations;
3rd edition, Wiley, 2016.
[3] https://tutorial.math.lamar.edu/Classes/DE/DE.aspx
References