In [2]:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.font_manager import FontProperties
def plot_approx(approximations,h_list,l):
font = {'family': 'monospace',
'color': '#39393d',
'weight': 'light',
'size': 14, }
fig = plt.figure(figsize=(10,6), facecolor="white")
ax1 = plt.subplot(221)
ax2 = plt.subplot(222)
ax3 = plt.subplot(223)
ax4 = plt.subplot(224)
# generating time and amplitude np.arrays based on differing time steps
t_varying_h = []
for i in range(len(h_list)):
h = h_list[i]
n = int(l/h)
time = []
for i in range(n):
time.append(i*h)
t_varying_h.append(time)
y_varying_h = approximations
t_varying_h = np.array(t_varying_h)
subplot(ax1,t_varying_h[0],y_varying_h[0][:,0],h_list[0])
subplot(ax2,t_varying_h[1],y_varying_h[1][:,0],h_list[1])
subplot(ax3,t_varying_h[2],y_varying_h[2][:,0],h_list[2])
subplot(ax4,t_varying_h[3],y_varying_h[3][:,0],h_list[3])
plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=1.0)
plt.show()
def subplot(ax,x,y,h,fontsize=12):
title = "Euler approx for step size = " + str(h)
ax.scatter(x,y,s=0.45,color="#004da0")
#ax.locator_params(nbins=3)
ax.set_xlabel('Time', fontsize=fontsize)
ax.set_ylabel('Amplitude', fontsize=fontsize)
ax.set_title(title, fontsize=fontsize)
In [3]:
##### Approximating amplitude of a simple harmonic oscillator using the Euler and Runga Kutta method
y0 = np.array([1,0])
h_list = [0.001,0.005,0.01,0.05]
l = 100*np.pi
def integrate_y0(y0):
return np.array([y0[1],-y0[0]])
def euler(y0,h,n):
approx_y_vectors = []
approx_y_vectors.append(y0)
for i in range(n-1):
y1 = y0 + h*integrate_y0(y0)
approx_y_vectors.append(y1)
y0 = y1
return np.array(approx_y_vectors)
def sample_euler(y0,h_list,l):
euler_samples = []
for i in range(len(h_list)):
h = h_list[i]
n = int(l/h)
approx = euler(y0,h,n)
euler_samples.append(approx)
return np.array(euler_samples)
euler_samples = sample_euler(y0,h_list,l)
plot_approx(euler_samples,h_list,l)
In [4]:
y0 = np.array([1,0])
def runga_kutta(y0,h,n):
approx_y_vectors = []
approx_y_vectors.append(y0)
for i in range(n-1):
k1 = h * integrate_y0(y0)
k2 = h * integrate_y0(np.array([k1[0]+h/2,y0[1]+k1[1]/2]))
k3 = h * integrate_y0(np.array([k2[0]+h/2,y0[1]+k2[1]/2]))
k4 = h * integrate_y0(np.array([k3[0]+h,y0[1]+k3[1]]))
y1 = y0 + k1/6 + k2/3 + k3/3 + k4/6
approx_y_vectors.append(y1)
y0 = y1
return array(approx_y_vectors)
def sample_runga_kutta(y0,h_list,l):
runga_kutta_samples = []
for i in range(len(h_list)):
h = h_list[i]
n = int(l/h)
approx = runga_kutta(y0,h,n)
runga_kutta_samples.append(approx)
return np.array(runga_kutta_samples)
runga_kutta_samples = sample_runga_kutta(y0,h_list,l)
plot_approx(runga_kutta_samples,h_list,l)
In [220]:
### calculating the errors for h = 0.005
def actual_function(h,n):
actual_values = []
for i in range(n):
value = np.cos(i*h)
actual_values.append(value)
return array(actual_values)
def compare_euler_runga_kutta(y0,h,l):
n = int(l/h)
time = []
euler_values = euler(y0,h,n)[:,0]
runga_kutta_values = runga_kutta(y0,h,n)[:,0]
actual_values = actual_function(h,n)
errors = []
###############Â check that you're extracting the correct element from the euler array
for i in range(n):
euler_error = actual_values[i] - euler_values[i]
runga_kutta_error = actual_values[i] - runga_kutta_values[i]
errors.append([euler_error,runga_kutta_error])
return np.array(errors)
def plot_errors(errors,h,l):
fig = plt.figure(figsize=(10,6), facecolor="white")
ax = fig.add_subplot(111)
n = int(l/h)
time = []
for i in range(n):
time.append(i*h)
# plot euler error
euler_line, = plt.plot(time,errors[:,0],label='Euler error')
# plot runga kutta error
runga_line, = plt.plot(time,errors[:,1], label="Runga Kutaa error")
plt.legend(handles=[euler_line,runga_line])
plt.show()
def average_error_for_different_step_sizes(h_list,l,y0):
average_errors = []
for i in range(0,10):
h = 0.001*(i+1)
errors = compare_euler_runga_kutta(y0,h,l)
n = int(l/h)
euler_sum = np.sum(errors[:,0])
runga_sum = np.sum(errors[:,1])
average_errors.append([euler_sum,runga_sum])
print np.array(average_errors)
h = 0.005
#errors = compare_euler_runga_kutta(y0,h,l)
#average_error_for_different_step_sizes(0.005,l,y0)
plot_errors(errors,h,l)
In [ ]:
# determining errors from the step size
def error_from_step_size():
approx_errors = []
for i in range(0,50):
h = (i+1)*0.001
euler_error_h = euler(y0,h,10)
euler_error_half_h = euler(y0,h/2,10)
euler_error = euler_error_h[:,0] - euler_error_half_h[:,0]
#runga_error_h = runga_kutta(y0,h,10)
#runga_error_half_h = runga_kutta(y0,h/2,10)
#print euler_error
return(array(approx_errors))
def plot_error_vs_h():
errors = error_from_step_size()
fig = plt.figure(figsize=(10,6), facecolor="white")
ax = fig.add_subplot(111)
# plot euler error
euler_line, = plt.plot(errors[:,0],errors[:,1],label='Euler error')
# plot runga kutta error
runga_line, = plt.plot(errors[:,0],errors[:,2], label="Runga Kutaa error")
plt.legend(handles=[euler_line,runga_line])
plt.show()
#plot_error_vs_h()
error_from_step_size()
In [ ]: