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?

\begin{equation} y(x+h) = y(x) + \frac{h}{2}[f(x,y(x)) + f(x+h, y(x) + hf(x,y(x)))] \end{equation}

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{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} \end{equation}

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.

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;

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'}});

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'}});

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.

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]
}

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'}});