I am trying to convert code in BASIC to python that numerically computes the solution of the 4th order IVP (using Runge-Kutta 4th Order) defined as :
and the python code I wrote looks like this:
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
However the actual analytical results don't match with the numerical ones. They are almost always off by around 1.
Here is the BASIC code which supposedly works:
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= .S: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 1:T = T + H:NEXT K: END
