# 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

In [5]:
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

Out[5]:
make_system (generic function with 1 method)


#### 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.

In [6]:
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

Out[6]:
euler (generic function with 1 method)


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.

In [7]:
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,:]')  Calculating 3142 steps elapsed time: 0.801462585 seconds (28475876 bytes allocated)  Out[7]: 1-element Array{Any,1}: PyObject <matplotlib.lines.Line2D object at 0x112d3b250>  With a step size of 0.0001 we finally see something more closely approximating the right answer, obviously as the expense of execution time. In [4]: 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,:]')

Calculating 3141593 steps
elapsed time: 1.590679661 seconds (50265840 bytes allocated)


Out[4]:
1-element Array{Any,1}:
PyObject <matplotlib.lines.Line2D object at 0x10ed96ed0>


#### 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.

In [8]:
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

Out[8]:
rk4_fixed_fast (generic function with 1 method)

In [14]:
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,:]')  Calculating 314160 steps elapsed time: 11.201457385 seconds (2018076248 bytes allocated)  Out[14]: 1-element Array{Any,1}: PyObject <matplotlib.lines.Line2D object at 0x112d56b90>  ### Comparing Euler and Runge-Kutta Here we'll compare the performance of the two algorithms for a variety of step sizes In [142]: da = [1, 2, 3, NaN, 5] filter(x -> !is(x, NaN), da)  Out[142]: 4-element Array{Float64,1}: 1.0 2.0 3.0 5.0  In [178]: 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)  Calculating 315 steps using Euler elapsed time: 0.000191948 seconds (5328 bytes allocated) Calculating 315 steps using RK4 elapsed time: 0.000544636 seconds (5712 bytes allocated) Calculating 3142 steps using Euler elapsed time: 0.001399866 seconds (50560 bytes allocated) Calculating 3142 steps using RK4 elapsed time: 0.005111929 seconds (50944 bytes allocated) Calculating 31416 steps using Euler elapsed time: 0.012768524 seconds (502944 bytes allocated) Calculating 31416 steps using RK4 elapsed time: 0.049974957 seconds (503328 bytes allocated) Calculating 314160 steps using Euler elapsed time: 0.125643699 seconds (5026848 bytes allocated) Calculating 314160 steps using RK4 elapsed time: 0.495246999 seconds (5027232 bytes allocated) Calculating 3141593 steps using Euler elapsed time: 1.299615605 seconds (50265776 bytes allocated) Calculating 3141593 steps using RK4 elapsed time: 4.887341487 seconds (50266160 bytes allocated)  ### Euler Error (%) stepavg_value_erravg_slope_errfinal_value_errfinal_slope_err 11.02.0148161167415956e471.2680345987231057e48-100.0-1.1519286028685861e50 20.11.5989170591796696e81.374381033507333e82.780520228473212e89.207798228018274e9 30.01149.17912042570063147.3189389418824380.9146456266811924.5132999228988 40.0018.2958918317640728.30952340405123417.00885014207598263.18449275239829 50.00010.78970402641744590.7896983070025681.58319783621921763.210788331515642 ### RK4 Error (%) stepavg_value_erravg_slope_errfinal_value_errfinal_slope_err 11.0106.97899365615065106.39149551606869-104.95681885479662-12.828661099133823 20.10.37639293840673610.05919084925115199-0.00372941356465829620.4373744930627087 30.011.075961364506886e-58.539312291995246e-6-4.607172555919578e-80.00028250878940007287 40.0019.977816279269523e-101.3011024797826422e-95.10702609308233e-129.694741279885758e-7 50.00015.873245571624272e-115.658716707043306e-11-3.474998074498982e-12-1.9662477687929266e-7 # 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 In [10]: 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

Out[10]:
step_rk4! (generic function with 1 method)

In [11]:
# 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))

elapsed time: 0.885063959 seconds (32087704 bytes allocated)


Out[11]:
1-element Array{Any,1}:
PyObject <matplotlib.lines.Line2D object at 0x111fea910>


# 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}$$

In [2]:
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

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

Out[2]:
step_rk4_f! (generic function with 1 method)

In [13]:
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]')

elapsed time: 0.379127041 seconds (68744384 bytes allocated)


Out[13]:
PyObject <matplotlib.collections.PathCollection object at 0x1136804d0>