WS07: Improved Euler’s methods - partial solutions

WS07: Improved Euler’s methods - partial solutions#

These exercises are indented to give you practice at using the material on numerical approximation and are intended to reinforce the material that was covered in lectures.

Please attempt the worksheet before your tutorial. Support is available in your tutorial or in the Class Team.

These are partial solutions. Please create Issues and Pull requests with your solutions.

Part a (pen and paper warm up)#

1. Midpoint method#

For a differential equation

\[ u^\prime(t) = f(t,u), \]

given a time step \(\Delta t = t_1 - t_0\) and \(u(t_0)=u_0\), \(u(t_1)=u_1\) may be computed by

  • the Euler’s formula

    \[ \frac{u_1-u_0}{\Delta t} = f(t_0,u_0) \]

  • or the midpoint formula

    \[ \frac{u_1-u_0}{\Delta t} = f(t_{1/2},u_{1/2}), \]

    where

    \[ t_{1/2} = \frac{t_0+t_1}{2} \]

    and \(u_{1/2}\)

    \[ u_{1/2} = u_0 + \frac{\Delta t}{2} f(t_0,u_0), \]

    which is computed by the Euler’s formula on a half-time step of \(\frac{\Delta t}{2}\).

Both formulas are approximations of the original differential equation. However, the midpoint formula has one more step to compute a midpoint \((t_{1/2}, u_{1/2})\), which is more accurate. Apply the midpoint formula to the same problem considered in WS06:

\[ d^\prime(t) = t^2 -3t + d \]

with \(\Delta t = 0.01\), and \(d(0) =1\). Use the midpoint method to compute \(d(0.01)\).

Part b (code implementations and testing)#

2. Implement the midpoint method#

import numpy as np

def midpoint_method(t0, d0, dt, n):
    d = np.zeros(n + 1)
    d[0] = d0

    t = np.zeros(n + 1)
    for i in range(n + 1):
        t[i] = t0 + i * dt

    for i in range(1, n + 1):
        d_half = d[i - 1] + 0.5 * dt * f(t[i - 1], d[i - 1])
        d[i] = d[i - 1] + dt * f((t[i - 1] + t[i]) / 2.0, d_half)

    return t, d


def Euler_method(t0, d0, dt, n):
    d = np.zeros(n + 1)
    d[0] = d0

    t = np.zeros(n + 1)
    for i in range(n + 1):
        t[i] = t0 + i * dt

    for i in range(1, n + 1):
        d[i] = d[i - 1] + dt * f(t[i - 1], d[i - 1])

    return t, d

3. Test your implementation#

Use the example in Q4 to test your implementation and compare with the exact solution: \(d(t) = -t^2 +t +1\)

def f(t, d):
    return t * t - 3.0 * t + d


def d_exact(t):
    return -t * t + t + 1.0


t0 = float(0.0)
d0 = float(1.0)
dt = float(0.01)
n = int(100)

t_Euler, d_Euler = Euler_method(t0, d0, dt, n)
t_midpoint, d_midpoint = midpoint_method(t0, d0, dt, n)

# test against your solution from part a
d_parta_Euler = np.array([1.0, 1.01, 1.019801])
np.testing.assert_allclose(d_Euler[0:3], d_parta_Euler)
d_parta_midpoint = np.array([1.0, 1.00990025, 1.0196005025125])
np.testing.assert_allclose(d_midpoint[0:3], d_parta_midpoint)

de = np.zeros(n + 1)
for i in range(n + 1):
    de[i] = d_exact(t_Euler[i])


import matplotlib.pyplot as plt

fig = plt.figure()
plt.plot(t_Euler, d_Euler, "-b", label="Euler")
plt.plot(t_Euler, de, "-r", label="Exact")
plt.plot(t_midpoint, d_midpoint, "--g", label="midpoint")
plt.legend(loc="upper left")
plt.xlabel("t")
plt.ylabel("d")
plt.grid()
plt.show()

4. Test a larger time step#

def f(t, d):
    return t * t - 3.0 * t + d


def d_exact(t):
    return -t * t + t + 1.0


t0 = float(0.0)
d0 = float(1.0)
dt = float(0.1)
n = int(10)

t, d_Euler = Euler_method(t0, d0, dt, n)
t, d_midpoint = midpoint_method(t0, d0, dt, n)

de = np.zeros(n + 1)
for i in range(n + 1):
    de[i] = d_exact(t[i])


import matplotlib.pyplot as plt

fig = plt.figure()
plt.plot(t, d_Euler, "-b", label="Euler")
plt.plot(t, de, "-r", label="Exact")
plt.plot(t, d_midpoint, "--g", label="midpoint")
plt.legend(loc="upper left")
plt.xlabel("t")
plt.ylabel("d")
plt.grid()
plt.show()

Part c: Extension#

Apply your Python code of Euler’s method to equation

\[ d^\prime(t) = 2 \sqrt d \quad \text{starting from} \quad d(0)=0. \]

The exact solution is \(d(t) = t^2\).

  • Check whether your code will fail.

  • Explain why it fails and how you can avoid this?

def f(t, d):
    return 2.0 * d**0.5


def d_exact(t):
    return t**2.0


# starting at the proper initial condition
t0 = float(0.0)
d0 = float(0.0)
dt = float(0.1)
n = int(10)

t, d_Euler = Euler_method(t0, d0, dt, n)
t, d_midpoint = midpoint_method(t0, d0, dt, n)

de = np.zeros(n + 1)
for i in range(n + 1):
    de[i] = d_exact(t[i])

import matplotlib.pyplot as plt

fig = plt.figure()
plt.plot(t, d_Euler, "-b", label="Euler")
plt.plot(t, de, "-r", label="Exact")
plt.plot(t, d_midpoint, "--g", label="midpoint")
plt.legend(loc="upper left")
plt.xlabel("t")
plt.ylabel("d")
plt.grid()
plt.show()

# starting at a tweaked initial condition
t0 = float(0.1)
d0 = float(0.01)
dt = float(0.1)
n = int(10)

t, d_Euler = Euler_method(t0, d0, dt, n)
t, d_midpoint = midpoint_method(t0, d0, dt, n)

de = np.zeros(n + 1)
for i in range(n + 1):
    de[i] = d_exact(t[i])

import matplotlib.pyplot as plt

fig = plt.figure()
plt.plot(t, d_Euler, "-b", label="Euler")
plt.plot(t, de, "-r", label="Exact")
plt.plot(t, d_midpoint, "--g", label="midpoint")
plt.legend(loc="upper left")
plt.xlabel("t")
plt.ylabel("d")
plt.grid()
plt.show()
Contents