DECLARE SUB RK4n (n%, x#, y#(), h#, y1#()) DECLARE SUB F (n%, x#, y#(), yp#()) '************************************************************************** '* Solve an ordinary system of first order differential equations (N<=10) * '* with initial conditions using a Runge-Kutta integration method. * '* ---------------------------------------------------------------------- * '* SAMPLE RUN: * '* Integrate first order DE system: * '* y1' = y2 * '* y2' = y3 * '* y3' = y4 * '* y4' = y5 * '* y5' = (45 * y3 * y4 * y5 - 40 * y4 * y4 * y4) / (9 * y3 * y3) * '* from x=0 to x=2 with initial conditions: y(0)..y(4) = 1 * '* (25 integration steps). * '* * '* X Y1 <------------------------> YN estimated * '* ---------------------------------------------------------------------- * '* 0.0000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 * '* 0.0800 1.0832871 1.0832862 1.0832443 1.0816193 1.0382909 * '* 0.1600 1.1735104 1.1734962 1.1731252 1.1657018 1.0606156 * '* 0.2400 1.2712454 1.2711671 1.2697821 1.2507235 1.0606138 * '* 0.3200 1.3771106 1.3768421 1.3732111 1.3346056 1.0307873 * '* 0.4000 1.4917679 1.4910565 1.4832168 1.4146287 0.9626611 * '* 0.4800 1.6159210 1.6143206 1.5993581 1.4873681 0.8471623 * '* 0.5600 1.7503129 1.7470979 1.7208913 1.5486703 0.6752862 * '* 0.6400 1.8957208 1.8897779 1.8467119 1.5936970 0.4391094 * '* 0.7200 2.0529492 2.0426461 1.9753055 1.6170637 0.1331626 * '* 0.8000 2.2228199 2.2058489 2.1047132 1.6131014 -0.2438996 * '* 0.8800 2.4061600 2.3793580 2.2325243 1.5762584 -0.6875083 * '* 0.9600 2.6037859 2.5629350 2.3559061 1.5016398 -1.1855881 * '* 1.0400 2.8164856 2.7561004 2.4716817 1.3856519 -1.7176006 * '* 1.1200 3.0449981 2.9581109 2.5764618 1.2266816 -2.2547153 * '* 1.2000 3.2899923 3.1679497 2.6668279 1.0257030 -2.7615081 * '* 1.2800 3.5520448 3.3843326 2.7395571 0.7866810 -3.1993222 * '* 1.3600 3.8316195 3.6057337 2.7918674 0.5166480 -3.5309960 * '* 1.4400 4.1290502 3.8304300 2.8216525 0.2253716 -3.7261918 * '* 1.5200 4.4445267 4.0565635 2.8276730 -0.0753912 -3.7662105 * '* 1.6000 4.7780875 4.2822163 2.8096752 -0.3729657 -3.6471428 * '* 1.6800 5.1296178 4.5051904 2.7684146 -0.6549877 -3.3805446 * '* 1.7600 5.4988544 4.7245867 2.7055842 -0.9105780 -2.9914571 * '* 1.8400 5.8853965 4.9378740 2.6236555 -1.1312684 -2.5142924 * '* 1.9200 6.2887215 5.1439429 2.5256617 -1.3115454 -1.9876112 * '* 2.0000 6.7082039 5.3416408 2.4149534 -1.4489720 -1.4489720 * '* * '* ---------------------------------------------------------------------- * '* 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_Rk4n DEFDBL A-H, O-Z DEFINT I-N OPTION BASE 0 ISIZE = 9 'maximum number of DEs DIM y(ISIZE), y1(ISIZE) DIM finesse AS INTEGER n = 5 ' size of DE system x0 = 0# ' starting x xl = 2# ' ending x kl = 25 ' number of tabulated values in x finesse = 20 h = (xl - x0) / kl / finesse' elementary integration step x = x0 FOR i = 0 TO n - 1 y(i) = 1# ' starting values NEXT i ' write header CLS PRINT PRINT " X Y1 <------------------------> YN estimated" PRINT " ----------------------------------------------------------------------" ' integration loop RK4n n, x, y(), h, y1() PRINT USING "####.####"; x; FOR i = 0 TO n - 1 PRINT USING "####.#######"; y(i); NEXT i PRINT FOR k = 1 TO kl * finesse x = x + h FOR i = 0 TO n - 1 y(i) = y1(i) NEXT i RK4n n, x, y(), h, y1() IF (k MOD finesse) = 0 THEN PRINT USING "####.####"; x; FOR i = 0 TO n - 1 PRINT USING "####.#######"; y(i); NEXT i PRINT END IF NEXT k PRINT " ----------------------------------------------------------------------" PRINT END ' DE system to integrate (n=5) SUB F (n, x, y(), yp()) yp(0) = y(1) yp(1) = y(2) yp(2) = y(3) yp(3) = y(4) yp(4) = (45# * y(2) * y(3) * y(4) - 40# * y(3) * y(3) * y(3)) / (9# * y(2) * y(2)) END SUB 'Integrate first order DE sytem from x to x+h using Runge-Kutta SUB RK4n (n, x, y(), h, y1()) DIM c1(n - 1), c2(n - 1), c3(n - 1), c4(n - 1), yy(n - 1) F n, x, y(), c1() h2 = h / 2# FOR i = 0 TO n - 1 yy(i) = y(i) + h2 * c1(i) NEXT i F n, x + h2, yy(), c2() FOR i = 0 TO n - 1 yy(i) = y(i) + h2 * c2(i) NEXT i F n, x + h2, yy(), c3() FOR i = 0 TO n - 1 yy(i) = y(i) + h * c3(i) NEXT i F n, x + h, yy(), c4() FOR i = 0 TO n - 1 y1(i) = y(i) + h * (c1(i) + 2# * c2(i) + 2# * c3(i) + c4(i)) / 6# NEXT i END SUB