# 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>'
))
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.
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.
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"))
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"))
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.
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"))
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)
# Frequency analysis for amplitute 1
freqs = [1, 2, 5, 10]
results = [evaluate_adaptive_2(runge_kutte_step_param_osc, 1., 1e-3, omega=2*np.pi*f) for f in freqs]
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"))
data = [Scatter(x=t, y=trajectory, mode='lines', name='f= '+str(n)) for (_, trajectory, t), n in zip(results,freqs)]
layout = Layout()
fig = Figure(data=data, layout=layout)
iplot(fig, show_link=False)
You can actually run javascript inside Jupyter with a Python kernel, just by using the %%javascript magic command. Then, you can use the following hack to retrieve your variable in python! THIS IS SO COOL
%%javascript
var kernel = IPython.notebook.kernel;
function benchmark(terms, threads) {
//
// Split pi computation by number of workers
//
var inds = [];
for (var k = 1; k < terms; k += (terms/threads)){
inds.push(k);}
inds.push(terms);
//
// Spawn workers
//
var finished = [];
var sum_pi = 0;
var st_tm = Date.now();
for (var i = 0; i < threads; i++) {
var blob = new Blob(['(' + worker.toString() + '())']);
var url = window.URL.createObjectURL(blob);
var webworker = new Worker(url);
webworker.addEventListener('message', function (evt) {
var pi = evt.data.pi;
finished.push('done');
sum_pi += pi;
if (finished.length == threads){
kernel.execute("val = " + sum_pi.toFixed(6));
var sec = (Date.now()-st_tm)/1000;
kernel.execute("sec = " + sec);
var mflops = 5 * terms / (sec * 1e6);
kernel.execute("mflops = " + mflops.toFixed(6));
sum_pi = 0;
finished = [];
}
this.terminate()
});
webworker.postMessage({start: inds[i],end: inds[i+1]})
}
window.URL.revokeObjectURL(url);
}
function worker() {
self.addEventListener('message', function (evt) {
var start = evt.data.start;
var end = evt.data.end;
var pi = 0;
var tstart = Date.now();
for (var term = start; term < end; ++term)
pi += 0.5 / ((term - 0.75) * (term - 0.25))
var tend = Date.now();
var dt = tend - tstart;
self.postMessage({pi: pi, dt: dt})
})
}
benchmark(1e8, 4)
print("Pi: {0}, Time: {1}, MFlops: {2}".format(val, sec, mflops))
print("This calculation happened in Javascript")