(71B) Implicit Trapezoidal Integration Method for ODE

[Namir]

Implicit Trapezoidal Integration Method for ODE
================================================

Given f(x,y) and (x0,y0) as the initial point, xf as the final point, h is the integration step, toler as the tolerance value.

The function to solve for y1 is:

F(y1) = y1 - y0 - h/2*(f(x0,y0) + f(x1,y1))

The loop to calculate y at x=xf is:

Code:

While x0 < xf
  x1 = x0 + h
  y1 = y0 + h*f(x0,y0)
  diff = 2*toler
  i = 1
  max = 100
  While abs(diff) >= toler and i <= max
    i += 1
    hx = 0.01*(1 + abs(y1)
    F0 = F(y1)
    Fp = F(y1+hx)
    Fm - F(y1-hx)
    diff = 2*hx*F0/(Fp - Fm)
    x -= diff
  end
  y0 = y1
  x0 = x1
end

Instructions
============

1) Define the differential equation at line 30. FNF should be a function of X and/or Y.
2) Edit the DATA statement in line 60 to contain the following values:
2.1) The value for the initial X.
2.2) The value for the initial Y.
2.3) The value for the final X.
2.4) The value for the integration step.
2.5) The tolerance value for solving the implicit step.
2.6) The number of results to skip storing in the results matrix R. Can be 0.
3) Press RUN. The program displays the current values of X. The program then displays the final value of Y. Finally the program displays the paired values of X and Y stored in the results matrix R.


Example
=======

given the ODE dy/dx = -2*y, x0=0, y0=1, xf=1, h=0.1, tolerance=1e-6. The matrix R stores the rows of (x, y) values, skipping every 0 values.

Press RUN. The program displays X values and then the final Y value as "Y = 0.13443063275". Press [f][Cont] to view the paired values of X and Y stored in matrix R. The last value is the final value of Y.


Program Listing
===============

Code:

10 REM  IMPLICIT TRAPEZOIDAL ODE
20 DESTROY ALL
30 DEF FNF(X,Y) = -2*Y
40 DEF FNG(X0,Y0,X1,Y1,H) = Y1 - Y0 - H/2*(F(X0,Y0) + F(X1,Y1))
50 READ X0,Y0,X9,H,T0,S
60 DATA 0,1,1,0.1,1E-6
70 S0 = 1 +INT((X9-X0)/H)
80 DIM R(S0)
90 C = S
100 R0 = 1
110 R(1,1) = X0 @ R(1,2) = Y0
120 'WHILE1':
130 IF X0 >=  X9 THEN GOTO 'END1'
140 X1 = X0 + H
150 REM INITILIAL GUESS FOR Y1
160 Y1 = Y0 + H*FX(X0, Y0) 
170 D = 2*T0
180 I = 1
190 M9 = 100
200 'WHILE2':
210 IF ABS(DIFF)< T0 OR I > M9 THEN GOTO 'END2'
220 I = I + 1
230 H0 = 0.01*(1+ABS(Y1))
240 F0 = G(X0,Y0,X1,Y1,H)
250 F1 = G(X0,Y0,X1,Y1+H0,H)
260 F2 = G(X0,Y0,X1,Y1-H0,H)
270 D = 2*H0*F0/(F1 - F2)
280 Y1 = Y1 - D
290 GOTO 'WHILE2'
300 'END2':
310 Y0 = Y1
320 X0 = X1
330 IF C > 0 THEN C = C - 1 @ GOTO 'END3'
340 R0 = R0 + 1
350 R(R0,1) = X0
360 R(R0,2) = Y0
370 C = S
380 'END3':
390 GOTO 'WHILE1'
400 'END1':
410 IF R(R0,1) = X9 THEN GOTO 'END4'
420 R0 = R0 + 1
430 R(R0,1) = X0
440 R(R0,2) = Y0
450 'END4':
460 DISP 'Y = ";Y0 @ PAUSE
470 REM YOU CAN REPLACE THE WAIT WITH PAUSE IN THE NEXT LOOP
480 FOR I=1 TO S0
490 DISP "X = ";R(I,1) @ WAIT 2
500 DISP "Y = ";R(I,2) @ WAIT 2
510 NEXT I
520 END