Runge-Kutta代码未与内置方法收敛 [英] Runge-Kutta code not converging with builtin method

问题描述

我正在尝试实现runge-kutta方法来解决Lotka-Volterra系统，但是代码(波纹管)无法正常工作.我遵循了在StackOverflow的其他主题中找到的建议，但是结果与内置的Runge-Kutta方法(例如Pylab中的rk4方法)不一致.有人可以帮助我吗?

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

---

我刚刚按照评论更新了代码.因此，函数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

尽管如此，结果如下:

在此处输入图像描述

当buitin方法rk4(来自Pylab)得出以下结果时:

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

在此处输入图片描述

因此，当然我的代码仍然不正确，因为其结果与内置的rk4方法不同.拜托，有人可以帮我吗?

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?

推荐答案

您正在做一个非常典型的错误，例如，参见 Python中RK4算法中的错误

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

是

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]不变，但不能同时使用两个变体.

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

更新:一个更隐蔽的错误是推断的初始值类型以及x向量数组的结果.按照原始定义，它们都是整数，因此

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)

将始终四舍五入为整数.并且由于存在两个值均低于1的相位，因此积分稳定在零.因此修改为

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

为避免该问题.

这篇关于Runge-Kutta代码未与内置方法收敛的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！