RK4法求解四阶IVP：BASIC转Python代码数值结果偏差问题

我正在将用于数值求解四阶初值问题（IVP）的Runge-Kutta 4阶（RK4）方法的BASIC代码转换为Python代码，该四阶IVP定义为：
$y^{(4)} = -4\sin t - y - y' - y'' + 4y''' + y'''(y''' - y'' + y' - y)$，其解析解为$y(t)=e^t+\sin t$。

我编写的Python代码 #

from math import sin as SIN, cos as COS, exp as EXP
A=[0]*5; G=[[0]*5]*5; W=[0]*5; X=[0]*5; W[1]=W[4]=1/6; W[2]=W[3]=1/3
def FNF(T, X1, X2, X3, X4): return -4*SIN(T)-X1-X2-X3-(4*X4)+X4*(X4-X3+X2-X1)
def FNQ(Z): return Z
print("t0,n,h,kmax"); T, N, H, KM = [float(_) for _ in input().split(" ")]
for I in range(1, int(N)): print(f"INPUT D{I}(X)"); X[int(N-I)] = float(input())
for K in range(1, int(KM)+1): 
    print(T, X[int(N)], EXP(T)+SIN(T))
    for I in range(1, 4+1):
        for J in range(1, int(N)+1):
            P = 0.5
            if I == 1: P = 0
            if I == 4: P = 1
            if J > 1: G[I][J] = FNQ(X[J - 1])
            X[J] += H*W[I]*G[I][J]
            A[J] = X[J] + P*H*G[I][J]
        G[I][1] = FNF(T + P*H, A[1], A[2], A[3], A[4])
    T += H

然而数值计算结果与解析结果始终存在约1的偏差。以下是经验证可用的BASIC代码：

经验证可用的BASIC代码 #

CLS:KEY OFF: "Program RK-n for nth order IVPs"
DIM A(4),G(4,4),W(4),X(4):W(1) = W(4) = 1/6:W(2)= W(3)=1/3
DEF FNF(T,X1,X2,X3,X4)= -4*SIN(T)- X1 - X2- X3 + 4*X4 + X4*(X4- X3 + X2 - X1)
DEF FNQ(Z) = Z
INPUT"t0,n,h,kmax";T,N,H,KM
FOR I=1 TO N-1: PRINT "INPUT D";I;"(X)":INPUT X(N -I):NEXT I
FOR K =1 TO KM
PRINT T;X(N);EXP(T) + SIN(T)
FOR I=1 TO 4: FOR J = 1 TO N
P= .5:IF I =1 THEN P= 0: IF I = 4 THEN P = 1
IF J > 1 THEN G(I,J) = FNQ(X(J-1))
X(J)=X(J) + H*W(I)*G(I,J)
A(J)=X(J)+ P*H*G(I,J)
NEXT J
G(1,1)=FNF(T+P*H,A(1),A(2),A(3),A(4))
NEXT I:T = T + H:NEXT K: END

问题排查与修正 #

对比两段代码，发现两处关键错误：

1. FNF函数符号错误
   BASIC中的FNF函数是-4*SIN(T)- X1 - X2- X3 + 4*X4 + ...，而Python代码中错误写成了-4*SIN(T)-X1-X2-X3-(4*X4)+...，将+4*X4误写为-4*X4，这是结果偏差的核心原因。

2. 二维数组初始化引用共享问题
   Python中G=[[0]*5]*5的写法会创建5个指向同一列表的引用，修改任意子列表都会影响所有子列表。正确的初始化方式应为G = [[0]*5 for _ in range(5)]，避免引用共享。

修正后的Python代码 #

from math import sin as SIN, cos as COS, exp as EXP
# 正确初始化二维数组，避免引用共享
A = [0] * 5
G = [[0] * 5 for _ in range(5)]
W = [0] * 5
X = [0] * 5
W[1] = W[4] = 1/6
W[2] = W[3] = 1/3

# 修正FNF函数的符号错误：+4*X4而非-4*X4
def FNF(T, X1, X2, X3, X4):
    return -4*SIN(T) - X1 - X2 - X3 + 4*X4 + X4*(X4 - X3 + X2 - X1)

def FNQ(Z):
    return Z

print("t0,n,h,kmax")
T, N, H, KM = [float(_) for _ in input().split(" ")]
N_int = int(N)
for I in range(1, N_int):
    print(f"INPUT D{I}(X)")
    X[N_int - I] = float(input())

for K in range(1, int(KM)+1): 
    print(T, X[N_int], EXP(T)+SIN(T))
    for I in range(1, 5):
        for J in range(1, N_int+1):
            P = 0.5
            if I == 1:
                P = 0
            if I == 4:
                P = 1
            if J > 1:
                G[I][J] = FNQ(X[J - 1])
            X[J] += H * W[I] * G[I][J]
            A[J] = X[J] + P * H * G[I][J]
        # 计算四阶导数的FNF值
        G[I][1] = FNF(T + P*H, A[1], A[2], A[3], A[4])
    T += H

修正后，数值计算结果将与解析解一致，不再存在约1的偏差。

内容的提问来源于stack exchange，提问作者Ajaykrishnan R