c * * * * * * * * * * * * * * * * * * * * * * * * * c This example is slightly modified version of example presented c by Assyr Abdulle and described in the book of c G.Wanner and E.Hairer Solving Ordinary Differential Equations II c * * * * * * * * * * * * * * * * * * * * * * * * * c c This driver shows how to use DUMKA3. It solves a c system of ODEs resulting from the 2-dimensional space c discretization of the Brusselator equations (u=u(x,y,t),v=v(x,y,t)): c c u_t=1+u^2*v-4.4*u+0.1*(u_{xx}+u_{yy})+f(x,y,t) c v_t=3.4*u-u^2*v+0.1*(v_{xx}+v_{yy}) for t>=0, 0<= x <= 1, 0<= y <= 1 c c with initial conditions c c u(x,y,0)=22*y*(1-y)^{3/2} v(x,y,0)=27*x*(1-x)^{3/2} c c and periodic boundary conditions c c u(x+1,y,t)=u(x,y,t), v(x,y+1,t)=v(x,y,t). c c The function f is defined by (inhomogeneity) c c f(x,y,t)=5 if (x-0.3)^2+(y-0.6)^2<= 0.1^2 and t>=1.1 c =0 else c c We discretize the space variables with c x_i=i/(N+1), y_i=i/(N+1) for i=0,1,...,N, c with N=128. We obtain a system of 32768 c equations. The spectral radius of the Jacobian c can be estimated with the Gershgorin theorem c (13200 is an estimation for it). Thus we c provide an external function RHO, giving c the spectral radius of the Jacobian matrix. c As output point we choose t_out=1.5 c c-------------------------------------------------------- include 'dumka3.f' implicit double precision (a-h,o-z) parameter (nsd=128,neqn=nsd*nsd*2) dimension y(nsd,nsd,2),z(nsd,nsd,2),z1(nsd,nsd,2),z2(nsd,nsd,2), & z3(nsd,nsd,2),z4(nsd,nsd,2),atol(1),rtol(1) integer iwork(12) external fbrus c --- common parameters for the problem ----- common/trans/alf,ns,nssq,nsm1sq c ----- file for solution ----- rewind 8 c ----- dimensions ----- ns=nsd nssq=ns*ns n=2*nssq alf=1.0d-1 alph=alf c-------------------------------------------------------- c Initialize iwork: c iwork(1)=0 DUMKA3 (!try to!) calculate spectral radius internally c (if iwork(1)=1 - rho is used to determine spectral radius) c iwork(2)=0 The Jacobian is not constant c (if iwork(2)=1 - rho is called only once) c iwork(3)=0 Calculate problem until t=tend c (if iwork(3)=1 - dumka exits after the first step) c iwork(4)=0 Atol and rtol are scalars. c-------------------------------------------------------- iwork(1)=0 iwork(2)=0 iwork(3)=0 iwork(4)=0 c ----- initial and end point of integration ----- t=0.0d0 tend=1.5d0 c ----- initial values ----- ans=ns do j=1,ns yy=(j-1)/ans do i=1,ns xx=(i-1)/ans y(j,i,1)=22.d0*yy*(1.d0-yy)**(1.5d0) y(j,i,2)=27.d0*xx*(1.d0-xx)**(1.5d0) end do end do c ----- required tolerance ----- tol = 0.01d0 rtol(1)=tol atol(1)=tol c ----- initial step size ----- h=1.0d-4 c ----- integration ----- write(6,*) 'integration of the 2-dim brusselator problem' call dumka3(neqn,t,tend,h,fbrus,y,z,z1,z2,z3, & z4,atol,rtol,iwork) open(8,file='solDumkaU') open(9,file='solDumkaV') c ----- print solution ----- do j=1,ns yy=(j-1)/ans do i=1,ns xx=(i-1)/ans write (8,300) xx,yy,y(j,i,1) write (9,300) xx,yy,y(j,i,2) end do end do close (8) close(9) c ----- print statistics ----- write(6,*) 'Solution is tabulated in file sol.out' write(6,*) 'Max estimation of the spectral radius=',iwork(11) write(6,*) 'Min estimation of the spectral radius=',iwork(12) write(6,*) 'Max number of stages used=',iwork(10) write(6,*) 'Number of f eval. for the spectr. radius=',iwork(9) write (6,91) iwork(5),iwork(6),iwork(7),iwork(8) 91 format(' Number of f evaluations=',i5,' steps=',i4, & ' accpt=',i4,' rejct=',i3) 300 format(1f22.16,' ',1f22.16,' ',1f22.16) c-------------------------------------------------------- c End of main program c-------------------------------------------------------- end c-------------------------------------------------------- c The subroutine RHO gives an estimation of the spectral c radius of the Jacobian matrix of the problem. This c is a bound for the whole interval and thus RHO is called c once. c-------------------------------------------------------- double precision function rho(neqn,t,y) implicit double precision (a-h,o-z) dimension y(ns,ns,2) common/trans/alf,ns,nssq,nsm1sq rho = 8.0d0*nssq*alf + 2.d0 return end c-------------------------------------------------------- c The subroutine fun compute the value of f(x,y) and c has to be declared as external. c-------------------------------------------------------- subroutine fbrus(neqn,x,y,f) c ----- brusselator with diffusion in 2 dim. space ----- implicit double precision (a-h,o-z) dimension y(ns,ns,2),f(ns,ns,2) common/trans/alf,ns,nssq,nsm1sq c ----- constants for inhomogenity ----- ans=ns radsq=0.1d0**2 if(x.ge.1.1d0)then bet=5.00d0 else bet=0.00d0 end if c ----- big loop ----- do j=1,ns do i=1,ns c ----- left neighbour ----- if(i.eq.1)then uleft=y(j,ns,1) vleft=y(j,ns,2) else uleft=y(j,i-1,1) vleft=y(j,i-1,2) end if c ----- right neighbour ----- if(i.eq.ns)then uright=y(j,1,1) vright=y(j,1,2) else uright=y(j,i+1,1) vright=y(j,i+1,2) end if c ----- lower neighbour ----- if(j.eq.1)then ulow=y(ns,i,1) vlow=y(ns,i,2) else ulow=y(j-1,i,1) vlow=y(j-1,i,2) end if c ----- upper neighbour ----- if(j.eq.ns)then uup=y(1,i,1) vup=y(1,i,2) else uup=y(j+1,i,1) vup=y(j+1,i,2) end if c ----- the derivative ----- uij=y(j,i,1) vij=y(j,i,2) f(j,i,1)=1.d0+uij*uij*vij-4.4d0*uij & +alf*nssq*(uleft+uright+ulow+uup-4.d0*uij) f(j,i,2)=3.4d0*uij - uij*uij*vij & +alf*nssq*(vleft+vright+vlow+vup-4.d0*vij) c ----- inhomogenity ----- yy=(j-1)/ans xx=(i-1)/ans if(((xx-0.3d0)**2+(yy-0.6d0)**2).le.radsq)then f(j,i,1)=f(j,i,1)+bet end if end do end do return end