!******************************************************************** !* 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 * !* 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 * !* ----------------------------------------------- * !* * !* F90 Release By J-P Moreau, Paris. * !* (www.jpmoreau.fr) * !******************************************************************** Program AdamsBashforth implicit real*8 a-h,o-z dimension B(0:3),X(0:3),Y(0:3) H = 0.05d0 !integration step !Initial conditions X(0)=0.d0 !starting abscissa X0 X1=1.d0 !ending abscissa X1 Y(0)=1.d0 !Initial value of Y at X0 !write header print *,' X Y Y exact Error ' print *,' -----------------------------------------------' !write initial line write(*,10) X(0), Y(0), Y(0),' 0' !use Runge-Kutta to start do K=0, 1 call RK4(K,H,X,Y); YEX=FX(X(K+1)); ER=DABS(YEX-Y(K+1)) write(*,11) X(K+1), Y(K+1), YEX, ER end do !main integration loop do while(X(2)