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

ответ

1

Короткий тест демонстрирует

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]) 

Таким образом, чтобы получить значимые результаты, либо установить t[n-1]=1 в начале или вычислить ошибку как

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

Вау, это просто легко исправить. Спасибо! – infinitylord

 Смежные вопросы

  • Нет связанных вопросов^_^