(71B) Implicit Trapezoidal Integration Method for ODE
02-23-2024, 10:06 PM (This post was last modified: 02-23-2024 10:07 PM by Namir.)
Post: #1
 Namir Senior Member Posts: 1,085 Joined: Dec 2013
(71B) Implicit Trapezoidal Integration Method for ODE
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
02-25-2024, 05:02 PM
Post: #2
 Valentin Albillo Senior Member Posts: 1,100 Joined: Feb 2015
RE: (71B) Implicit Trapezoidal Integration Method for ODE
.
Your 1,000th post, Namir, congratulations !

V.

All My Articles & other Materials here:  Valentin Albillo's HP Collection

02-25-2024, 10:57 PM
Post: #3
 Namir Senior Member Posts: 1,085 Joined: Dec 2013
RE: (71B) Implicit Trapezoidal Integration Method for ODE
(02-25-2024 05:02 PM)Valentin Albillo Wrote:  .
Your 1,000th post, Namir, congratulations !

V.

Thank you .. thank you ... <bowns down>

:-)

Namir
 « Next Oldest | Next Newest »

User(s) browsing this thread: 1 Guest(s)