0

This is the BASIC code for solving a nth order IVP using the RK 4 method:

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(I,1)=FNF(T+P*H,A(1),A(2),A(3),A(4))
NEXT I:T = T + H:NEXT K: END

The problem is how is X(1) or A(1) calculated i.e. what happens when J = 1 ? This part :

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)

only gets run when J > 1 and this part:

G(I,1)=FNF(T+P*H,A(1),A(2),A(3),A(4))

requires A(1) to be calculated.

4
  • "This part :.... (code)....only gets run when J > 1". I think/know that the IF statement ends at the end of the line, so the X(J)=... and A(J)=... lines are executed for every value of J.
    – Luuk
    Commented Feb 11 at 18:45
  • @Luuk yeah i thought about that but the G(I, K) is not initialized right ? the initialization part i.e. THEN G(I,J) = FNQ(X(J-1)) is run only when J > 1 Commented Feb 11 at 18:58
  • Does the value not default to 0 ? (It's been a while since I used BASIC, and this might even depend on the implementation of it..... )
    – Luuk
    Commented Feb 11 at 19:17
  • @Luuk Oh actually i am new to BASIC but that probably is the case Commented Feb 11 at 19:49

0

Browse other questions tagged or ask your own question.