Код Рунге-Кутта не сходится со встроенным методом

Question

Я пытаюсь реализовать метод runge-kutta для решения системы Lotka-Volterra, но код (ниже) работает неправильно. Я следовал рекомендациям, которые я нашел в других темах StackOverflow, но результаты не сходятся со встроенным методом Runge-Kutta, например, с использованием метода rk4, доступного в Pylab. Кто-то может мне помочь?

import matplotlib.pyplot as plt 
import numpy as np 
from pylab import * 

def meurk4(f, x0, t): 
    n = len(t) 
    x = np.array([ x0 ] * n)  

    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 * h * k1, t[i] + 0.5 * h) 
     k3 = h * f(x[i] + 0.5 * h * k2, t[i] + 0.5 * h) 
     k4 = h * f(x[i] + h * k3, t[i] + h) 

     x[i+1] = x[i] + (k1 + 2 * (k2 + k3) + k4) * 6**-1 

    return x 

def model(state,t): 

    x,y = state  

    a = 0.8 
    b = 0.02 
    c = 0.2 
    d = 0.004 
    k = 600 

    return np.array([ x*(a*(1-x*k**-1)-b*y) , -y*(c - d*x) ]) # corresponds to [dx/dt, dy/dt] 

# initial conditions for the system 
x0 = 500 
y0 = 200 

# vector of time 
t = np.linspace(0, 50, 100) 

result = meurk4(model, [x0,y0], t) 
print result 

plt.plot(t,result) 

plt.xlabel('Time') 
plt.ylabel('Population Size') 
plt.legend(('x (prey)','y (predator)')) 
plt.title('Lotka-Volterra Model') 
plt.show() 

---

Я просто обновил код следующие замечания. Таким образом, функция meurk4:

def meurk4(f, x0, t): 
     n = len(t) 
     x = np.array([ x0 ] * n)  

     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 * h * k1, t[i] + 0.5 * h) 
      k3 = h * f(x[i] + 0.5 * h * k2, t[i] + 0.5 * h) 
      k4 = h * f(x[i] + h * k3, t[i] + h) 

      x[i+1] = x[i] + (k1 + 2 * (k2 + k3) + k4) * 6**-1 

     return x 

теперь становится (исправлено):

def meurk4(f, x0, t): 
    n = len(t) 
    x = np.array([ x0 ] * n)  

    for i in range(n - 1): 

     h = t[i+1] - t[i] 

     k1 = f(x[i], t[i]) 
     k2 = f(x[i] + 0.5 * h * k1, t[i] + 0.5 * h) 
     k3 = f(x[i] + 0.5 * h * k2, t[i] + 0.5 * h) 
     k4 = f(x[i] + h * k3, t[i] + h) 

     x[i+1] = x[i] + (k1 + 2 * (k2 + k3) + k4) * (h/6) 

    return x 

Тем не менее, результаты заключается в следующем:

enter image description here

В то время как метод buitin RK4 (от Pylab):

enter image description here

Так что, конечно, мой код по-прежнему неверен, поскольку его результаты не совпадают с встроенным методом rk4. Пожалуйста, кто-нибудь может мне помочь?

Answer 1

Вы делаете очень типичную ошибку, смотри, например, How to pass a hard coded differential equation through Runge-Kutta 4 или здесь Error in RK4 algorithm in Python

Это либо

k2 = f(x+0.5*h*k1, t+0.5*h) 
... 
x[i+1]=x[i]+(k1+2*(k2+k3)+k4)*(h/6) 

или

k2 = h*f(x+0.5*k1, t+0.5*h) 

и так далее, с x[i+1], как это было , но не оба варианта одновременно.

---

Обновление: Более коварная ошибка выведенный тип начальных значений и, как следствие этого массива x векторов. По первоначальному определению оба являются целыми числами, и поэтому

x = np.array([ x0 ] * n)  

создает список целых векторов. Таким образом, шаг обновления

x[i+1] = x[i] + (k1 + 2 * (k2 + k3) + k4) * (h/6) 

всегда будет округлено до целого. А так как существует фаза, где оба значения падают ниже 1, интеграция стабилизируется в нуле. Таким образом, до

# initial conditions for the system 
x0 = 500.0 
y0 = 200.0 

, чтобы избежать этой проблемы.