Dr.

Markus Faustmann
Paul Dunhofer
Institute of Analysis and Scientific Computing
Summer term 2025
due on 14.03.2025

Scientific programming for interdisciplinary mathematics - Worksheet 1

Topics: convergence order, numerical di!erentiation, numerical methods for linear ODEs

Convergence order

We consider an abstract numerical method that approximates a quantity yω → R by {yh }h>0 ↑

R such that
yh ↓ yω for h ↓ 0.
Here h > 0 is a discretization parameter that controls the accuracy of the procedure, i.e.,
smaller values of h lead to better approximations yh . The method is said to have order of
convergence p > 0 (or: convergence rate p > 0) if there exists a constant C > 0 such that

e h ↔ C hp for all h > 0, where eh := |yω ↗ yh |.

With the ansatz eh = C hp , the order of convergence can be expressed by the sequence (ph )h>0
with
log(eh /eh/2 ) log(eh ) ↗ log(eh/2 )
ph = = for h > 0.
log(h/(h/2)) log(2)
Note that the calculation of ph requires the availability of two iterates yh , yh/2 and, in
particular, of the unavailable solution yω . However, one can show that the di!erence e!h :=
|yh/2 ↗ yh | behaves (asymptotically) like the error eh . This motivates the formula for the
experimental convergence order

eh ) ↗ log(!
log(! eh/2 )
p!h = for h > 0,
log(2)

in which the (unknown) error eh is approximated by e!h . Note that the calculation of p!h
requires the availability of three iterates yh , yh/2 , and yh/4 , but not that of the exact solution
yω .
In so-called double logarithmic plots (log-log plots) a logarithmic scaling is used for both
coordinate axes, i.e., log(x) is plotted on the x-axis and log(y) on the y-axis.

Problem 1:
How are power functions of the form eh = C hp represented in a log-log plot? How can you directly
read the order p and the constant C > 0 from a log-log plot of eh = C hp ?

1
Problem 2:
Consider a su”ciently smooth function f → C k [x0 ↗ ω, x0 + ω] for x0 → R, ω > 0 and k → N. Since
the first derivative f → is often unknown, one must resort to numerical methods to approximate the
value f → (x0 ) at the point x0 . One possibility is to use the one-sided di!erence quotient
f (x0 + h) ↗ f (x0 )
Dh1 f (x0 ) := ↘ f → (x0 ) for 0 < h < ω.
h
Write a Matlab script implementing the one-sided di!erence quotient for h → {2↑n : n = 0, . . . , N }
with N = 10. What rate of convergence do you observe for the error eh := |f → (x0 ) ↗ Dh1 f (x0 )| for
f1 (x) := exp(x), f2 (x) := sin(x), and f3 (x) := x3/2 at x0 := 0 for h ↓ 0? Represent the error
graphically in a log-log plot, where you plot eh over h.

Numerics of ODEs

The task is to compute a function y : [0, T ] ↓ Rd from a given initial value y0 → Rd and a
function f → C([0, T ] ≃ Rd ; Rd ) determining the derivative y → . The corresponding initial value
problem seeks y → C 1 ([0, T ]; Rd ) such that
" #
y → (t) = f t, y(t) for all t → [0, T ] and y(0) = y0 . (1)

A numerical solution (yε )L ε=1 approximates y(tε ) ↘ yε in certain nodes 0 = t0 < t1 < · · · <
tL = T , e.g., in equidistant nodes, tε+1 = tε + h for all ε = 0, . . . , L ↗ 1 and h := T /L. Taylor’s
theorem applied to the unknown function y and the ODE (1) shows, for r, s → [0, T ] with
|s ↗ r| ⇐ 1 , that
" #
y(s) = y(r) + (s ↗ r) y → (r) + O(|s ↗ r|2 ) ↘ y(r) + (s ↗ r) f r, y(r) . (2)

The choice s = t + h and r = t leads us to the explicit Euler method

" #
yε+1 = yε + h f tε , yε for all ε = 0, . . . , L ↗ 1. (3)

Alternatively, the choice s = t and r = t + h provides the implicit Euler method

" #
yε+1 = yε + h f tε+1 , yε+1 for all ε = 0, . . . , L ↗ 1. (4)

Moreover, one may consider the implicit midpoint rule defined by

$t + t %
ε ε+1 yε + yε+1
yε+1 = yε + h f , for all ε = 0, . . . , L ↗ 1. (5)
2 2
While yε+1 in (3) can always be computed, (4)–(5) provide only an implicit determination
of yε+1 (since it is involved on both sides of the equality). For usual problems, the Banach
fixed-point theorem proves that (4)–(5) have a unique solution yε+1 if h is su”ciently small.
A general one-step method is determined by the increment function # : [0, T ] ≃ Rd ≃ R+ ↓ Rd
in
yε+1 := yε + h #(tε , yε , h). (6)
All these methods can be considered for variable step sizes hε := tε+1 ↗ tε instead of a uniform
step size hε = h.

2
Problem 3:
For d = 1, the scalar initial value problem (1) is called linear if

f (t, y) = f1 (t) + f2 (t)y with scalar functions f1 , f2 : R ↓ R. (7)

Implement Matlab functions for the one-step methods (3)–(5) to solve a linear scalar initial value
problem. Adhere to the following signatures:

(a) y = explicitEuler(t, f1, f2, y0)

(b) y = implicitEuler(t, f1, f2, y0)

(c) y = implicitMidpoint(t, f1, f2, y0)

The arguments consist of the row vector of time steps t = (t0 , . . . , tL ) (not necessarily equidistant),
the function handles f1, f2 determining the right-hand side f , and the initial value y0. The output
y = (y0 , . . . , yL ) is the row vector containing the approximations to the solution y in the time steps
t, i.e., yε ↘ y(tε ).

Problem 4:
Consider the scalar initial value problem

y → (t) = ↗y(t) for all t → [0, 1] and y(0) = y0 := 1. (8)

Note that (8) is linear (7) with f1 (t) := 0 and f2 (t) := ↗1 for all t → [0, 1]. Compute approximations
to the solution with all three functions y = method(t, f1, f2, y0) from Problem 3. Use equidistant
time steps with step size h = 1/L. The exact solution to (8) reads y(t) = exp(↗t) and enables the
computation of the exact error eh := max |y(tε ) ↗ yε | in the maximum norm.
ε=1,...,L

(a) Plot the exact solution y and the three approximations (yε )ε in dependence of t into one plot.

(b) Plot the error eh in dependence of the step size h → {2↑n : n = 1, . . . , N } with N = 10 and
compare the convergence rates for the di!erent methods.

Problem 5:
Modify your three functions from Problem 3 for the solution of systems of linear initial value
problems (1) with f satisfying (7) for f1 : R ↓ Rd and f2 : R ↓ Rd↓d . The signatures of the functions
should remain the same as in Problem 3. For the row vector t = (t0 , . . . , tL ) → RL+1 as input, the
output y = (y0 , . . . , yL ) → Rd↓(L+1) should provide the columns y(:, ε + 1) = yε ↘ y(tε ) → Rd . Ahead
of the implementation, derive mathematically what has to be done. Deriving the mathematical
formula for the implementation, what is the di!erence in the application of the explicit and implicit
methods?

Problem 6:
Consider the linear system of ordinary di!erential equations
& ' & '
y1→ = ↗y2 y1 (0) = 1
with initial values (9)
y2→ = y1 y2 (0) = 0

3
Apply the three methods (3)–(5) using the functions from Problem 5 to solve the system (9). What
are the functions f1 and f2 in this example? The exact solution to (9) reads y(t) := (cos(t), sin(t))↔
and enables the computation of the exact error eh as in Problem 4 (with the absolute value | · |
replaced by the maximum norm ⇒ · ⇒↗ on Rd ).

(a) Plot the error eh := maxε=1,...,L ⇒y(tε ) ↗ yε ⇒↗ in dependence of the step size h → {2↑n : n =
1, . . . , N } with N = 10 and compare the convergence rates for the di!erent methods.

(b) Plot the exact and the discrete solution in a phase field plot, i.e., the component yε,1 (on the
y-axis) in dependence of yε,2 (on the x-axis). What does the exact solution look like in the
phase field plot? What is the di!erence between the three one-step methods? Use a larger
step size, e.g., h = 0.1, to see the di!erences more clearly.

Runge–Kutta methods

The most important class of one-step methods are the so-called Runge–Kutta methods. Let
A → Rm↓m and b, c → Rm with 0 ↔ c1 ↔ c2 ↔ · · · ↔ cm ↔ 1. An m-stage Runge–Kutta method
is a one-step method in the sense of (6) with the increment function
m
(
#(t, y, h) := bj k j , (10)
j=1

where the so-called stages satisfy the implicit conditions

$ m
( %
kj = f t + cj h, y + h Ajε kε for all j = 1, . . . , m. (11)
ε=1

The method is called explicit if A is a strictly lower triangular matrix, i.e., Ajε = 0 for all
j ↔ ε. Otherwise, it is called implicit. Usually, a Runge–Kutta method is written in terms of
c A
its Butcher tableau , where zero entries of A are usually neglected. As for the implicit
b
methods (4)–(5), one can show that for usual problems the implicit formulation (11) admits
a unique solution k = (k1 , . . . , km ) → Rm by means of the Banach fixed-point theorem.

Problem 7:
Implement a Matlab function

y = explicitRungeKutta(A, b, c, t, f1, f2, y0)

for an explicit m-step Runge–Kutta method given by the strictly lower triangular matrix A → Rm↓m
and b, c → Rm , applied to a scalar and linear initial value problem (see Problem 4).

(a) Which of the three methods (3)–(5) are Runge–Kutta methods and which can be implemented
by your explicitRungeKutta function? What are the respective Butcher tableaus?

(b) Compare the new function explicitRungeKutta.m to the direct implementation from Prob-
lem 3.

4
Problem 8:
Extend your implementation of explicit Runge–Kutta methods from Problem 7 to general linear
initial value problems, i.e., general d ⇑ 1; cf. Problem 5.

Problem 9:
Use your implementation from Problem 8 to solve the initial value problem (9) from Problem 6.
Test the classical RK4 method given by the Butcher tableau

0
1/2 1/2
c A
b
= 1/2 0 1/2
1 0 0 1
1/6 1/3 1/3 1/6

(a) Compute the error and compare the convergence rate with the results from Problem 6 (a).
What convergence rate do you observe?

(b) Plot the solution into the phase field plot from Problem 6 (b).

Problem 10:
Implement a Matlab function y = implicitRungeKutta(A, b, c, t, f1, f2, y0) for arbitrary (implicit)
Runge–Kutta methods for a scalar and linear initial value problem (see Problem 4).
Hint: Proceed as for the implicit Euler method (4) in Problem 3 and write the stage-conditions
(11) as an m ≃ m linear system with solution k = (k1 , . . . , km )↔ → Rm .

Problem 11:
Consider the 2-stage and 3-stage Gauss–Runge–Kutta methods
⇓ ⇓
(3 ↗ ⇓3)/6 1/4
⇓ (3 ↗ 2 3)/12
c A
b
= (3 + 3)/6 (3 + 2 3)/12 1/4
1/2 1/2
⇓ ⇓ ⇓
1/2 ↗ 15/10 5/36
⇓ 2/9 ↗ 15/15 5/36 ↗ ⇓15/30
c A
b
= 1/2
⇓ 5/36 + ⇓15/24 2/9
1/2 + 15/10 5/36 + 15/30 2/9 + 15/15
⇓ 5/36 ↗ 15/24
5/36
5/18 4/9 5/18

What convergence rates do you observe for the initial value problem from Problem 4.

5
Take-home message

Explicit Runge–Kutta methods (like the explicit Euler method) to solve ODEs

y → = f (t, y(t))

can easily be implemented for arbitrary (nonlinear) f and any dimension d ⇑ 1. Implicit
Runge–Kutta methods (like the implicit Euler method or the midpoint scheme) lead to
nonlinear systems of equations that must be solved to compute the stages. If f is linear in
the sense of (7), then

1-step Runge–Kutta methods for arbitrary d ⇑ 1,

m-step Runge–Kutta methods for d = 1

require the solution of one linear system per time step. Actually, the same holds indeed for
arbitrary m ⇑ 1 and d ⇑ 1.