'******************************************************************** '* Solve Y' = F(X,Y) with initial conditions using the Adams- * '* Moulton Prediction-Correction Method * '* ---------------------------------------------------------------- * '* SAMPLE RUN * '* (Integrate Y' = -Y + X/(1+X)^2 from X=0 to X=1 with initial * '* condition Y(0) = 1 ) * '* * '* X Y Y True |Y-Y True| * '* --------------------------------------------- * '* 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.869569 0.869565 0.000004 1 * '* 0.200000 0.833340 0.833333 0.000006 1 * '* 0.250000 0.800008 0.800000 0.000009 1 * '* 0.300000 0.769241 0.769231 0.000010 1 * '* 0.350000 0.740752 0.740741 0.000011 1 * '* 0.400000 0.714298 0.714286 0.000012 1 * '* 0.450000 0.689668 0.689655 0.000012 1 * '* 0.500000 0.666679 0.666667 0.000013 1 * '* 0.550000 0.645174 0.645161 0.000013 1 * '* 0.600000 0.625013 0.625000 0.000013 1 * '* 0.650000 0.606073 0.606061 0.000013 1 * '* 0.700000 0.588248 0.588235 0.000013 1 * '* 0.750000 0.571441 0.571429 0.000013 1 * '* 0.800000 0.555568 0.555556 0.000012 1 * '* 0.850000 0.540553 0.540541 0.000012 1 * '* 0.900000 0.526327 0.526316 0.000012 1 * '* 0.950000 0.512832 0.512821 0.000011 1 * '* 1.000000 0.500011 0.500000 0.000011 1 * '* ---------------------------------------------------------------- * '* REFERENCE: "Méthode de calcul numérique- Tome 2 - Programmes en * '* Basic et en Pascal By Claude Nowakowski, Edition du * '* P.S.I., 1984". * '* * '* Quick Basic Release By J-P Moreau, Paris. * '* (www.jpmoreau.fr) * '******************************************************************** 'See explanation file adambash.txt '--------------------------------- DEFINT I-N DEFDBL A-H, O-Z OPTION BASE 0 'index begins from zero DIM X(3), Y(3) H = .05 'integration step X(0) = 0#: Y(0) = 1# 'Initial conditions EC = .000001# 'Precision F\$ = " ##.###### ##.###### ##.###### ##.###### ##" CLS PRINT " X Y Y True |Y-Y True| " PRINT " ---------------------------------------------" PRINT USING F\$; X(0); Y(0); Y(0); X(0) 'Start with Runge-Kutta FOR K = 0 TO 1 XX = X(K): YY = Y(K): GOSUB 1000: C1 = C XX = X(K) + H / 2#: YY = Y(K) + H / 2# * C1: GOSUB 1000: C2 = C YY = Y(K) + H / 2# * C2: GOSUB 1000: C3 = C X(K + 1) = X(K) + H XX = X(K + 1): YY = Y(K) + H * C3: GOSUB 1000: C4 = C Y(K + 1) = Y(K) + H * (C1 + 2 * C2 + 2 * C3 + C4) / 6# XX = X(K + 1): GOSUB 2000: ER = ABS(F - Y(K + 1)) PRINT USING F\$; X(K + 1); Y(K + 1); F; ER NEXT K 100 K = 2 XX = X(K): YY = Y(K): GOSUB 1000: C1 = C XX = X(K - 1): YY = Y(K - 1): GOSUB 1000: C2 = C XX = X(K - 2): YY = Y(K - 2): GOSUB 1000: C3 = C X(K + 1) = X(K) + H YP = Y(K) + H / 12# * (23 * C1 - 16 * C2 + 5 * C3) L = 0 200 XX = X(K + 1): YY = YP: GOSUB 1000: C1 = C XX = X(K): YY = Y(K): GOSUB 1000: C2 = C XX = X(K - 1): YY = Y(K - 1): GOSUB 1000: C3 = C YC = Y(K) + H / 12# * (5 * C1 + 8 * C2 - C3) 'PRINT YC IF ABS(YP - YC) > EC THEN YP = YC: L = L + 1: GOTO 200 END IF Y(K + 1) = YC: XX = X(K + 1): GOSUB 2000: ER = ABS(F - Y(K + 1)) PRINT USING F\$; X(K + 1); Y(K + 1); F; ER; L FOR K = 0 TO 2 X(K) = X(K + 1): Y(K) = Y(K + 1) NEXT K IF X(2) < 1# THEN GOTO 100 END IF END 1000 C = -YY + XX / ((1# + XX) ^ 2) RETURN 2000 F = 1# / (1# + XX) RETURN