2
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()
當我運行的代碼,它會導致與產生零誤差在DT的許多價值觀,具有正誤差之間的尖峯鋸齒情節龍格 - 庫塔四階誤差逼近。 (抱歉,我無法嵌入圖表的圖像)。這些大的尖峯不應該發生,因爲隨着時間步長dt的減小,誤差應該會持續下降。但是,我不知道如何解決這個問題,也不知道錯誤來自哪裏。我嘗試通過在while循環中添加尖峯消除尖峯,希望只有在dt時的誤差大於dt + delta時的誤差才能使其添加到我的dtime數組中,但它導致了完全相同的圖形。對於不同的時間步長
哇,這只是一個簡單的解決。謝謝! – infinitylord