/******************************************************************* * 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.800009 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". * * * * C++ Release By J-P Moreau, Paris. * * (www.jpmoreau.fr) * ******************************************************************** See explanation file adambash.txt --------------------------------- */ #include #include //Labels: e100, e200 double X[4], Y[4]; double C1,C2,C3,C4, EC, ER, F, H, XX, YY, YC, YP; int K, L; double Sqr(double x) { return x*x; } void P1000(double XX, double YY, double *C) { *C = -YY + XX / (Sqr(1.0 + XX)); } void P2000(double XX, double *F) { *F = 1.0 / (1.0 + XX); } void main() { H = 0.05; // integration step X[0] = 0.0; Y[0] = 1.0; // Initial conditions EC = 0.000001; // Precision printf("\n X Y Y True |Y-Y True|\n"); printf(" ---------------------------------------------\n"); printf(" %9.6f %9.6f %9.6f %9.6f\n", X[0], Y[0], Y[0], X[0]); // Start with Runge-Kutta for (K=0; K<2; K++) { XX = X[K]; YY = Y[K]; P1000(XX,YY,&C1); XX = X[K] + H / 2.0; YY = Y[K] + H / 2.0 * C1; P1000(XX,YY,&C2); YY = Y[K] + H / 2.0 * C2; P1000(XX,YY,&C3); X[K + 1] = X[K] + H; XX = X[K + 1]; YY = Y[K] + H * C3; P1000(XX,YY,&C4); Y[K + 1] = Y[K] + H * (C1 + 2 * C2 + 2 * C3 + C4) / 6.0; XX = X[K + 1]; P2000(XX,&F); ER = fabs(F - Y[K + 1]); printf(" %9.6f %9.6f %9.6f %9.6f\n", X[K+1], Y[K+1], F, ER); } e100: K = 2; XX = X[K]; YY = Y[K]; P1000(XX,YY,&C1); XX = X[K - 1]; YY = Y[K - 1]; P1000(XX,YY,&C2); XX = X[K - 2]; YY = Y[K - 2]; P1000(XX,YY,&C3); X[K + 1] = X[K] + H; YP = Y[K] + H / 12.0 * (23.0 * C1 - 16.0 * C2 + 5.0 * C3); L = 0; e200: XX = X[K + 1]; YY = YP; P1000(XX,YY,&C1); XX = X[K]; YY = Y[K]; P1000(XX,YY,&C2); XX = X[K - 1]; YY = Y[K - 1]; P1000(XX,YY,&C3); YC = Y[K] + H / 12.0 * (5.0 * C1 + 8.0 * C2 - C3); if (fabs(YP - YC) > EC) { YP = YC; L++; goto e200; } Y[K + 1] = YC; XX = X[K + 1]; P2000(XX,&F); ER = fabs(F - Y[K + 1]); printf(" %9.6f %9.6f %9.6f %9.6f %d\n", X[K+1], Y[K+1], F, ER, L); for (K=0; K<3; K++) { X[K] = X[K + 1]; Y[K] = Y[K + 1]; } if (X[2] < 1.0) goto e100; } // end of file Adammoul.cpp