Chapter 7 - Finite Differences: Ordinary Differential Equations

In [14]:
# Enabling interactive ploting and other pre-ample
from plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot
from IPython.display import display, HTML
from plotly.graph_objs import *
init_notebook_mode(connected=True)
display(HTML(
    '<script>'
        'var waitForPlotly = setInterval( function() {'
            'if( typeof(window.Plotly) !== "undefined" ){'
                'MathJax.Hub.Config({ SVG: { font: "STIX-Web" }, displayAlign: "center" });'
                'MathJax.Hub.Queue(["setRenderer", MathJax.Hub, "SVG"]);'
                'clearInterval(waitForPlotly);'
            '}}, 250 );'
    '</script>'
))

(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?

The Heun method equation is:   $y(x+h) = y(x) + \frac{h}{2}\big\{f(x,y(x)) + f[x+h, y(x) + hf(x,y(x))]\big\}$

For bookkeeping, let's substitute with the variables $u = x+h$ and $v = y(x) + hf(x,y(x))$. We will need their derivatives with respect to $h$, which are the following $$ \begin{align} \frac{\partial u}{\partial h} &= 1 \\ \frac{\partial v}{\partial h} &= f \end{align} $$

Now, let's expand the slope average in the RHS of the first equation (Taylor expansion as a function of $h$): $$\begin{align} y(x+h) &= y(x) + \frac{h}{2}\big\{f(x,y(x)) + f(u, v)\big\} \\ &= y(x) + \frac{h}{2}\big\{f(x,y(x)) + f(x,y(x)) + h\frac{d}{dh}f(u, v) + \mathcal{O}(h^2)\big\} \\ &= y(x) + \frac{h}{2}\big\{2f(x,y(x)) + h\big(\frac{\partial f}{\partial u}\frac{\partial u}{\partial h} + \frac{\partial f}{\partial v}\frac{\partial v}{\partial h}\big) + \mathcal{O}(h^2) \big\} \\ &= y(x) + \frac{h}{2}\big\{2f(x,y(x)) + h\big[\frac{\partial f}{\partial x} + \frac{\partial f}{\partial h}f + \mathcal{O}(h^2) \big]\big\} \\ &= y(x) + hf(x,y(x)) + \frac{h^2}{2}\big[ \frac{\partial f}{\partial x} + \frac{\partial f}{\partial h}f \big] + \mathcal{O}(h^3) \end{align}$$

Conclusively, this method is a third-order approximation, and thus has no second-order approximation error.

(7.2) For a simple harmonic oscillator $\ddot{y}+\dot{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.

For the simple harmonic oscillator, we have the following system of equations: $$\begin{align} \frac{dy}{dt} = f && \frac{df}{dt} = -y \end{align}$$

thus to find $y(t)$ with an Euler method, we need to iteratively solve $$\begin{align} y(t+h) &= y(t) + h\frac{dy}{dt} = y(t) + hf\\ f(t+h) &= f(t) + h\frac{df}{dt} = f(t) -hy \end{align}$$

and respectively for a fourth-order Runge-Kutte method $$\begin{align} y(t+h) &= y(t) + h\big(\frac{k_{1,y} + 2k_{2,y} + 2k_{3,y} + k_{4,y}}{6}\big) \\ f(t+h) &= f(t) + h\big(\frac{k_{1,f} + 2k_{2,f} + 2k_{3,f} + k_{4,f}}{6}\big) \end{align}$$

with (performing half-stepping trick recursively) $$\begin{align} k_{1,y} &= f(t) && k_{1,f} = \dot{f}(t) = -y(t) \\ k_{2,y} &= f(t) + \frac{h}{2}k_{1,f} && k_{2,f} = -y(t) + \frac{h}{2}k_{1,y} \\ k_{3,y} &= f(t) + \frac{h}{2}k_{2,f} && k_{2,f} = -y(t) + \frac{h}{2}k_{2,y} \\ k_{4,y} &= f(t) + hk_{3,f} && k_{2,f} = -y(t) + hk_{3,y} \\ \end{align}$$

To calculate the errors, we know that the actual solution to the differential equation is $y(t)=cos(t).$ We use Python to solve iteratively the above equations and plot the error for $5$ different step sizes.

In [131]:
import numpy as np
import time
from tabulate import tabulate

def euler_step(h, y, f):
    return (y+h*f, f-h*y)

def runge_kutte_step(h, y, f):
    ky=np.zeros(4); kf = np.zeros(4)
    ky[0] = h*f
    kf[0] = -h*y
    ky[1] = h*(f + kf[0]/2.0)
    kf[1] = h*(-y - ky[0]/2.0)
    ky[2] = h*(f + kf[1]/2.0)
    kf[2] = h*(-y - ky[1]/2.0)
    ky[3] = h*(f + kf[2])
    kf[3] = h*(-y - ky[2])
    return (y + np.dot(ky,[1.,2.,2.,1.])/6., f + np.dot(kf,[1.,2.,2.,1.])/6. )

def evaluate(step_method, step_sizes):
    start = time.time()
    avg_err = np.zeros_like(step_sizes)
    avg_loc_err = np.zeros_like(step_sizes)
    val_err = np.zeros_like(step_sizes)
    slope_err = np.zeros_like(step_sizes)
    durations = np.zeros_like(step_sizes)
    for i, h in enumerate(step_sizes):
        t = 0
        steps = int(100*np.pi/h)
        err = 0
        loc_err = 0
        y = 1
        f = 0
        for s in range(steps):
            (y_full,f_full) = step_method(h, y, f)
            (y_half, f_half) = step_method(h/2., y, f)
            (y_2half, f_2half) = step_method(h/2., y_half, f_half)
            err += np.abs(np.cos(t) - y_full)
            loc_err += np.abs(y_full - y_2half)
            t += h
            y = y_full
            f = f_full
        avg_err[i] = err / steps
        avg_loc_err[i] = loc_err / steps
        val_err[i] = np.abs(np.cos(t) - y)
        slope_err[i] = np.abs(np.sin(t) - f)
        durations[i] = time.time() - start
    return np.array([step_sizes, avg_loc_err, avg_err, val_err, slope_err, durations]).T.tolist() 

step_sizes = [1., .1, .01, .001, .0001]
euler_results = evaluate(euler_step, step_sizes)
print("Euler method results:")
print(tabulate(euler_results, headers=["Step", "Avg. Local Error", "Avg. Error", "Final Value Error", "Slope Error", "Duration"], 
                     tablefmt="rst", numalign="center", stralign="center", floatfmt=".2e"))
runge_results = evaluate(runge_kutte_step, step_sizes)
print("Runge Kutte method results:")
print(tabulate(runge_results, headers=["Step", "Avg. Local Error", "Avg. Error", "Final Value Error", "Slope Error", "Duration"], 
                     tablefmt="rst", numalign="center", stralign="center", floatfmt=".2e"))
Euler method results:
========  ==================  ============  ===================  =============  ==========
  Step     Avg. Local Error    Avg. Error    Final Value Error    Slope Error    Duration
========  ==================  ============  ===================  =============  ==========
1.00e+00       2.42e+44         9.70e+44         9.87e-01          1.83e+47      2.28e-03
1.00e-01       6.13e+02         2.46e+05         2.78e+06          5.45e+06      1.45e-02
1.00e-02       3.86e-05         9.08e-01         3.81e+00          1.04e-01      1.36e-01
1.00e-03       1.72e-07         5.27e-02         1.70e-01          6.98e-04      1.27e+00
1.00e-04       1.60e-09         5.03e-03         1.58e-02          1.33e-04      1.28e+01
========  ==================  ============  ===================  =============  ==========
Runge Kutte method results:
========  ==================  ============  ===================  =============  ==========
  Step     Avg. Local Error    Avg. Error    Final Value Error    Slope Error    Duration
========  ==================  ============  ===================  =============  ==========
1.00e+00       2.20e-03         4.53e-01         1.04e+00          2.97e-01      1.08e-02
1.00e-01       4.97e-08         6.36e-02         3.72e-05          1.19e-01      1.15e-01
1.00e-02       4.97e-13         6.37e-03         4.59e-10          1.85e-02      1.18e+00
1.00e-03       2.24e-17         6.37e-04         1.39e-13          5.31e-04      1.17e+01
1.00e-04       2.16e-17         6.37e-05         4.75e-13          1.31e-04      1.22e+02
========  ==================  ============  ===================  =============  ==========

(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 [150]:
def evaluate_adaptive(step_method, initial_step, desired_error):
    start = time.time()
    h = initial_step
    h_avg = 0
    t = 0
    steps = 0
    err = 0
    loc_err = 0
    y = 1
    f = 0
    while t < 100*np.pi/h:
        (y_full,f_full) = step_method(h, y, f)
        (y_half, f_half) = step_method(h/2., y, f)
        (y_2half, f_2half) = step_method(h/2., y_half, f_half)
        steps += 1
        if np.abs(y_full - y_2half) > desired_error:  # If we are above the desired error
            h = h/np.sqrt(2)  # Reduce the step size
            continue
        y = y_2half
        f = f_2half
        loc_err += np.abs(y_full - y_2half)
        err += np.abs(np.cos(t) - y)
        t += h
        h_avg += h
        if np.abs(y_full - y_2half) < (desired_error/2.):  # If we are lower than the desired error
            h = h*np.sqrt(2)  # Increase the step size
            
    avg_err = err / steps
    avg_loc_err = loc_err / steps
    h_avg = h_avg / steps
    val_err = np.abs(np.cos(t) - y)
    slope_err = np.abs(np.sin(t) - f)
    duration = time.time() - start
    return desired_error, h_avg, avg_loc_err, avg_err, val_err, slope_err, duration

des_err = [1., .1, .01, .001, .0001]
adaptive_results = [evaluate_adaptive(runge_kutte_step, 1., err) for err in des_err]

print("Adaptive Runge Kutte method results, Initial step 1.:")
print(tabulate(adaptive_results, headers=["Desired Error", "Average Step", "Avg. Local Error", "Avg. Error", "Final Value Error", "Slope Error", "Duration (s)"], 
                     tablefmt="rst", numalign="center", stralign="center", floatfmt=".2g"))
adaptive_results = [evaluate_adaptive(runge_kutte_step, 0.1, err) for err in des_err]
print("Adaptive Runge Kutte method results, Initial step 0.1:")
print(tabulate(adaptive_results, headers=["Desired Error", "Average Step", "Avg. Local Error", "Avg. Error", "Final Value Error", "Slope Error", "Duration (s)"], 
                     tablefmt="rst", numalign="center", stralign="center", floatfmt=".2g"))
Adaptive Runge Kutte method results, Initial step 1.:
===============  ==============  ==================  ============  ===================  =============  ==============
 Desired Error    Average Step    Avg. Local Error    Avg. Error    Final Value Error    Slope Error    Duration (s)
===============  ==============  ==================  ============  ===================  =============  ==============
       1              2.2               0.22             0.77             0.35              0.91           0.0016
      0.1             0.89             0.022             0.43             0.38               1.6           0.015
     0.01             0.69             0.0033            0.4              0.052              1.5           0.016
     0.001            0.35            0.00023            0.21            0.0075               1            0.038
    0.0001            0.21            2.2e-05            0.13            0.0012             0.71           0.075
===============  ==============  ==================  ============  ===================  =============  ==============
Adaptive Runge Kutte method results, Initial step 0.1:
===============  ==============  ==================  ============  ===================  =============  ==============
 Desired Error    Average Step    Avg. Local Error    Avg. Error    Final Value Error    Slope Error    Duration (s)
===============  ==============  ==================  ============  ===================  =============  ==============
       1              1.7               0.2              0.45             0.61               0.8           0.0012
      0.1              1               0.029             0.5              0.21               1.6           0.0043
     0.01             0.66             0.0033            0.36             0.096              1.9           0.012
     0.001            0.32             0.0002            0.19             0.015              1.5           0.032
    0.0001            0.32            4.3e-05            0.2             0.0028              1.1           0.041
===============  ==============  ==================  ============  ===================  =============  ==============

(7.4) Numerically solve the differential equation $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 apmlitute and frequency of the excitation.

We will solve the system for $z(t) = A\cos{\omega t}$, assuming the following initial conditions (for $t=0$): $$\begin{align}z(0)=A && \dot{z}(0) = 0 && \theta(0)=0 && \dot{\theta}(0)=0\end {align} $$

We transform the second order differential equation into a system of first order equations: $$\begin{align} \frac{d\theta}{dt} &= f \\ \frac{df}{dt} &= \ddot{\theta} = -\frac{g-\omega^2A\cos{\omega t}}{l}\sin{\theta} \end{align}$$

The rest of the analytical form of a fourth order Runge-Kutta is demonstrated in Problem (7.2). We continue with the numerical evaluation using Python. We will use the adaptive Runge-Kutta as it runs very fast.

In [191]:
def runge_kutte_step_param_osc(h, t, y, f, g=9.8, l=10, A=1, omega=1):
    ky=np.zeros(4); kf = np.zeros(4)
    ky[0] = h*f
    kf[0] = h*(-(g-omega*omega*A*np.cos(omega*t))/l)*np.sin(y)
    ky[1] = h*(f + kf[0]/2.0)
    kf[1] = h*((-(g-omega*omega*A*np.cos(omega*t))/l)*np.sin(y) - ky[0]/2.0)
    ky[2] = h*(f + kf[1]/2.0)
    kf[2] = h*((-(g-omega*omega*A*np.cos(omega*t))/l)*np.sin(y) - ky[1]/2.0)
    ky[3] = h*(f + kf[2])
    kf[3] = h*((-(g-omega*omega*A*np.cos(omega*t))/l)*np.sin(y) - ky[2])
    return (y + np.dot(ky,[1.,2.,2.,1.])/6., f + np.dot(kf,[1.,2.,2.,1.])/6. )

def evaluate_adaptive_2(step_method, initial_step, desired_error, **options):
    start = time.time()
    h = initial_step
    h_avg = 0
    t = [0,]
    steps = 0
    err = 0
    loc_err = 0
    y = 1
    f = 0
    trajectory = [0,]
    while t[-1] < 50*np.pi/h:
        (y_full,f_full) = step_method(h, t[-1], y, f, **options)
        (y_half, f_half) = step_method(h/2., t[-1], y, f, **options)
        (y_2half, f_2half) = step_method(h/2., t[-1], y_half, f_half, **options)
        steps += 1
        if np.abs(y_full - y_2half) > desired_error:  # If we are above the desired error
            h = h/np.sqrt(2)  # Reduce the step size
            continue
        y = y_2half
        f = f_2half
        trajectory.append(y)
        loc_err += np.abs(y_full - y_2half)
        err += np.abs(np.cos(t[-1]) - y)
        t.append(t[-1] + h)
        h_avg += h
        if np.abs(y_full - y_2half) < (desired_error/2.):  # If we are lower than the desired error
            h = h*np.sqrt(2)  # Increase the step size
            
    avg_err = err / steps
    avg_loc_err = loc_err / steps
    h_avg = h_avg / steps
    val_err = np.abs(np.cos(t[-1]) - y)
    slope_err = np.abs(np.sin(t[-1]) - f)
    duration = time.time() - start
    return [desired_error, h_avg, avg_loc_err, avg_err, val_err, slope_err, duration], trajectory,t

# Amplitute analysis
amps = [1, 10, 20, 30]
results = [evaluate_adaptive_2(runge_kutte_step_param_osc, .001, 1e-5, A=a) for a in amps]
for adaptive_results, trajectory, t in results:
    print("Adaptive Runge Kutte method results, Initial step .01:")
    print(tabulate([adaptive_results], headers=["Desired Error", "Average Step", "Avg. Local Error", "Avg. Error", "Final Value Error", "Slope Error", "Duration (s)"], 
                         tablefmt="rst", numalign="center", stralign="center", floatfmt=".2g"))
Adaptive Runge Kutte method results, Initial step .01:
===============  ==============  ==================  ============  ===================  =============  ==============
 Desired Error    Average Step    Avg. Local Error    Avg. Error    Final Value Error    Slope Error    Duration (s)
===============  ==============  ==================  ============  ===================  =============  ==============
     1e-05           0.084            3.9e-06            0.52             0.49                1             0.38
===============  ==============  ==================  ============  ===================  =============  ==============
Adaptive Runge Kutte method results, Initial step .01:
===============  ==============  ==================  ============  ===================  =============  ==============
 Desired Error    Average Step    Avg. Local Error    Avg. Error    Final Value Error    Slope Error    Duration (s)
===============  ==============  ==================  ============  ===================  =============  ==============
     1e-05           0.043            4.1e-06            4.8               6.3              0.83            0.72
===============  ==============  ==================  ============  ===================  =============  ==============
Adaptive Runge Kutte method results, Initial step .01:
===============  ==============  ==================  ============  ===================  =============  ==============
 Desired Error    Average Step    Avg. Local Error    Avg. Error    Final Value Error    Slope Error    Duration (s)
===============  ==============  ==================  ============  ===================  =============  ==============
     1e-05           0.033            4.4e-06            9.3               13               0.85            1.3
===============  ==============  ==================  ============  ===================  =============  ==============
Adaptive Runge Kutte method results, Initial step .01:
===============  ==============  ==================  ============  ===================  =============  ==============
 Desired Error    Average Step    Avg. Local Error    Avg. Error    Final Value Error    Slope Error    Duration (s)
===============  ==============  ==================  ============  ===================  =============  ==============
     1e-05            0.03            4.3e-06             30               1.2              0.43            1.4
===============  ==============  ==================  ============  ===================  =============  ==============
In [195]:
data = [Scatter(x=t, y=trajectory, mode='lines', name='A= '+str(n)) for (_, trajectory, t), n in zip(results,amps)]
layout = Layout()
fig = Figure(data=data, layout=layout)
iplot(fig, show_link=False)