C * * * * * * * * * * * * * * * * * * * * * * * * * C --- DRIVER FOR MEBDF ON THE ANDREW SQUEEZING PROBLEM C * * * * * * * * * * * * * * * * * * * * * * * * * IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (ND=27,LWORK=(41+3*ND)*ND+3,LIWORK=ND+14) DIMENSION Y(ND),WORK(LWORK),IWORK(LIWORK),MBND(4),true(nd) + ,error(nd),massbnd(4) REAL*4 TARRAY(2) external f,pderv,mas write(6,2050) 2050 format(1x,'andrews squeezing problem by mebdf') C --- LOOP FOR DIFFERENT TOLERANCES NTOLMN=4 NTOLMX=10 NTOLDF=2 NRLOOP=(NTOLMX-NTOLMN)*NTOLDF+1 TOLST=0.1D0**NTOLMN TOLFC=0.1D0**(1.D0/NTOLDF) iwork(14)=100000 work(1)=0.0D+0 DO 30 NTOL=1,NRLOOP summ=0.0D+0 sumh=0.0D+0 C --- DIMENSION OF THE SYSTEM N=27 NEQN = N C --- SET DEFAULT VALUES MF=21 INDEX=1 LOUT=6 C --- INITIAL VALUES X=0.0D0 rtol=tolst ATOL=RTOL ITOL=2 MAXDER=6 C --- MAXIMAL NUMBER OF STEPS c IWORK(6)=10000 C --- ENDPOINT OF INTEGRATION y(1 ) = -0.0617138900142764496358948458001d0 y(2 ) = 0d0 y(3 ) = 0.455279819163070380255912382449d0 y(4 ) = 0.222668390165885884674473185609d0 y(5 ) = 0.487364979543842550225598953530d0 y(6 ) = -0.222668390165885884674473185609d0 y(7 ) = 1.23054744454982119249735015568d0 y(8 ) = 0d0 y(9 ) = 0d0 y(10) = 0d0 y(11) = 0d0 y(12) = 0d0 y(13) = 0d0 y(14) = 0d0 y(15) = 14222.4439199541138705911625887d0 y(16) = -10666.8329399655854029433719415d0 y(17) = 0d0 y(18) = 0d0 y(19) = 0d0 y(20) = 0d0 y(21) = 0d0 y(22) = 98.5668703962410896057654982170d0 y(23) = -6.12268834425566265503114393122d0 y(24) = 0d0 y(25) = 0d0 y(26) = 0d0 y(27) = 0d0 H=RTOL CCC CALL DTIME(TARRAY) it1=mclock() XOUT=3.0d-2 xend=xout 220 CONTINUE C --- CALL OF THE SUBROUTINE iwork(1)=7 iwork(2)=7 iwork(3)=13 massbnd(1)=1 massbnd(2)=0 massbnd(3)=0 massbnd(4)=1 mbnd(1)=neqn mbnd(2)=neqn mbnd(3)=neqn mbnd(4)=neqn CALL MEBDF(N,X,H,Y,XOUT,XEND,MF,INDEX,LOUT,LWORK, + WORK,LIWORK,IWORK,MBND,massbnd,MAXDER,ITOL,RTOL,ATOL, + rpar,ipar,f,pderv,mas,ier) IF (INDEX.EQ.1) THEN INDEX=0 GOTO 220 ELSE XOUT=XOUT+1.D0 END IF sum=0.0D+0 do 566 i=1,n true(i)=y(i) 566 continue y( 1) = 0.1581077d+2 y( 2) = -0.1575637d+2 y( 3) = 0.4082224d-1 y( 4) = -0.5347301d+0 y( 5) = 0.5244100d+0 y( 6) = 0.5347301d+0 y( 7) = 0.1048081d+1 y( 8) = 0.1139920d+4 y( 9) = -0.1424379d+4 y(10) = 0.1103278d+2 y(11) = 0.1929315d+2 y(12) = 0.5735633d+0 y(13) = -0.1929315d+2 y(14) = 0.3231753d+0 y(15) = -0.2463131d+5 y(16) = 0.5184964d+5 y(17) = 0.3241026d+6 y(18) = 0.5667495d+6 y(19) = 0.1674364d+5 y(20) = -0.5667495d+6 y(21) = 0.9826496d+4 y(22) = 0.1991753d+3 y(23) = -0.2975527d+2 y(24) = 0.2306654d+2 y(25) = 0.3145273d+2 y(26) = 0.2264249d+2 y(27) = 0.1161737d+2 do 270 k=1,n error(k)=dabs(true(k)-y(k)) sum=sum+((error(k)/(rtol*dabs(true(k))+atol))**2) 270 continue sum=dsqrt(sum/dble(nd)) 20 continue it2=mclock() tarray(1)=(it2-it1)/100.0D+0 summ=sum*atol us=-log10(summ) write(6,9358) us 9358 format(1x,'number of correct digits =',f8.4) NFE=IWORK(7) NJE=IWORK(8) NSTEP=IWORK(5) NDEC=IWORK(9) NBSOL=IWORK(10) NFAIL=IWORK(6) WRITE(6,901)NFE,NJE,NSTEP,NDEC,NBSOL,NFAIL 901 format(1x,6i12) WRITE(6,*)' ***** TOL=',RTOL,' ELAPSED TIME=',TARRAY(1),' ****' C -------- NEW TOLERANCE --- 25 TOLST=TOLST*TOLFC 30 CONTINUE STOP END c----------------------------------------------------------------------- subroutine f(neqn,t,y,dy,ipar,rpar,ier) integer neqn double precision t,y(neqn),dy(neqn),ipar(*),rpar(*) integer i,j double precision m1,m2,m3,m4,m5,m6,m7,xa,ya,xb,yb,xc,yc,c0, + i1,i2,i3,i4,i5,i6,i7,d,da,e,ea,rr,ra,l0, + ss,sa,sb,sc,sd,ta,tb,u,ua,ub,zf,zt,fa,mom, + sibe,sith,siga,siph,side,siom,siep, + cobe,coth,coga,coph,code,coom,coep, + sibeth,siphde,siomep,cobeth,cophde,coomep, + bep,thp,php,dep,omp,epp, + m(7,7),ff(7),gp(6,7),g(6),xd,yd,lang,force,fx,fy parameter (m1=.04325d0,m2=.00365d0,m3=.02373d0,m4=.00706d0, + m5=.07050d0,m6=.00706d0,m7=.05498d0, + xa=-.06934d0,ya=-.00227d0, + xb=-0.03635d0,yb=.03273d0, + xc=.014d0,yc=.072d0,c0=4530d0, + i1=2.194d-6,i2=4.410d-7,i3=5.255d-6,i4=5.667d-7, + i5=1.169d-5,i6=5.667d-7,i7=1.912d-5, + d=28d-3,da=115d-4,e=2d-2,ea=1421d-5, + rr=7d-3,ra=92d-5,l0=7785d-5, + ss=35d-3,sa=1874d-5,sb=1043d-5,sc=18d-3,sd=2d-2, + ta=2308d-5,tb=916d-5,u=4d-2,ua=1228d-5,ub=449d-5, + zf=2d-2,zt=4d-2,fa=1421d-5,mom=33d-3) sibe = sin(y(1)) sith = sin(y(2)) siga = sin(y(3)) siph = sin(y(4)) side = sin(y(5)) siom = sin(y(6)) siep = sin(y(7)) c cobe = cos(y(1)) coth = cos(y(2)) coga = cos(y(3)) coph = cos(y(4)) code = cos(y(5)) coom = cos(y(6)) coep = cos(y(7)) c sibeth = sin(y(1)+y(2)) siphde = sin(y(4)+y(5)) siomep = sin(y(6)+y(7)) c cobeth = cos(y(1)+y(2)) cophde = cos(y(4)+y(5)) coomep = cos(y(6)+y(7)) c bep = y(8) thp = y(9) php = y(11) dep = y(12) omp = y(13) epp = y(14) c do 20 j = 1,7 do 10 i = 1,7 m(i,j) = 0d0 10 continue 20 continue c m(1,1) = m1*ra**2 + m2*(rr**2-2*da*rr*coth+da**2) + i1 + i2 m(2,1) = m2*(da**2-da*rr*coth) + i2 m(2,2) = m2*da**2 + i2 m(3,3) = m3*(sa**2+sb**2) + i3 m(4,4) = m4*(e-ea)**2 + i4 m(5,4) = m4*((e-ea)**2+zt*(e-ea)*siph) + i4 m(5,5) = m4*(zt**2+2*zt*(e-ea)*siph+(e-ea)**2) + m5*(ta**2+tb**2) + + i4 + i5 m(6,6) = m6*(zf-fa)**2 + i6 m(7,6) = m6*((zf-fa)**2-u*(zf-fa)*siom) + i6 m(7,7) = m6*((zf-fa)**2-2*u*(zf-fa)*siom+u**2) + m7*(ua**2+ub**2) + + i6 + i7 do 40 j=2,7 do 30 i=1,j-1 m(i,j) = m(j,i) 30 continue 40 continue c xd = sd*coga + sc*siga + xb yd = sd*siga - sc*coga + yb lang = sqrt ((xd-xc)**2 + (yd-yc)**2) force = - c0 * (lang - l0)/lang fx = force * (xd-xc) fy = force * (yd-yc) ff(1) = mom - m2*da*rr*thp*(thp+2*bep)*sith ff(2) = m2*da*rr*bep**2*sith ff(3) = fx*(sc*coga - sd*siga) + fy*(sd*coga + sc*siga) ff(4) = m4*zt*(e-ea)*dep**2*coph ff(5) = - m4*zt*(e-ea)*php*(php+2*dep)*coph ff(6) = - m6*u*(zf-fa)*epp**2*coom ff(7) = m6*u*(zf-fa)*omp*(omp+2*epp)*coom c do 60 j=1,7 do 50 i=1,6 gp(i,j) = 0d0 50 continue 60 continue c gp(1,1) = - rr*sibe + d*sibeth gp(1,2) = d*sibeth gp(1,3) = - ss*coga gp(2,1) = rr*cobe - d*cobeth gp(2,2) = - d*cobeth gp(2,3) = - ss*siga gp(3,1) = - rr*sibe + d*sibeth gp(3,2) = d*sibeth gp(3,4) = - e*cophde gp(3,5) = - e*cophde + zt*side gp(4,1) = rr*cobe - d*cobeth gp(4,2) = - d*cobeth gp(4,4) = - e*siphde gp(4,5) = - e*siphde - zt*code gp(5,1) = - rr*sibe + d*sibeth gp(5,2) = d*sibeth gp(5,6) = zf*siomep gp(5,7) = zf*siomep - u*coep gp(6,1) = rr*cobe - d*cobeth gp(6,2) = - d*cobeth gp(6,6) = - zf*coomep gp(6,7) = - zf*coomep - u*siep c g(1) = rr*cobe - d*cobeth - ss*siga - xb g(2) = rr*sibe - d*sibeth + ss*coga - yb g(3) = rr*cobe - d*cobeth - e*siphde - zt*code - xa g(4) = rr*sibe - d*sibeth + e*cophde - zt*side - ya g(5) = rr*cobe - d*cobeth - zf*coomep - u*siep - xa g(6) = rr*sibe - d*sibeth - zf*siomep + u*coep - ya c do 70 i=1,14 dy(i) = y(i+7) 70 continue do 100 i=15,21 dy(i) = -ff(i-14) do 80 j=1,7 dy(i) = dy(i)+m(i-14,j)*y(j+14) 80 continue do 90 j=1,6 dy(i) = dy(i)+gp(j,i-14)*y(j+21) 90 continue 100 continue do 110 i=22,27 dy(i) = g(i-21) 110 continue return end c----------------------------------------------------------------------- subroutine pderv(t,y,jac,neqn,ldim,ipar,rpar,ier) integer neqn,ldim double precision t,y(neqn),jac(ldim,neqn),ipar(*),rpar(*) c----------------------------------------------------------------------- c the Jacobian computed here is an approximation, see p. 540 of c Hairer & Wanner `solving ordinary differential equations II' c----------------------------------------------------------------------- integer i,j double precision m1,m2,m3,m4,m5,m6,m7, + i1,i2,i3,i4,i5,i6,i7,d,da,e,ea,rr,ra, + ss,sa,sb,ta,tb,u,ua,ub,zf,zt,fa, + sibe,siga,siph,side,siom,siep, + cobe,coth,coga,code,coep, + sibeth,siphde,siomep,cobeth,cophde,coomep, + m(7,7),gp(6,7) parameter (m1=.04325d0,m2=.00365d0,m3=.02373d0,m4=.00706d0, + m5=.07050d0,m6=.00706d0,m7=.05498d0, + i1=2.194d-6,i2=4.410d-7,i3=5.255d-6,i4=5.667d-7, + i5=1.169d-5,i6=5.667d-7,i7=1.912d-5, + d=28d-3,da=115d-4,e=2d-2,ea=1421d-5, + rr=7d-3,ra=92d-5, + ss=35d-3,sa=1874d-5,sb=1043d-5, + ta=2308d-5,tb=916d-5,u=4d-2,ua=1228d-5,ub=449d-5, + zf=2d-2,zt=4d-2,fa=1421d-5) sibe = sin(y(1)) siga = sin(y(3)) siph = sin(y(4)) side = sin(y(5)) siom = sin(y(6)) siep = sin(y(7)) c cobe = cos(y(1)) coth = cos(y(2)) coga = cos(y(3)) code = cos(y(5)) coep = cos(y(7)) c sibeth = sin(y(1)+y(2)) siphde = sin(y(4)+y(5)) siomep = sin(y(6)+y(7)) c cobeth = cos(y(1)+y(2)) cophde = cos(y(4)+y(5)) coomep = cos(y(6)+y(7)) c do 20 j = 1,7 do 10 i = 1,7 m(i,j) = 0d0 10 continue 20 continue c m(1,1) = m1*ra**2 + m2*(rr**2-2*da*rr*coth+da**2) + i1 +i2 m(2,1) = m2*(da**2-da*rr*coth) + i2 m(2,2) = m2*da**2 + i2 m(3,3) = m3*(sa**2+sb**2) + i3 m(4,4) = m4*(e-ea)**2 + i4 m(5,4) = m4*((e-ea)**2+zt*(e-ea)*siph) + i4 m(5,5) = m4*(zt**2+2*zt*(e-ea)*siph+(e-ea)**2) + m5*(ta**2+tb**2) + + i4 + i5 m(6,6) = m6*(zf-fa)**2 + i6 m(7,6) = m6*((zf-fa)**2-u*(zf-fa)*siom) + i6 m(7,7) = m6*((zf-fa)**2-2*u*(zf-fa)*siom+u**2) + m7*(ua**2+ub**2) + + i6 + i7 do 40 j=2,7 do 30 i=1,j-1 m(i,j) = m(j,i) 30 continue 40 continue c do 60 j=1,7 do 50 i=1,6 gp(i,j) = 0d0 50 continue 60 continue c gp(1,1) = - rr*sibe + d*sibeth gp(1,2) = d*sibeth gp(1,3) = - ss*coga gp(2,1) = rr*cobe - d*cobeth gp(2,2) = - d*cobeth gp(2,3) = - ss*siga gp(3,1) = - rr*sibe + d*sibeth gp(3,2) = d*sibeth gp(3,4) = - e*cophde gp(3,5) = - e*cophde + zt*side gp(4,1) = rr*cobe - d*cobeth gp(4,2) = - d*cobeth gp(4,4) = - e*siphde gp(4,5) = - e*siphde - zt*code gp(5,1) = - rr*sibe + d*sibeth gp(5,2) = d*sibeth gp(5,6) = zf*siomep gp(5,7) = zf*siomep - u*coep gp(6,1) = rr*cobe - d*cobeth gp(6,2) = - d*cobeth gp(6,6) = - zf*coomep gp(6,7) = - zf*coomep - u*siep c do 80 j=1,neqn do 70 i=1,neqn jac(i,j) = 0d0 70 continue 80 continue do 90 i=1,14 jac(i,i+7) = 1d0 90 continue do 110 i=1,7 do 100 j=1,7 jac(14+j,14+i) = m(j,i) 100 continue 110 continue do 130 i=1,6 do 120 j=1,7 jac(14+j,21+i) = gp(i,j) 120 continue 130 continue do 150 i=1,7 do 140 j=1,6 jac(21+j,i) = gp(j,i) 140 continue 150 continue return end c----------------------------------------------------------------------- c----------------------------------------------------------------------- subroutine mas(neqn,am,ldim,ipar,rpar,ier) integer neqn,ldim double precision am(ldim,neqn),ipar(*),rpar(*) integer i do 10 i=1,neqn am(1,i) = 0d0 10 continue do 20 i=1,14 am(1,i) = 1d0 20 continue return end c----------------------------------------------------------------------- c-----------------------------------------------------------------------