from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
def f(x, t): #function
return -x
def exact(t): #exact solution
return np.exp(-t)
def Rk4(x0, t0, dt): #Runge-Kutta Fourth Order Error
t = np.arange(0, 1+dt, dt)
n = len(t)
x = np.array([x0]*n)
x[0],t[0] = x0,t0
for i in range(n-1):
h = t[i+1]-t[i]
k1 = h*f(x[i], t[i])
k2 = h*f(x[i]+0.5*k1, t[i]+0.5*h)
k3 = h*f(x[i]+0.5*k2, t[i]+0.5*h)
k4 = h*f(x[i]+k3, t[i+1])
x[i+1] = x[i]+(k1+2.0*(k2+k3)+k4 )/6.0
E = abs(x[n-1]-exact(1))
return E
vecRk4 = np.vectorize(Rk4)
dt = 10e-4
dtime = []
delta = 10e-4
while dt < 1:
if Rk4(1.0,0.0,dt) < Rk4(1.0,0.0,dt+delta):
dtime.append(dt)
S = vecRk4(1.0,0.0,dtime)
dt = dt + delta
plt.plot(dtime,S)
plt.xlabel("dt (s)")
plt.ylabel("Error")
plt.show()
When I run the code, it results in a jagged plot with spikes that yield zero error at many values of dt, with positive error in-between. (sorry, I can't embed an image of the graph). These large spikes should not be occurring, as there should be a continuous decrease in error as the time-step dt is decreased. However, I'm not sure how to fix this nor do I know where the error results from. I tried eliminating the spikes by adding in the while loop, hoping to have it only add points to my dtime array if the error at dt is larger than the error at dt+delta, but it resulted in precisely the same graph.
