'******************************************************************** '* Solve Y' = F(X,Y) with Initial Condition Y(X0)=Y0 using the * '* Adams-Bashforth Method * '* ---------------------------------------------------------------- * '* REFERENCE: "Méthode de calcul numérique- Tome 2 - Programmes en * '* Basic et en Pascal By Claude Nowakowski, Edition du * '* P.S.I., 1984" [4]. * '* ---------------------------------------------------------------- * '* SAMPLE RUN: * '* (Solve Y' = -Y + X/((1+X)*(1+X)) with Y(0)=1 for X0=0 up to * '* X1=1.0, exact solution is Y = 1 / (1+X). * '* * '* X Y Y exact Error * '* ----------------------------------------------- * '* 0.000000 1.000000 1.000000 0 * '* 0.000000 1.000000 1.000000 0.000000 * '* 0.050000 0.952381 0.952381 0.000000 * '* 0.100000 0.909091 0.909091 0.000000 * '* 0.150000 0.869525 0.869565 0.000040 * '* 0.200000 0.833265 0.833333 0.000068 * '* 0.250000 0.799910 0.800000 0.000090 * '* 0.300000 0.769125 0.769231 0.000106 * '* 0.350000 0.740623 0.740741 0.000117 * '* 0.400000 0.714160 0.714286 0.000125 * '* 0.450000 0.689525 0.689655 0.000131 * '* 0.500000 0.666533 0.666667 0.000134 * '* 0.550000 0.645026 0.645161 0.000135 * '* 0.600000 0.624865 0.625000 0.000135 * '* 0.650000 0.605926 0.606061 0.000134 * '* 0.700000 0.588103 0.588235 0.000133 * '* 0.750000 0.571298 0.571429 0.000131 * '* 0.800000 0.555428 0.555556 0.000128 * '* 0.850000 0.540416 0.540541 0.000125 * '* 0.900000 0.526194 0.526316 0.000121 * '* 0.950000 0.512703 0.512821 0.000118 * '* 1.000000 0.499886 0.500000 0.000114 * '* ----------------------------------------------- * '* * '* Basic Release By J-P Moreau, Paris. * '* (www.jpmoreau.fr) * '******************************************************************** DEFDBL A-H, O-Z DEFINT I-N OPTION BASE 0 H = .05 'integration step DIM B(4), X(4), Y(4) 'Initial conditions X(0) = 0# 'starting X X1 = 1# 'ending X Y(0) = 1# 'initial Y 'write header CLS PRINT PRINT " X Y Y exact Error " PRINT " ------------------------------------------------" 'write initial line F\$ = "#####.######" PRINT USING F\$; X(0); PRINT USING F\$; Y(0); PRINT USING F\$; Y(0); PRINT " 0" 'use Runge-Kutta to start FOR K = 0 TO 1 GOSUB 1000 'call RK4 XX = X(K + 1): GOSUB 600: ER = ABS(FX - Y(K + 1)) PRINT USING F\$; X(K + 1); PRINT USING F\$; Y(K + 1); PRINT USING F\$; FX; PRINT USING F\$; ER NEXT K 'main integration loop WHILE X(2) < X1 FOR I = 1 TO 3 XX = X(3 - I): YY = Y(3 - I): GOSUB 500 B(I) = F NEXT I X(3) = X(2) + H Y(3) = Y(2) + H * (23# * B(1) - 16# * B(2) + 5# * B(3)) / 12# XX = X(3): GOSUB 600: ER = ABS(Y(3) - FX) PRINT USING F\$; X(3); PRINT USING F\$; Y(3); PRINT USING F\$; FX; PRINT USING F\$; ER FOR K = 0 TO 2 X(K) = X(K + 1): Y(K) = Y(K + 1) NEXT K WEND PRINT " ------------------------------------------------" END 'User defined function Y'=F(XX,YY) 500 F = -YY + XX / ((1# + XX) * (1# + XX)) RETURN 'Exact solution Y=FX(XX) 600 FX = 1# / (1# + XX) RETURN 'Runge-Kutta method to calculate first points only 1000 XX = X(K): YY = Y(K): GOSUB 500: C1 = F XX = X(K) + H / 2#: YY = Y(K) + H / 2# * C1: GOSUB 500: C2 = F XX = X(K) + H / 2#: YY = Y(K) + H / 2# * C2: GOSUB 500: C3 = F XX = X(K) + H: YY = Y(K) + H * C3: GOSUB 500: C4 = F X(K + 1) = X(K) + H Y(K + 1) = Y(K) + H * (C1 + 2 * C2 + 2 * C3 + C4) / 6# RETURN 'end of file adambash.bas