我的 rk4 实现的加权平均值中的不同近似值都是相同的

The distinct approximatons in the weighted average for my rk4 implementation are all the same

好的,所以我正在尝试编写一个四阶龙格库塔方法的实现来实现微分方程的数值逼近,这只是我数学课程的一部分,但我也想学习一些编程,但问题是它使用这些系数在每个近似步骤中。所以我们从 x 值和 y 值开始,它想为后面的 x 值找到 y 值,我给它一个步长,通常是 0.1,它在 x 中向上移动 0.1,并给了我一个新的近似值y 值,并且在这些步骤中的每一步中,它都执行 4 个近似值,我称之为 k1、k2、k3 和 k4。然后它采用这 4 个近似值的加权平均值,并给我最终的近似值,这应该是非常准确的。问题是对于每一步我都得到 k1=k2=k3,但是 k4 是不同的。这似乎不对。我在一个函数上测试它,给它 f(x,y)=2x-3y+1,使用步长 h,我的 k 如下:

k1=f(x,y)
k2=f(x+.5h,y+.5hk1)
k3=f(x+.5h,y+.5hk2)
k4=f(x+h,y+hk3)

所以只是用代数方法检查这个有 k1=k2 我最终得到 9y=1-6x,但我给它的初始值是 (x0,y0)=(1,5),显然是 45 !=1-6

所以这一切似乎都是错误的。根据教科书我没有得到正确的答案,而且这些 k 真的不应该相同。所以无论如何这是我的代码。 k 之间的打印语句只是作为测试,以查看 k 在我们循环 n 时实际上在变化。

import numpy

def rk4(x0,y0,xf,h,f):
    y=[]
    x=numpy.linspace(x0,xf,(xf-x0)/h+1)
    y.insert(0,y0)
    for n in range(len(x)-1):
        k1=f(x[n],y[n])
        print k1
        k2=f(x[n]+(1/2)*h,y[n]+(1/2)*h*k1)
        print k2
        k3=f(x[n]+(1/2)*h,y[n]+(1/2)*h*k2)
        print k3
        k4=f(x[n]+h,y[n]+h*k3)
        print k4
        y.insert(n+1,y[n]+(h/6)*(k1+2*k2+2*k3+k4))
        print y[n]+(h/6)*(k1+2*k2+2*k3+k4)

    print x
    print y

def twoxminusthreeyplusone(x,y):
    return 2*x-3*y+1

rk4(1,5,1.5,0.1,twoxminusthreeyplusone)

然后我得到这个输出

/usr/bin/python2.7 /home/t/PycharmProjects/untitled/chunk.py
-12.0
-12.0
-12.0
-8.2
3.86333333333
-8.39
-8.39
-8.39
-5.673
3.06961666667
-5.80885
-5.80885
-5.80885
-3.866195
2.52110925
-3.96332775
-3.96332775
-3.96332775
-2.574329425
2.14792644708
-2.64377934125
-2.64377934125
-2.64377934125
-1.65064553888
1.900100743
[ 1.   1.1  1.2  1.3  1.4  1.5]
[5, 3.8633333333333333, 3.0696166666666667, 2.5211092500000003, 2.1479264470833335, 1.9001007429979169]

Process finished with exit code 0

所以我不明白。 k 相同似乎是主要问题。我已经尝试了几种不同的方法,就像之前一样,我让每个 k 都像 y 一样作为一个空列表开始,而不是仅仅将 k 定义为 x[n],y[n] 的函数,并在使用 .insert(n,x) 或 .insert(n-1,x) 循环遍历 n 取决于我如何定义每个步骤,但这似乎不必要地复杂,我认为实际上最终遇到了同样的问题

您在通常意想不到的地方遇到了典型问题。

将所有因素 (1/2) 替换为 0.5(或将 (1/2)*h 替换为 h/2,但尽可能避免除法)。整数除法结果为 0,从而导致观察到的行为。


但请注意,对于 python 3,您的代码可以正常工作,但必须强制执行整数行为。