DECLARE SUB RK4 (X#, Y#, Z#, h#, x1#, y1#, z1#) DECLARE FUNCTION F# (X#, Y#, Z#) DECLARE FUNCTION G# (X#, Y#, Z#) '************************************************************************** '* Solve a first order DE system (N=2) of the form: * '* y' = F(x,y,z), z'=G(x,y,z) with conditions: y(x0)=a and z(x0)=b from * '* x=x0 to x=x1 using a Runge-Kutta integration method. * '* ---------------------------------------------------------------------- * '* SAMPLE RUN: * '* Integrate first order DE systemy: * '* y' = y*z + cos(x) - 0.5*sin(2*x) * '* z' = y*y + z*z - (1 + sin(x)) * '* from x=0 to x=2 with initial conditions: y(0)=0 and z(0)=0 * '* (25 integration steps). * '* * '* X Y estimated Z estimated * '* ------------------------------------ * '* 0.0000 0.0000000 0.0000000 * '* 0.0800 0.0765518 -0.0828580 * '* 0.1600 0.1452904 -0.1700792 * '* 0.2400 0.2050556 -0.2597237 * '* 0.3200 0.2550329 -0.3500533 * '* 0.4000 0.2948004 -0.4396008 * '* 0.4800 0.3243379 -0.5271915 * '* 0.5600 0.3440043 -0.6119278 * '* 0.6400 0.3544874 -0.6931461 * '* 0.7200 0.3567373 -0.7703644 * '* 0.8000 0.3518908 -0.8432315 * '* 0.8800 0.3411943 -0.9114874 * '* 0.9600 0.3259299 -0.9749377 * '* 1.0400 0.3073506 -1.0334409 * '* 1.1200 0.2866242 -1.0869048 * '* 1.2000 0.2647908 -1.1352886 * '* 1.2800 0.2427308 -1.1786041 * '* 1.3600 0.2211465 -1.2169165 * '* 1.4400 0.2005537 -1.2503410 * '* 1.5200 0.1812830 -1.2790362 * '* 1.6000 0.1634888 -1.3031952 * '* 1.6800 0.1471644 -1.3230351 * '* 1.7600 0.1321612 -1.3387869 * '* 1.8400 0.1182104 -1.3506867 * '* 1.9200 0.1049465 -1.3589674 * '* 2.0000 0.0919314 -1.3638530 * '* * '* ---------------------------------------------------------------------- * '* REFERENCE: "Methode de calcul numerique- Tome 2 - Programmes en Basic * '* et en Pascal By Claude Nowakowski, Edition du P.S.I., 1984"* '* [BIBLI 04]. * '* * '* Basic Release By J-P Moreau, Paris. * '* (www.jpmoreau.fr) * '************************************************************************** 'Program Test_Rk4 DEFDBL A-H, O-Z DEFINT I-N OPTION BASE 0 ISIZE = 100 DIM X(ISIZE), Y(ISIZE), Z(ISIZE) CLS x0 = 0# 'starting x xl = 2# 'ending x kl = 25 'number of steps in x h = (xl - x0) / kl 'integration step X(0) = x0 Y(0) = 0#: Z(0) = 0# 'starting values 'integration loop FOR k = 0 TO kl - 1 RK4 X(k), Y(k), Z(k), h, X(k + 1), Y(k + 1), Z(k + 1) NEXT k 'write header PRINT PRINT " X Y estimated Z estimated " PRINT " -------------------------------------" 'write kl+1 result lines FF\$ = " ###.#######" FOR k = 0 TO kl PRINT USING " #.##"; X(k); PRINT USING FF\$; Y(k); PRINT USING FF\$; Z(k) NEXT k PRINT " -------------------------------------" PRINT END 'of main program 'y'=yz + cos(x) - 0.5sin(2x) FUNCTION F (X, Y, Z) F = Y * Z + COS(X) - .5# * SIN(2# * X) END FUNCTION 'z'=yy + zz -(1+sin(x)) FUNCTION G (X, Y, Z) G = Y * Y + Z * Z - (1# + SIN(X)) END FUNCTION 'Integrate sytem from x to x+h using Runge-Kutta SUB RK4 (X, Y, Z, h, x1, y1, z1) c1 = F(X, Y, Z) d1 = G(X, Y, Z) h2 = h / 2# c2 = F(X + h2, Y + h2 * c1, Z + h2 * d1) d2 = G(X + h2, Y + h2 * c1, Z + h2 * d1) c3 = F(X + h2, Y + h2 * c2, Z + h2 * d2) d3 = G(X + h2, Y + h2 * c2, Z + h2 * d2) c4 = F(X + h, Y + h * c3, Z + h * d3) d4 = G(X + h, Y + h * c3, Z + h * d3) x1 = X + h y1 = Y + h * (c1 + 2# * c2 + 2# * c3 + c4) / 6# z1 = Z + h * (d1 + 2# * d2 + 2# * d3 + d4) / 6# END SUB