7.1 - test
What is the 2nd-order approximation error of the Heun method, which averages the slope at the beginning and the end of the interval?
Solution
Starting with the basic form of the solve step:
$$ y(x+h) = y(x) + \frac{h}{2}\bigg(f\big(x, y(x)\big) + f\big(x + h, y(x) + hf(x, y(x))\big)\bigg) $$
and taylor-expanding out the calculation of the slope at the endpoint: $$ \begin{array} \\ f\big(x + h, y(x) + hf(x, y(x))\big) &= f\big(x, y(x)\big) + h\frac{d}{dh} f\big[x + h, y(x) + hf(x, y(x))\big] + \mathcal{O}\left(h^2\right) \\ &= f\big(x, y(x)\big) + h\left[\frac{\partial f}{\partial x} + \frac{\partial f}{\partial y}f(x, y(x))\right] + \mathcal{O}\left(h^2\right) \\ \end{array} $$
Substituting back into the original equation gives us: $$ \begin{array} \\ y(x+h) &= y(x) + \frac{h}{2}\big(f(x, y(x)\big) + \frac{h}{2}\big(f(x, y(x)\big) + \frac{h^2}{2}\left[\frac{\partial f}{\partial x} + \frac{\partial f}{\partial y}f(x, y(x))\right] + \mathcal{O}\left(h^3\right) \\ &= y(x) + h\big(f(x, y(x)\big) + \frac{h^2}{2}\left[\frac{\partial f}{\partial x} + \frac{\partial f}{\partial y}f(x, y(x))\right] + \mathcal{O}\left(h^3\right) \end{array} $$
Which we can see matches the first two terms of the tylor expansion, so the 2nd-order error is zero.
7.2
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.
Solution
Start by converting the 2nd-order ODE to a system of 1st-order ODEs. If we set $\dot y = y_2$ and $y = y_1$, then our system is:
$$ \begin{array} \\ \dot y_1 = y_2 \\ \dot y_2 = -y_1 \\ \end{array} $$
Or in matrix form:
$$ \left[ \begin{array} \\ \dot y_1 \\ \dot y_2 \end{array} \right] = \left[ \begin{array}{cc} \\ 0 & 1 \\ -1 & 0 \end{array} \right] \left[ \begin{array} \\ y_1 \\ y_2 \end{array} \right] $$
Implementation
make_system() is a convenience function that takes a single higher-order ODE and builds the appropriate matrix representing a system of 1st-order ODEs
function make_system(coefs)
# Takes in a vector of coeficients V representing a high order ODE, e.g. x''' = v1 x + v2 x' + v3 x'' and
# returns a matrix representing the corresponding system of 1st-order ODEs
ord = length(coefs)
top = hcat(zeros(ord-1), eye(ord-1, ord-1))
return vcat(top, coefs')
end
Euler's Method
euler() implements Euler's algorithm for solving a system of 1st-order ODEs. It takes a system matrix A, an initial value vector init, and a set of timesteps steps. The code is pretty direct, with a couple of tweaks to optimize speed and memory usage, namely pre-allocating temporary arrays outside the loop, and copying into them element-by-element instead of using slices which currently in julia cause a copy to be allocated for the slice. There's active work going on in the Julia community to use views for slices, which would allow the below code to be much cleaner.
function colcopy!(v, M, col)
# copies the col-th column of M into v without creating an intermediate slice
for i in 1:size(M, 1)
v[i] = M[i, col]
end
nothing
end
function euler(A, init, steps)
N = length(steps)
M = length(init)
y = zeros(M, N)
# pre-allocate outside the loop
slope = zeros(M, 1)
current = zeros(M, 1)
# set the initial conditions
y[:, 1] = init
for i in 1:(N-1)
colcopy!(current, y, i)
A_mul_B!(slope, A, current)
h = steps[i+1] - steps[i]
for j in 1:M
y[j, i+1] = y[j, i] + h * slope[j]
end
end
return y
end
Here we see an example solving with the Euler method. With a step size of 0.1 we clearly see the solution blowing up, as the system is right on the edge of stability.
using PyPlot
step = 0.1
steps = 0:step:100pi
# set up the system of equations
A = make_system([-1, 0])
init = [1, 0]
println("Calculating $(length(steps)) steps")
@time y = euler(A, init, steps)
plot(steps, y[1,:]')
#plot(steps, y[2,:]')
With a step size of 0.0001 we finally see something more closely approximating the right answer, obviously as the expense of execution time.
step = 0.0001
steps = 0:step:100pi
# set up the system of equations
A = make_system([-1, 0])
init = [1, 0]
println("Calculating $(length(steps)) steps")
@time y = euler(A, init, steps)
plot(steps, y[1,:]')
Fixed-Step Runge-Kutta
Now we'll implement a 4th-order Runge-Kutta solver using a fixed step size, as above with Euler. Because this can handle much larger step sizes it was less important to do the optimizations I did in the Euler implementation, so the code is a lot easier to read.
function rk4_fixed(A, init, steps)
N = length(steps)
M = length(init)
y = zeros(M, N)
# pre-allocate outside the loop
k1 = zeros(M, 1)
k2 = zeros(M, 1)
k3 = zeros(M, 1)
k4 = zeros(M, 1)
# set the initial conditions
y[:, 1] = init
for i in 1:(N-1)
h = steps[i+1] - steps[i]
k1[:] = h * A * y[:, i]
k2[:] = h * A * (y[:, i] .+ k1 ./ 2)
k3[:] = h * A * (y[:, i] .+ k2 ./ 2)
k4[:] = h * A * (y[:, i] .+ k3)
y[:, i+1] = y[:, i] + k1/6 + k2/3 + k3/3 + k4/6
end
return y
end
function rk4_fixed_fast(A, init, steps)
N = length(steps)
M = length(init)
y = zeros(M, N)
# pre-allocate outside the loop
k1 = zeros(M, 1)
k2 = zeros(M, 1)
k3 = zeros(M, 1)
k4 = zeros(M, 1)
current = zeros(M, 1)
temp = zeros(M, 1)
# set the initial conditions
y[:, 1] = init
for i in 1:(N-1)
h = steps[i+1] - steps[i]
colcopy!(current, y, i)
# calc k1
for j in 1:M
temp[j] = h * current[j]
end
A_mul_B!(k1, A, temp)
# calc k2
for j in 1:M
temp[j] = h * (current[j] + k1[j] / 2)
end
A_mul_B!(k2, A, temp)
# calc k3
for j in 1:M
temp[j] = h * (current[j] + k2[j] / 2)
end
A_mul_B!(k3, A, temp)
# calc k4
for j in 1:M
temp[j] = h * (current[j] + k3[j])
end
A_mul_B!(k4, A, temp)
for j in 1:M
y[j, i+1] = current[j] + k1[j]/6 + k2[j]/3 + k3[j]/3 + k4[j]/6
end
end
return y
end
step = 0.001
steps = 0:step:100pi
# set up the system of equations
A = make_system([-1, 0])
init = [1, 0]
println("Calculating $(length(steps)) steps")
@time y = rk4_fixed(A, init, steps)
plot(steps, y[1,:]')
Comparing Euler and Runge-Kutta
Here we'll compare the performance of the two algorithms for a variety of step sizes
da = [1, 2, 3, NaN, 5]
filter(x -> !is(x, NaN), da)
using DataFrames
function dropnan(v)
return filter(x -> !is(x, NaN), v)
end
# set up the dataframe we'll be using to store the statistics
step_sizes = [1, 0.1, 0.01, 0.001, 0.0001]
euler_data = DataFrame()
rk4_data = DataFrame()
euler_data[:step] = step_sizes
rk4_data[:step] = step_sizes
for column in [:avg_value_err, :avg_slope_err, :final_value_err, :final_slope_err]
euler_data[column] = zeros(length(step_sizes))
rk4_data[column] = zeros(length(step_sizes))
end
# set up the actual system we'll be simulating
A = make_system([-1, 0])
init = [1, 0]
for i in 1:length(step_sizes)
steps = 0:step_sizes[i]:100pi
println("Calculating $(length(steps)) steps using Euler")
@time euler_y = euler(A, init, steps)
println("Calculating $(length(steps)) steps using RK4")
@time rk4_y = rk4_fixed_fast(A, init, steps)
truth = vcat(cos(steps)', -sin(steps)')
euler_err = (euler_y - truth) ./ truth * 100
rk4_err = (rk4_y - truth) ./ truth * 100
euler_data[:final_value_err][i] = euler_err[1, end]
euler_data[:final_slope_err][i] = euler_err[2, end]
euler_data[:avg_value_err][i] = mean(dropnan(abs(euler_err[1, :])))
euler_data[:avg_slope_err][i] = mean(dropnan(abs(euler_err[2, :])))
rk4_data[:final_value_err][i] = rk4_err[1, end]
rk4_data[:final_slope_err][i] = rk4_err[2, end]
rk4_data[:avg_value_err][i] = mean(dropnan(abs(rk4_err[1, :])))
rk4_data[:avg_slope_err][i] = mean(dropnan(abs(rk4_err[2, :])))
end
# wait for the printf displays to catch up
sleep(0.1)
display("text/html", "<h3>Euler Error (%)</h3>")
display(euler_data)
display("text/html", "<h3>RK4 Error (%)</h3>")
display(rk4_data)
7.3
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
function rk4_adaptive(A, init, err_min, err_max, t_min, t_max, init_step=0.0001, downstep=0.5, upstep=1.4)
@assert all(err_min .< err_max)
@assert t_min < t_max
# we'll start with an initial N to allocate arrays
N = 10000
steps = Array(Float64, N)
M = length(init)
y = Array(Float64, M, N)
# pre-allocate outside the loop
k1 = Array(Float64, M, 1)
k2 = Array(Float64, M, 1)
k3 = Array(Float64, M, 1)
k4 = Array(Float64, M, 1)
y_whole = Array(Float32, M, 1)
y_half = Array(Float32, M, 1)
# set the initial conditions
i = 1
y[:, 1] = init
t = t_min
steps[1] = t
step_size = init_step
while t <= t_max
# check to see whether we need to resize our data arrays
if i >= N
info("resizing arrays to $(2*N)")
resize!(steps, N * 2)
y = hcat(y, Array(Float64, M, N))
N *= 2
end
err = infs(Float64, M)
while any(err .> err_max)
# calculate in two half-steps
step_rk4!(y_half, y[:, i], A, step_size / 2, k1, k2, k3, k4)
step_rk4!(y_half, y_half, A, step_size / 2, k1, k2, k3, k4)
# calculate in one whole-step
step_rk4!(y_whole, y[:, i], A, step_size, k1, k2, k3, k4)
err = y_whole - y_half
if any(err .> err_max)
step_size *= downstep
end
end
y[:, i+1] = y_half + (y_half - y_whole) / 15
i += 1
t += step_size
if all(err .< err_min)
step_size *= upstep
end
steps[i] = t
end
return y, steps, i
end
function step_rk4!(y_next, y_current, A, h, k1, k2, k3, k4)
k1[:] = h .* A * y_current
k2[:] = h .* A * (y_current .+ k1 ./ 2)
k3[:] = h .* A * (y_current .+ k2 ./ 2)
k4[:] = h .* A * (y_current .+ k3)
y_next[:] = y_current + k1/6 + k2/3 + k3/3 + k4/6
end
# set up the system of equations
A = make_system([-1, 0])
init = [1, 0]
@time y, steps, N = rk4_adaptive(A, init, [1e-3, 1e-3], [1e-2, 1e-2], 0, 20pi)
plot(steps[1:N], y[1,1:N]')
scatter(steps[1:N], y[1,1:N]')
plot(0:0.0001:20pi, cos(0:0.0001:20pi))
7.4
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.
Solution
Because $z$ is periodic, we'll consider it to be of the form $z(t) = A \cos(\omega t)$, so $\ddot z(t) = -A \omega^2 \cos(\omega t)$.
Now, as before we turn this 2nd-order system into a set of 1st order equations by setting:
$$ \begin{array} \\ \dot\theta_2 = \frac{-g - \ddot z}{l} \sin(\theta_1) \\ \dot\theta_1 = \theta_2 \\ \end{array} $$
Substituting in our expression for $\ddot z$ we get the function we evaluate at each step of the solver:
$$ \begin{array} \\ \dot\theta_2 = \frac{A \omega^2 cos(\omega t) - g}{l} \sin(\theta_1) \\ \dot\theta_1 = \theta_2 \\ \end{array} $$
function rk4_adaptive_f(f, init, err_min, err_max, t_min, t_max, init_step=0.0001, downstep=0.5, upstep=1.4)
@assert all(err_min .< err_max)
@assert t_min < t_max
# we'll start with an initial N to allocate arrays
N = 10000
steps = Array(Float64, N)
M = length(init)
y = Array(Float64, M, N)
# pre-allocate outside the loop
k1 = Array(Float64, M, 1)
k2 = Array(Float64, M, 1)
k3 = Array(Float64, M, 1)
k4 = Array(Float64, M, 1)
y_whole = Array(Float32, M, 1)
y_half = Array(Float32, M, 1)
# set the initial conditions
i = 1
y[:, 1] = init
t = t_min
steps[1] = t
step_size = init_step
while t <= t_max
# check to see whether we need to resize our data arrays
if i >= N
info("resizing arrays to $(2*N)")
resize!(steps, N * 2)
y = hcat(y, Array(Float64, M, N))
N *= 2
end
err = infs(Float64, M)
while any(err .> err_max)
# calculate in two half-steps
step_rk4_f!(f, y_half, y[:, i], t, step_size / 2, k1, k2, k3, k4)
step_rk4_f!(f, y_half, y_half, t + step_size / 2, step_size / 2, k1, k2, k3, k4)
# calculate in one whole-step
step_rk4_f!(f, y_whole, y[:, i], t, step_size, k1, k2, k3, k4)
err = y_whole - y_half
if any(err .> err_max)
step_size *= downstep
end
end
y[:, i+1] = y_half + (y_half - y_whole) / 15
i += 1
t += step_size
if all(err .< err_min)
step_size *= upstep
end
steps[i] = t
end
return y, steps, i
end
function step_rk4_f!(f, y_next, y_current, t, h, k1, k2, k3, k4)
k1[:] = h .* f(y_current, t)
k2[:] = h .* f(y_current .+ k1 ./ 2, t + h/2)
k3[:] = h .* f(y_current .+ k2 ./ 2, t + h/2)
k4[:] = h .* f(y_current .+ k3, t + h)
y_next[:] = y_current + k1/6 + k2/3 + k3/3 + k4/6
end
function platform_pendulum!(A, omega, g, l, state, t)
return [state[2],
(A * omega^2 * cos(omega * t) - g) / l * sin(state[1])]
end
A = 1
omega = 2
g = 9.8
l = 8
step_func!(state, t) = platform_pendulum!(A, omega, g, l, state, t)
init = [0.1, 0]
@time y, steps, N = rk4_adaptive_f(step_func!, init, [1e-9, 1e-9], [1e-5, 1e-5], 0, 100pi)
plot(steps[1:N], y[1,1:N]')
scatter(steps[1:N], y[1,1:N]')
