DECLARE SUB RK4 (l%, x#, y#, z#, h#, x1#, y1#, z1#) DECLARE FUNCTION F# (x#, y#, z#) DECLARE FUNCTION G# (l%, x#, y#, z#) DECLARE FUNCTION FX# (x#) '************************************************************************** '* Solve a boundary value problem for a second order DE of the form: * '* y" - p(x)*y = f(x) with conditions: y(x0)=a and y(x1)=b (x0<=x<=x1) * '* using a Runge-Kutta integration method. * '* ---------------------------------------------------------------------- * '* Explanations: if y0(x) is a particular solution of y"-p(x)y = f(x) (1) * '* and y1(x), y2(x) two lineary independant solutions of y"-p(x)y = 0 (2) * '* (right hand side removed), then the general solution of (1) has the * '* form: y(x) = y0(x) + c1*y1(x) + c2*y2(x), the constants c1, c2 are * '* obtained from boundary conditions. * '* A more efficient way to proceed is: 1) calculate a particular solution * '* y1(x) of (2) such as y1(x0)=0 and 2) calculate a particular solution * '* y2(x) of (1) such as y2(x0)=a, then the general solution of (1) has the* '* form: y(x) = y1(x) + C * y2(x), constant C is determined by: * '* C * y2(x1) + y1(x1) = b * '* ---------------------------------------------------------------------- * '* SAMPLE RUN: * '* Integrate y" - (1/x)y' - (3/xý)y = 5xý from x0=1 to x1=2 with boundary * '* conditions: y(x0)=11 and y(x1)=36. The exact solution is: * '* y(x) = 8/x + 2xý + x4 * '* The initial 2nd order problem is replaced by the following first order * '* DE systems: * '* 1st system 2nd system * '* y' = z y' = z * '* z' = z/xý + 3y/xý z' = z/xý + 3y/xý + 5xý * '* with z'(x0)=0 with z'(x0)=11 * '* * '* X Y estimated Y true Absolute Error * '* ---------------------------------------------------- * '* 1.0000 11.00000000 11.00000000 0.00000000 * '* 1.0500 11.14980444 11.14980387 0.00000057 * '* 1.1000 11.39882827 11.39882727 0.00000100 * '* 1.1500 11.74727930 11.74727799 0.00000131 * '* 1.2000 12.19626820 12.19626667 0.00000153 * '* 1.2500 12.74765794 12.74765625 0.00000169 * '* 1.3000 13.40394794 13.40394615 0.00000178 * '* 1.3500 14.16818401 14.16818218 0.00000184 * '* 1.4000 15.04388756 15.04388571 0.00000185 * '* 1.4500 16.03499946 16.03499763 0.00000183 * '* 1.5000 17.14583511 17.14583333 0.00000178 * '* 1.5500 18.38104827 18.38104657 0.00000170 * '* 1.6000 19.74560159 19.74560000 0.00000159 * '* 1.6500 21.24474257 21.24474110 0.00000147 * '* 1.7000 22.88398367 22.88398235 0.00000132 * '* 1.7500 24.66908597 24.66908482 0.00000115 * '* 1.8000 26.60604540 26.60604444 0.00000096 * '* 1.8500 28.70108132 28.70108057 0.00000075 * '* 1.9000 30.96062683 30.96062632 0.00000052 * '* 1.9500 33.39132062 33.39132035 0.00000027 * '* 2.0000 36.00000000 36.00000000 0.00000000 * '* * '* ---------------------------------------------------------------------- * '* REFERENCE: "Methode de calcul numerique- Tome 2 - Programmes en Basic * '* et en Pascal By Claude Nowakowski, Edition du P.S.I., 1984"* '* [BIBLI 04]. * '* * '* Quick Basic Release By J-P Moreau, Paris. * '* (www.jpmoreau.fr) * '************************************************************************** 'Program Limits DEFDBL A-H, O-Z DEFINT I-N OPTION BASE 0 ISIZE = 100 DIM x(ISIZE), y(3, ISIZE), z(3, ISIZE) x0 = 1# 'starting x xl = 2# 'ending x kl = 20 'number of steps in x h = (xl - x0) / kl 'integration step ya = 11# 'starting value for y yb = 36# 'final value for y'' x(0) = x0 'integration loop FOR l = 1 TO 2 y(l, 0) = (l - 1) * ya: z(l, 0) = 1# FOR k = 0 TO kl - 1 RK4 l, x(k), y(l, k), z(l, k), h, x(k + 1), y(l, k + 1), z(l, k + 1) NEXT k NEXT l 'solution is y0=y2+cc*y1, cc is such as y(x1)=yb cc = (yb - y(2, kl)) / y(1, kl) 'write header CLS PRINT PRINT " X Y estimated Y true Absolute Error " PRINT " ----------------------------------------------------" 'write kl+1 result lines FF\$ = " ###.########" FOR k = 0 TO kl y(0, k) = y(2, k) + cc * y(1, k) yex = FX(x(k)) 'exact value er = ABS(yex - y(0, k))'current error PRINT USING " #.##"; x(k); PRINT USING FF\$; y(0, k); PRINT USING FF\$; yex; PRINT USING FF\$; er NEXT k PRINT " ----------------------------------------------------" END 'of main program 'y'=z FUNCTION F (x, y, z) F = z END FUNCTION 'exact solution: y=8/x + 2x3 + x4 FUNCTION FX (x) FX = 8# / x + x * x * x * (x + 2#) END FUNCTION 'l=1: z'=z/x²+3y/x² or l=2: z'=z/x²+3y/x²+5x² FUNCTION G (l, x, y, z) gg = z / x + 3! * y / (x * x) IF l = 2 THEN gg = gg + 5# * x * x G = gg END FUNCTION 'Integrate sytem from x to x+h using Runge-Kutta SUB RK4 (l, x, y, z, h, x1, y1, z1) c1 = F(x, y, z) d1 = G(l, x, y, z) h2 = h / 2# c2 = F(x + h2, y + h2 * c1, z + h2 * d1) d2 = G(l, x + h2, y + h2 * c1, z + h2 * d1) c3 = F(x + h2, y + h2 * c2, z + h2 * d2) d3 = G(l, x + h2, y + h2 * c2, z + h2 * d2) c4 = F(x + h, y + h * c3, z + h * d3) d4 = G(l, 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