What is the second-order approximation error of the Heun method, which averages the slope at the beginning and the end of the interval?
$y(x+h) = y(x) + \frac{h}{2}\{f(x,y(x))+f[x+h,y(x)+hf(x,y(x))]\}$
$y(x+h) = y(x) + \frac{h}{2}f(x,y(x)) + \frac{h}{2}f[x+h, y(x) + hf(x,y(x))]$
$\frac{h}{2}f[x+h, y(x) + hf(x,y(x))] = f(x,y(x)) + h \frac{d}{dh} f[x+h, y(x)+hf(x,y(x))]|_{h=0} + \mathcal{O}(h^2)$
$= f(x,y(x)) + h[1\frac{\partial f}{\partial x} + f(x, y(x)) \frac{\partial f}{\partial y}] + \mathcal{O}(h^2)$
$\implies y(x+h) = y(x) + \frac{h}{2}f(x,y(x)) + \frac{h}{2}\{f(x,y(x)) + h[\frac{\partial f}{\partial x} + f(x,y(x))\frac{\partial f}{\partial y}] + \mathcal{O}(h^2)\}$
$y(x+h) = y(x) + hf(x,y(x)) + \frac{h^2}{2}[\frac{\partial f}{\partial x} + f(x,y(x))\frac{\partial f}{\partial y}] + \mathcal{O}(h^3)$
We can see that this is identical to the actual Taylor expansion of the problem, and therefore, there is no difference in the first terms. The error lies in the other higher terms, giving us an error $\mathcal{O}(h^3)$.
For a simple harmonic oscillator $\ddot{y}+y=0$, with initial conditions $y(0)=1,\ \dot{y}(0)=0$, find $y(t)$ from $t=0$ to $100\pi$. Use an Euler method and a fixed-step fourth-order Runge-Kutta method. For each method check how the average local error, and the error in the final value and slope, depend on the step size.
# Adapted from Neil's euler.c code
import numpy as np
steps = [1, 0.1, 0.01, 0.001, 0.0001, 0.00001]
for h in steps:
t_max = 100*np.pi
t = 0
N = int(t_max / h)
y0 = 1
y1 = 0
error = 0
while (t < t_max):
new_y0 = y0 + h*y1
new_y1 = y1 - h*y0
y0 = new_y0
y1 = new_y1
t = t + h
error = error + np.absolute(np.cos(t) - y0)
error = error/N
final_err = np.absolute(np.cos(t) - y0)
print("Step Size: %f Error: %e Final Error: %e Slope: %e\n" % (h, error, final_err, y1))
# Adapted from Neil's rk.c code
import numpy as np
steps = [1, 0.1, 0.01, 0.001, 0.0001, 0.00001]
for h in steps:
t_max = 100*np.pi
t = 0
N = int(t_max / h)
y0 = 1
y1 = 0
error = 0
while (t < t_max):
A0 = h*y1
A1 = -h*y0
B0 = h*(y1 + A1/2.0)
B1 = -h*(y0 + A0/2.0)
C0 = h*(y1 + B1/2.0)
C1 = -h*(y0 + B0/2.0)
D0 = h*(y1 + C1)
D1 = -h*(y0 + C0)
y0 = y0 + (A0 + 2.0*B0 + 2.0*C0 + D0)/6.0
y1 = y1 + (A1 + 2.0*B1 + 2.0*C1 + D1)/6.0
t = t + h
error = error + np.absolute(np.cos(t) - y0)
error = error / N
final_err = np.absolute(np.cos(t) - y0)
print("Step Size: %f Error: %e Final Error: %e Slope: %e\n" % (h, error, final_err, y1))
Write a fourth-order Runge-Kutta adaptive stepper for the preceding problem, and check how the average step size that it finds depends on the desired local error.
# Adapted from Neil's rkstep.c code
import numpy as np
def step(y0, y1, h):
A0 = h*y1
A1 = -h*y0
B0 = h*(y1 + A1/2.0)
B1 = -h*(y0 + A0/2.0)
C0 = h*(y1 + B1/2.0)
C1 = -h*(y0 + B0/2.0)
D0 = h*(y1 + C1)
D1 = -h*(y0 + C0)
y0_final = y0 + (A0 + 2.0*B0 + 2.0*C0 + D0)/6.0
y1_final = y1 + (A1 + 2.0*B1 + 2.0*C1 + D1)/6.0
return y0_final, y1_final
target_error = [1, 0.1, 0.01, 0.001, 0.0001, 0.00001, 0.000001, 0.0000001, 0.00000001, 0.000000001, 0.0000000001]
for target in target_error:
t_max = 100*np.pi
t = 0
N = 0
y0 = 1
y1 = 0
error = 0
err_avg = 0
h = 0.1
h_avg = 0
Y1_0 = 0
Y1_1 = 0
Y2_0 = 0
Y2_1 = 0
Y3_0 = 0
Y3_1 = 0
while (t < t_max):
N = N + 1
Y1_0, Y1_1 = step(y0, y1, (h/2.0))
Y2_0, Y2_1 = step(Y1_0, Y1_1, (h/2.0))
Y3_0, Y3_1 = step(y0, y1, h)
local_err = Y3_0 - Y2_0
if (np.absolute(local_err) > target):
h = h / 1.3
continue
y0 = Y2_0
y1 = Y2_1
t = t + h
err_avg = err_avg + np.absolute(local_err)
error = error + np.absolute(np.cos(t) - y0)
h_avg = h_avg + h
if (np.absolute(local_err) < (target/10.0)):
h = 1.2 * h
error = error / N
h_avg = h_avg / N
err_avg = err_avg / N
final_err = np.absolute(np.cos(t) - y0)
print("Target Error: %e Avg. Step Size: %f Avg. Error: %e \n Final Error: %e Avg. Local Error: %e\n" % (target, h_avg, error, final_err, err_avg))
Numerically solve the differential equation found in Problem 5.3:
$l\ddot{\theta}+(g+\ddot{z})sin\theta=0$
Take the motion of the platform to be periodic, and interactively explore the dynamics of the pendulum as a function of the amplitude and frequency of the excitation.
See Neil's pendulum.html example of this in motion.