# Week 3 - Finite Differences: ODE¶

## Problem 7.1:¶

#### 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)))]$$

Now we can use the Taylor expansion for $f(x+h, y(x) + hf(x,y(x)))$:

$$\begin{split} 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))]_0 + O(h^2) \\ &= f(x,y(x)) + h[\frac{\partial f}{\partial x} + f(x,y(x) \frac{\partial f}{\partial y}] + O(h^2) \end{split}$$

Plugging back into the original \textit{Heun} equation:

$$\begin{split} y(x+h) &= y(x) + \frac{h}{2}[f(x,y(x)) + h[\frac{\partial f}{\partial x} + f(x,y(x) \frac{\partial f}{\partial y}] + O(h^2)] \ &= y(x) + hf(x,y(x)) + \frac{h^2}{2} [\frac{\partial f}{\partial x} + f(x,y(x) \frac{\partial f}{\partial y}] + O(h^3) \end{split}$$

And the error is $O(h^3)$

## Problem 7.2:¶

#### For a simple harmonic oscillator y¨+y = 0, with initial conditions y(0) = 1, y˙(0) = 0, find y(t) from t = 0 to 100π. 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.¶

In [31]:
function euler(start, end, init1, init2, step){
var y = [init1]
var z = [init2]

for (i = start; i <= end-step; i += step){
y.push(y[y.length-1] + step*z[z.length-1])
z.push(z[z.length-1] - step*y[y.length-1])
}

return y
}

function runge_kutta(start, end, init1, init2, step){
var y = [init1]
var z = [init2]

for (i = start; i <= end-step; i += step){
var k_1 = step*z[z.length-1]
var l_1 = -step*y[y.length-1]
var k_2 = step*(z[z.length-1] + 1/2*l_1)
var l_2 = -step*(y[y.length-1] + 1/2*k_1)
var k_3 = step*(z[z.length-1] + 1/2*l_2)
var l_3 = -step*(y[y.length-1] + 1/2*k_2)
var k_4 = step*(z[z.length-1] + l_3)
var l_4 = -step*(y[y.length-1] + k_3)

y.push(y[y.length-1] + 1/6*(k_1 + 2*k_2 + 2*k_3 + k_4))
z.push(z[z.length-1] + 1/6*(l_1 + 2*l_2 + 2*l_3 + l_4))
}

return y
}

function real(start, end, step){
y = []

for (i = start; i <= end; i += step){
y.push(Math.cos(i))
}

return y
}

var start = 0
var end = 100*Math.PI

var eu_end_err = []
var rk_end_err = []
var eu_mean_err = []
var rk_mean_err = []
var steps = []

for (step = 0.01; step <= 1; step += 0.01){
var eu = euler(start, end, 1, 0 , step)
var rk = runge_kutta(start, end, 1, 0 , step)
var re = real(start, end, step)

eu_end_err.push(eu[eu.length-1] - re[re.length-1])
rk_end_err.push(rk[rk.length-1] - re[re.length-1])

var s1 = 0
var s2 = 0
for (k = 0; k < eu.length; k++){
s1 += Math.abs(eu[k] - re[k])
s2 += Math.abs(rk[k] - re[k])
}

eu_mean_err.push(s1 / eu.length)
rk_mean_err.push(s2 / eu.length)

steps.push(step)
}

.html('<script>require(["https://cdn.plot.ly/plotly-latest.min.js"], function(Plotly) { if (Plotly) window.Plotly = Plotly; })</script>');

var Plot = require('plotly-notebook-js');
var NotebookPlot = Plot.createPlot().constructor;
NotebookPlot.prototype._toHtml = NotebookPlot.prototype.render;

Out[31]:
In [32]:
Plot.createPlot([{ x: steps, y: eu_end_err, name: 'Euler method' },
{ x: steps, y: rk_end_err, name: 'Runge-Kutta method'}],
{title: 'error at end point VS. step size',
xaxis: {title: 'step size'},
yaxis: {title: 'error'}});

Out[32]:
In [33]:
Plot.createPlot([{ x: steps, y: eu_mean_err, name: 'Euler method' },
{ x: steps, y: rk_mean_err, name: 'Runge-Kutta method'}],
{title: 'mean local error VS. step size',
xaxis: {title: 'step size'},
yaxis: {title: 'error'}});

Out[33]:

## Problem 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 [ ]:
function single_step_rk(init1, init2, step){
var y = init1
var z = init2

var k_1 = step*z
var l_1 = -step*y
var k_2 = step*(z + 1/2*l_1)
var l_2 = -step*(y + 1/2*k_1)
var k_3 = step*(z + 1/2*l_2)
var l_3 = -step*(y + 1/2*k_2)
var k_4 = step*(z + l_3)
var l_4 = -step*(y + k_3)

var new_y = y + 1/6*(k_1 + 2*k_2 + 2*k_3 + k_4)
var new_z = z + 1/6*(l_1 + 2*l_2 + 2*l_3 + l_4)

return [new_y, new_z]
}

function adaptive_runge_kutta(start, end, init1, init2, step, err){
var y = [init1]
var z = [init2]
var steps = []

for (i = start+step; i <= end; i += step){
var res = single_step_rk(y[y.length-1], z[z.length-1], step)

var val_res = single_step_rk(y[y.length-1], z[z.length-1], step/2)
val_res = single_step_rk(val_res[0], val_res[1], step/2)

if (Math.abs(res[0] - val_res[0]) < err){
y.push(res[0])
z.push(res[1])
steps.push(step)
step *= 1.2
} else {
step = Math.max(0.001, step / 1.5)
i -= step
}
}

return [y, steps]
}

In [30]:
var init_step = 0.1
var mean_step = []
var errors = []

for (error = 0.000001; error <= .0001; error += 0.00001){
var ark = adaptive_runge_kutta(start, end, 1, 0 , init_step, error)
var steps = ark[1]

var s = 0
for (k = 0; k < steps.length; k++){
s += steps[k]
}

mean_step.push(s / steps.length)
errors.push(error)
}

Plot.createPlot([{ x: errors, y: mean_step }],
{title: 'mean step size VS. required error',
xaxis: {title: 'required error'},
yaxis: {title: 'mean step size'}});

Out[30]:
In [ ]: