Runge-Kutta code not converging with builtin method

Question

I am trying to implement the runge-kutta method to solve a Lotka-Volterra systtem, but the code (bellow) is not working properly. I followed the recomendations that I found in other topics of the StackOverflow, but the results do not converge with the builtin Runge-Kutta method, like rk4 method available in Pylab, for example. Someone could help me?

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

---

I just updated the code following the comments. So, the function 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

Becomes now (corrected):

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

Nevertheless, the results is the following:

While the buitin method rk4 (from Pylab) results the following:

So, certainly my code still is not correct, as its results are not the same of the builtin rk4 method. Please, someone can help me?

Inconsonant

with that modification, variant 1, you need to modify the factor in x[i+1] to (h/6). – Vervain

You are right, my mistake. Nevertheless, I just implemented this, and it still discrepant in relation to the Pylab builtin method. – Inconsonant

Are you absolutely sure that the constants and initial values are the same? – Vervain

I just changed result = meurk4( model, [x0,y0], t ) by result = rk4( model, [x0,y0], t ) on the code. – Inconsonant

One can only guess at the internals of pylab/matplotlib rk4. I'd assume that it has some kind of step size control. A step size of 0.5 is too large. 0.01 or 0.001 would be appropriate. But does not seem to help. – Vervain

I am a little bit frustrated, but I'll keep trying here. Anything else, I'll post here. Thanks LutzL for you contribuitions with my question. – Inconsonant

See update, unsuspected type inference often leads to strange numerical errors in python. – Vervain

Answer 1

You are doing a very typical error,see for instance How to pass a hard coded differential equation through Runge-Kutta 4 or here Error in RK4 algorithm in Python

It is either

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

or

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

and so on, with x[i+1] as it was, but not both variants at the same time.

---

Update: A more insidious error is the inferred type of the initial values and in consequence of the array of x vectors. By the original definition, both are integers, and thus

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

creates a list of integer vectors. Thus the update step

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

will always round to integer. And since there is a phase where both values fall below 1, the integration stabilizes at zero. Thus modify to

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

to avoid that problem.

Vervain

Thanks for correcting my mistake about using x.* in numpy, as you suggested I was thinking of Matlab. Your answer is right anyways, my brain is fried from a long day of work. – Deepen

I just removed the h factor before the k, i.e., variant 2, and the program runs. Even if it converges to zero. (Perhaps choose the parameters so that the equilibrium point is in the first quadrant) -- Please add a section to the question on what you corrected and a description of the new error. – Vervain

Thanks LutzL. I just updated the question with a new section with the code modified according with your sugestions. – Inconsonant