c Example A c c This example shows how to use RKC. It solves a system of ODEs that c arise from semi-discretization of the reaction-diffusion equation c c U = U + (1 - U)*U**2 for t >= 0, 0 <= x <= 10 c t xx c c Dirichlet boundary conditions specify U(0,t) and U(10,t) for all t >= 0 c and the initial values U(x,0) are specified. These values are taken from c an analytical solution that is evaluated in sol(x,t) so that the numerical c solution can be compared to a known solution. c c A semi-discretization of the PDE is obtained by choosing a set of points c {x_i} in [0, 10] and approximating U(x_i,t) by a function y_i(t). Here c neqn+2 equally spaced points x_i are used for neqn = 99. When the second c partial derivative of U with respect to x is approximated by central c differences, a system of neqn ODEs is obtained for the y_i(t). The c initial values y_i(0) are given by U(x_i,0) = sol(x_i,0). c c A common way to present the computed results is to plot approximations c to U(x,tout) on [0, 10] for a selection of times tout. This example c shows how to compute approximations to the y_i(t) at these specific times. c They are written to an output file for plotting; the most convenient way c to do this will depend on the system and the plotting package used. c c Because an analytical solution U(x,t) is available, the maximum error c of the approximation to U(x,tend) is computed and displayed. Here c tend = 15. It should be appreciated that this error has two parts, c one the error made by RKC in the time integration and the other from c the spacial discretization. Some statistics about the integration are c also displayed. c integer neqn,nout parameter (neqn=99,nout=4) integer info(4),idid double precision t,tend,rtol,atol double precision y(neqn),work(8+5*neqn) integer i,next double precision dx,delta,sol,tout(nout),yout(neqn), & truey,error integer nfe,nsteps,naccpt,nrejct,nfesig,maxm common /rkcdid/ nfe,nsteps,naccpt,nrejct,nfesig,maxm external f c------------------------------------------------------------ c Specify the interval of integration in time. The initial c values of the solution at mesh points are provided by the c analytical solution sol(x,t). c------------------------------------------------------------ t = 0d0 tend = 15d0 dx = 10d0/(neqn+1) do 10 i = 1, neqn y(i) = sol(i*dx,t) 10 continue c---------------------------------------------------------------------- c Initialize the output: Define the times at which solutions are to c be reported for plotting and output the number of these times. A c solution is computed on a mesh of equally spaced points in [0, 10]. c Output the number of points in the mesh and then the mesh itself. c Because tout(1) = t, output the initial values for the neqn solution c components along with the values given at 0 and 10. c---------------------------------------------------------------------- delta = (tend - t)/(nout-1) do 20 i = 1,nout tout(i) = t + (i-1)*delta 20 continue open(10,file='exaout') write(10,'(i10)') nout,neqn+2 do 30 i = 0,neqn+1 write(10,'(e10.4)') i*dx 30 continue next = 1 write(10,'(e10.4)') sol(0d0,tout(next)) write(10,'(e10.4)') y write(10,'(e10.4)') sol(10d0,tout(next)) next = next + 1 c-------------------------------------------------- c To compute results at specific times, the code c must return after each step. Common choices for c info(*) have value 0. c info(1) = 0 -- return after each step. c info(2) = 0 -- RKC computes the spectral radius. c info(3) = 0 -- the Jacobian may not be constant. c info(4) = 0 -- ATOL is a scalar. c-------------------------------------------------- info(1) = 0 info(2) = 0 info(3) = 0 info(4) = 0 c------------------------- c Specify the tolerances. c------------------------- rtol = 1d-4 atol = rtol c----------------------------- c Initialize the integration. c----------------------------- idid = 0 c--------------------- c Take a single step: c--------------------- 40 continue call rkc(neqn,f,y,t,tend,rtol,atol,info,work,idid) c------------------------------------------------------------- c Was the step successful? If not, quit with an explanation. c------------------------------------------------------------- if(idid .gt. 2) then write(*,*) ' Failed at t = ',t,' with idid = ',idid stop endif c-------------------------------------------------------------- c To get output at specific points, step towards TEND with RKC c until the integration passes the next output point. Compute c a result at the point using RKCINT. There might be several c output points in the span of a single step by RKC. c-------------------------------------------------------------- 50 continue if(t .ge. tout(next)) then call rkcint(work,tout(next),yout) write(10,'(e10.4)') sol(0d0,tout(next)) write(10,'(e10.4)') yout write(10,'(e10.4)') sol(10d0,tout(next)) next = next + 1 if(next .le. nout) goto 50 endif c-------------------------------------- c Monitor the cost of the integration. c-------------------------------------- if(nsteps .ge. 5000) then write(*,*) ' Quit because of too much work.' endif c------------------------------------- c If not done yet, take another step. c------------------------------------- if(idid .eq. 2) goto 40 c------------------------------------------------------ c Done. Compute the error and report some statistics. c------------------------------------------------------ error = 0d0 do 60 i = 1, neqn truey = sol(i*dx,t) error = max(error,abs(y(i) - truey)) 60 continue write(*,'(/a,d8.1,a,f6.1,a,d8.2)') ' With rtol = atol =',rtol, & ', the maximum error at tend =',tend,' was',error write(*,'(a,i5,a)') ' The integration cost',nfe, & ' function evaluations.' write(*,'(a,i4,a,i3,a)') ' There were',nsteps,' steps (', & nrejct,' rejected).' write(*,'(a,i4)') ' The maximum number of stages used was',maxm end double precision function sol(x,t) c------------------------------------------------------------ c An analytical solution to the reaction-diffusion equation. c------------------------------------------------------------ double precision x,t double precision v,z v = sqrt(0.5d0) z = x - v*t sol = 1d0/(1d0 + exp(v*z)) return end subroutine f(neqn,t,y,dy) c--------------------------------------------------------------- c Semi-discretization of reaction-diffusion equation by central c differences. The analytical solution sol(x,t) is used for c Dirichlet boundary conditions at x = 0 and x = 10. c--------------------------------------------------------------- integer neqn double precision t,y(neqn),dy(neqn) integer i double precision dx,dxsq,sol c dx = 10d0/(neqn+1) dxsq = dx**2 dy(1) = (sol(0d0,t)- 2d0*y(1) + y(2))/dxsq + & (1d0 - y(1))*y(1)**2 do 10 i = 2,neqn-1 dy(i) = (y(i-1) - 2d0*y(i) + y(i+1))/dxsq + & (1d0 - y(i))*y(i)**2 10 continue dy(neqn) = (y(neqn-1)- 2d0*y(neqn) + sol(10d0,t))/dxsq + & (1d0 - y(neqn))*y(neqn)**2 return end double precision function spcrad(neqn,t,y) c-------------------------- c This is a dummy routine. c-------------------------- integer neqn double precision t,y(neqn) spcrad = 0d0 return end