Ich versuche, die Runge-Kutta-Methode zu implementieren, um ein Lotka-Volterra-System zu lösen, aber der Code (siehe unten) funktioniert nicht richtig. Ich habe die Empfehlungen befolgt, die ich in anderen Themen des StackOverflow gefunden habe, aber die Ergebnisse stimmen nicht mit der integrierten Runge-Kutta-Methode überein, wie zum Beispiel der in Pylab verfügbaren rk4-Methode. Könnte mir jemand helfen?

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

Ich habe gerade den Code nach den Kommentaren aktualisiert. Also, die Funktionmeurk4:

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

Jetzt werden (korrigiert):

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

Trotzdem lautet das Ergebnis wie folgt:

Bildbeschreibung hier eingeben

Während die buitin-Methode rk4 (von Pylab) Folgendes ergibt:

Bildbeschreibung hier eingeben

Also, sicherlich ist mein Code immer noch nicht korrekt, da seine Ergebnisse nicht mit denen der in rk4 eingebauten Methode übereinstimmen. Kann mir bitte jemand helfen?