# Chapter 7 - Finite Differences: Ordinary Differential Equations¶

In [14]:
# Enabling interactive ploting and other pre-ample
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]

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"))
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. )

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)