Runge-Kutta Methods
CS/SE 4X03
Ned Nedialkov McMaster University March 24, 2021
Outline
Trapezoid
Implicit trapezoidal method Explicit trapezoidal method
Midpoint
Implicit midpoint method Explicit midpoint method
4th order Runge-Kutta Stepsize control
Trapezoid Midpoint 4th order Runge-Kutta Stepsize control Implicit trapezoidal method
• Consider y′(t) = f(t,y), y(ti) = yi
• Fromy(ti+1)=y(ti)+t f(s,y(s))ds,
f(s, y(s))ds
y(ti+1) = y(ti) +
ti
• Use the trapezoidal rule for the integral
ti
ti+1
ti+1 ti
y(ti+1) = y(ti) +
≈ y(ti) + h[f(ti, y(ti)) + f(ti+1, y(ti+1))]
f(s, y(s))ds
2
Copyright © 2021 N. Nedialkov. All rights reserved. 3/12
Trapezoid Midpoint 4th order Runge-Kutta Stepsize control
• From
y(ti+1) ≈ y(ti) + h[f(ti, y(ti)) + f(ti+1, y(ti+1))]
write
2
yi+1 =yi+h[f(ti,yi)+f(ti+1,yi+1)] 2
This is the implicit trapezoidal method
• We have to solve a nonlinear system in general for yi+1
Copyright © 2021 N. Nedialkov. All rights reserved. 4/12
Trapezoid Midpoint 4th order Runge-Kutta Stepsize control
• Local truncation error is
di = y(ti+1)−y(ti) − 1[f(ti,y(ti))+f(ti+1,y(ti+1))]
• di = O(h2)
h2
Copyright © 2021 N. Nedialkov. All rights reserved. 5/12
Trapezoid Midpoint 4th order Runge-Kutta Stepsize control Explicit trapezoidal method
• In the implicit trapezoidal rule, we need to solve for yi+1
• We can approximate y(ti+1) first using forward Euler:
Y =yi+hf(ti,yi)
• Then plug Y into the formula for the implicit trapezoidal
method
• This is a two-stage explicit Runge-Kutta method
yi+1 =yi+h[f(ti,yi)+f(ti+1,Y)] 2
• Local truncation error is
di = y(ti+1)−y(ti) − 1[f(ti,y(ti))+f(ti+1,y(ti)+hf(ti,y(ti)))]
h2
di = O(h2), a bit involved to derive it
Copyright © 2021 N. Nedialkov. All rights reserved. 6/12
Trapezoid Midpoint 4th order Runge-Kutta Stepsize control Implicit midpoint
• Use the midpoint quadrature rule:
yi+1 = yi + hf ti+1/2, yi+1/2
=yi +hf(ti +h/2,(yi +yi+1)/2)
• That is, we solve for yi+1 • Order is 2
Copyright © 2021 N. Nedialkov. All rights reserved. 7/12
Trapezoid Midpoint 4th order Runge-Kutta Stepsize control Explicit midpoint method
• Take a step of size h/2 with forward Euler Y =yi+hf(ti,yi)
2
• Plug into the formula from the midpoint quadrature rule: yi+1 =yi +hf(ti +h/2,Y),
• This is a two-stage explicit Runge-Kutta method • Order is 2
Copyright © 2021 N. Nedialkov. All rights reserved. 8/12
Trapezoid Midpoint 4th order Runge-Kutta Stepsize control Classical 4th order Runge-Kutta
• Based on Simpson’s quadrature rule • 4 stages
• Order 4, O(h4) accuracy
Y1 = yi
Y2 = yi + hf(ti,Y1)
2
Y3 =yi+hf(ti+h/2,Y2) 2
Y4 =yi +hf(ti +h/2,Y3)
yi+1 =yi+hf(ti,Y1)+2f(ti+h/2,Y2)+2f(ti+h/2,Y3)+f(ti+1,Y4)
6
Copyright © 2021 N. Nedialkov. All rights reserved. 9/12
Trapezoid Midpoint 4th order Runge-Kutta Stepsize control Stepsize control
• Estimate the error: Runge-Kutta pair (details omitted)
• Let hi be the current stepsize
• Local error is of the form ei = chq+1 i
• Assume an estimate for ei is computed
• We require ei = chq+1 ≤ tol
• If the tolerance is satisfied, we accept the step and predict stepsize for the next step
• Otherwise, reject the step and repeat with smaller h ̄i
i
Copyright © 2021 N. Nedialkov. All rights reserved.
10/12
Trapezoid Midpoint 4th order Runge-Kutta Stepsize control
• From chq+1 = ei, i
and
From which
c= ei hq+1
hi+1 = hi ei • Since tol ≥ ei, hi+1 ≥ hi
chq+1 = i+1
i
ei hq+1 = tol hq+1 i+1
i
tol 1/(q+1)
Copyright © 2021 N. Nedialkov. All rights reserved.
11/12
Trapezoid Midpoint 4th order Runge-Kutta Stepsize control
• If ei > tol , the stepsize is rejected • Repeat the step with
̄ tol 1/(q+1) hi = hi ei
• For “safety”, typically new stepsize is computed by 0.5 tol 1/(q+1)
0.9hi ei
Copyright © 2021 N. Nedialkov. All rights reserved.
12/12