infinitylord - 8 months ago 80

Python Question

`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.

Answer

A short test demonstrates

```
In [104]: import numpy as np
In [105]: dt = 0.6
In [106]: np.arange(0, 1+dt, dt)
Out[106]: array([ 0. , 0.6, 1.2])
```

Thus to get meaningful results, either set `t[n-1]=1`

at the start or compute the error as

```
E = abs(x[n-1]-exact(t[n-1]))
```

Source (Stackoverflow)