/******************************************************************* * 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" [BIBLI 04]. * * ---------------------------------------------------------------- * * 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.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 * * ----------------------------------------------- * * * * C++ Release By J-P Moreau, Paris. * * (www.jpmoreau.fr) * *******************************************************************/ #include #include const double H=0.05; // integration step double B[4],X[4],Y[4]; int I,K; double C1,C2,C3,C4,ER,X1,YEX; // User defined function Y'=F(X,Y) double F(double X, double Y) { return -Y + X/((1.0+X)*(1.0+X)); } // Exact solution Y=FX(X) double FX(double X) { return 1.0/(1.0+X); } // Runge-Kutta method to calculate first points only void RK4() { C1=F(X[K],Y[K]); C2=F(X[K]+H/2,Y[K]+H/2*C1); C3=F(X[K]+H/2,Y[K]+H/2*C2); C4=F(X[K]+H,Y[K]+H*C3); X[K+1]=X[K]+H; Y[K+1]=Y[K]+H*(C1+2*C2+2*C3+C4)/6; } void main() { // Initial conditions X[0]=0.0; //starting X X1=1.0; //ending X Y[0]=1.0; //initial Y // write header printf(" X Y Y exact Error \n"); printf(" -----------------------------------------------\n"); // write initial line printf("%12.6f%12.6f%12.6f 0.000000\n", X[0],Y[0],Y[0]); // use Runge-Kutta to start for (K=0; K<2; K++) { RK4(); YEX=FX(X[K+1]); ER=fabs(YEX-Y[K+1]); printf("%12.6f%12.6f%12.6f %f\n", X[K+1],Y[K+1],YEX,ER); } // main integration loop while (X[2]