!************************************************************************ !* Test program for the subroutine awp from the module AWP to solve a * !* first order ordinary system of differential equations using a Runge-* !* Kutta method with step size and precision controls * !* * !* Scope of program: * !* ================= * !* The program reads the input data from the file stdin and writes * !* output onto the file stdout. * !* * !* After reading the input, it is put out for control purposes; then * !* the adaptive initial value solver is executed and the computed * !* results are put out. * !* * !* To solve a differential equation, please proceed as in examples in * !* file t_dgls.f90. * !* * !* Construction of input data files: * !* ================================= * !* bspnummer Number of DE system in t_dgls.cpp * !* epsabs desired absolute error bound * !* epsrel desired relative error bound * !* x0 initial x-value * !* y0[0] \ initial y-value at x0 * !* ... > * !* y0[n-1] / * !* h starting step size * !* fmax maximal number of calls of right hand side * !* xend final desired x-value * !* * !* The number n of differential equations follows from the number of * !* the DE system chosen; it is stored in conjunction with the right * !* hand side function in t_dgls.cpp. * !* -------------------------------------------------------------------- * !* SAMPLE RUNS: * !* * !* Example: 0 * !* absolute error bound epsabs: 1d-7 * !* relative error bound:epsrel: 1d-9 * !* initial x-value x0: 0 * !* fonction value y0(0) at x0: 0 * !* fonction value y0(1) at x0: 0 * !* starting step size h: 0.05 * !* desired final x-value: 2 * !* maximal number of calls for r.h. side: 500 * !* * !* System of DEs: * !* -------------- * !* y1' = y1 * y2 + cos(x) - 0.5 * sin(2.0*x) * !* y2' = y1 * y1 + y2 * y2 - (1 + sin(x)) * !* * !* x = 2.000000000000000 * !* approximate solution y1(x) = 0.91931625E-01 * !* approximate solution y2(x) = -0.13638550E+01 * !* * !* Number of function calls: 415 * !* Final step size ........: 3.056893386511712E-002 * !* Error code .............: 0 * !* * !* Example: 1 * !* absolute error bound epsabs: 1e-8 * !* relative error bound:epsrel: 1e-10 * !* initial x-value x0: 0 * !* fonction value y0(0) at x0: 1 * !* starting step size h: 0.01 * !* desired final x-value: 1 * !* maximal number of calls for r.h. side: 500 * !* * !* System of DEs: * !* y' = -y + x/((1+x)*(1+x)) * !* * !* x = 1.000000000000000 * !* approximate solution y1(x) = 0.500000000e+00 * !* * !* Number of function calls: 175 * !* Final step size ........: 4.554846520992836E-002 * !* Error code .............: 0 * !* * !* Example: 3 * !* absolute error bound epsabs: 1e-6 * !* relative error bound:epsrel: 1e-8 * !* initial x-value x0: 0 * !* fonction value y0(0) at x0: 1 * !* fonction value y0(1) at x0: 1 * !* fonction value y0(2) at x0: 1 * !* fonction value y0(3) at x0: 1 * !* fonction value y0(4) at x0: 1 * !* starting step size h: 0.05 * !* desired final x-value: 2 * !* maximal number of calls for r.h. side: 2000 * !* * !* System of DEs: * !* -------------- * !* y1' = y2 * !* y2' = y3 * !* y3' = y4 * !* y4' = y5 * !* y5' = (45 * y3 * y4 * y5 - 40 * y4 * y4 * y4) / (9 * y3 * y3) * !* * !* x = 2.000000000000000 * !* approximate solution y1(x) = 0.67082039E+01 * !* approximate solution y2(x) = 0.53416408E+01 * !* approximate solution y3(x) = 0.24149534E+01 * !* approximate solution y4(x) = -0.14489721E+01 * !* approximate solution y5(x) = -0.14489720E+01 * !* * !* Number of function calls: 505 * !* Final step size ........: 3.157943229233413E-002 * !* Error code .............: 0 * !* * !* -------------------------------------------------------------------- * !* REF.: "Numerical Algorithms with C, By Gisela Engeln-Muellges * !* and Frank Uhlig, Springer-Verlag, 1996" [BIBLI 11]. * !* * !* F90 Release 1.1 By J-P Moreau, Paris. * !* (www.jpmoreau.fr) * !* * !* Release 1.1 (12/19/2006): added two more examples. * !************************************************************************ ! To link with files awp.f90 and t_dgls.f90. ! ----------------------------------------- Program Test_AWP real*8 xbegin,xend,epsabs,epsrel,h integer bspn,method,n,fmax,ncalls,error real*8 Y(0:9) character*70 dgltxt(0:4) !auxiliary vectors real*8 yhilf(0:9), k1(0:9), k2(0:9), k3(0:9), k4(0:9), k5(0:9), k6(0:9) write(*,10,advance='no'); read *, bspn !# example in t_dgls call settxt(bspn,n,dgltxt) !initialize n and dgltxt write(*,30,advance='no'); read *, epsabs write(*,40,advance='no'); read *, epsrel write(*,50,advance='no'); read *, xbegin do i=0, n-1 write(*,60,advance='no') i; read *, Y(i) end do write(*,70,advance='no'); read *, h write(*,80,advance='no'); read *, xend write(*,90,advance='no'); read *, fmax method=6 !England embedding formulas !call integration procedure with step size control call awp(xbegin,xend,bspn,n,Y,epsabs,epsrel,h,method,fmax,ncalls,error, & yhilf, k1, k2, k3, k4, k5, k6) write(*,*) ' ' write(*,*) ' System of DEs:' do i=0, n-1 write(*,*) dgltxt(i) end do if (error.ne.0) then print *,' ' print *,' Error in awp:', error else !write results print *,' ' print *,' x= ', xend do i=0, n-1 write(*,100) i+1, Y(i) end do print *,' ' print *,' Number of function calls: ', ncalls print *,' Final step size ........: ', h print *,' Error code .............: ', error print *,' ' end if 10 format(/' Example: ') 30 format(' absolute error bound epsabs: ') 40 format(' relative error bound epsrel: ') 50 format(' initial x-value x0: ') 60 format(' fonction value y0(',I1,') at x0: ') 70 format(' starting step size h: ') 80 format(' desired final x-value: ') 90 format(' maximal number of calls for r.h. side: ') 100 format(' approximate solution y',I1,'(x) = ',E15.8) END !end of file tawp.f90