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, надеясь, что он добавит только точки к моему массиву dtime, если ошибка в dt больше, чем ошибка в dt + delta, но это привело к точному тому же графику.Приближение ошибки 4-го порядка Рунге-Кутта для различных временных этапов
Вау, это просто легко исправить. Спасибо! – infinitylord