c----------------------------------------------------------------------- c----------------------------------------------------------------------- c c This file is part of the Test Set for IVP solvers c http://www.dm.uniba.it/~testset/ c c Andrews'' squeezing mechanism (in index 3 formulation) c index 3 DAE of dimension 27 c c DISCLAIMER: see c http://www.dm.uniba.it/~testset/disclaimer.php c c The most recent version of this source file can be found at c http://www.dm.uniba.it/~testset/src/problems/andrews.f c c This is revision c $Id: andrews.F,v 1.2 2006/10/02 10:29:13 testset Exp $ c c----------------------------------------------------------------------- integer function pidate() pidate = 20060828 return end c----------------------------------------------------------------------- subroutine prob(fullnm,problm,type, + neqn,ndisc,t, + numjac,mljac,mujac, + nummas,mlmas,mumas, + ind) character*(*) fullnm, problm, type integer neqn,ndisc,mljac,mujac,mlmas,mumas,ind(*) double precision t(0:*) logical numjac, nummas integer i fullnm = 'Andrews'' squeezing mechanism' problm = 'andrews' type = 'DAE' neqn = 27 ndisc = 0 t(0) = 0d0 t(1) = 3d-2 numjac = .false. mljac = neqn mujac = neqn mlmas = 0 mumas = 0 do 10 i=1,7 ind(i) = 1 10 continue do 20 i=8,14 ind(i) = 2 20 continue do 30 i=15,27 ind(i) = 3 30 continue return end c----------------------------------------------------------------------- subroutine init(neqn,t,y,yprime,consis) integer neqn double precision t,y(neqn),yprime(neqn) logical consis integer k 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 do 10 k=1,14 yprime(k) = y(k+7) 10 continue do 20 k=15,27 yprime(k) = 0d0 20 continue consis = .true. return end c----------------------------------------------------------------------- subroutine settolerances(neqn,rtol,atol,tolvec) integer neqn logical tolvec double precision rtol(neqn), atol(neqn) tolvec = .false. return end c----------------------------------------------------------------------- subroutine setoutput(neqn,solref,printsolout, + nindsol,indsol) logical solref, printsolout integer neqn, nindsol integer indsol(neqn) c the reference solution is available solref = .true. c output file is required printsolout = .true. c default values if printsolout is .true. nindsol = 7 c only nindsol component of indsol are referenced do i=1,nindsol indsol(i) = i end do return end c----------------------------------------------------------------------- subroutine feval(neqn,t,y,yprime,ff,ierr,rpar,ipar) integer neqn,ierr,ipar(*) double precision t,y(neqn),yprime(neqn),ff(neqn),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),f(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) f(1) = mom - m2*da*rr*thp*(thp+2*bep)*sith f(2) = m2*da*rr*bep**2*sith f(3) = fx*(sc*coga - sd*siga) + fy*(sd*coga + sc*siga) f(4) = m4*zt*(e-ea)*dep**2*coph f(5) = - m4*zt*(e-ea)*php*(php+2*dep)*coph f(6) = - m6*u*(zf-fa)*epp**2*coom f(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 ff(i) = y(i+7) 70 continue do 100 i=15,21 ff(i) = -f(i-14) do 80 j=1,7 ff(i) = ff(i)+m(i-14,j)*y(j+14) 80 continue do 90 j=1,6 ff(i) = ff(i)+gp(j,i-14)*y(j+21) 90 continue 100 continue do 110 i=22,27 ff(i) = g(i-21) 110 continue return end c----------------------------------------------------------------------- subroutine jeval(ldim,neqn,t,y,yprime,dfdy,ierr,rpar,ipar) integer ldim,neqn,ierr,ipar(*) double precision t,y(neqn),yprime(neqn),dfdy(ldim,neqn),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 dfdy(i,j) = 0d0 70 continue 80 continue do 90 i=1,14 dfdy(i,i+7) = 1d0 90 continue do 110 i=1,7 do 100 j=1,7 dfdy(14+j,14+i) = m(j,i) 100 continue 110 continue do 130 i=1,6 do 120 j=1,7 dfdy(14+j,21+i) = gp(i,j) 120 continue 130 continue do 150 i=1,7 do 140 j=1,6 dfdy(21+j,i) = gp(j,i) 140 continue 150 continue return end c----------------------------------------------------------------------- subroutine meval(ldim,neqn,t,y,yprime,dfddy,ierr,rpar,ipar) integer ldim,neqn,ierr,ipar(*) double precision t,y(neqn),yprime(neqn),dfddy(ldim,neqn),rpar(*) integer i do 10 i=1,neqn dfddy(1,i) = 0d0 10 continue do 20 i=1,14 dfddy(1,i) = 1d0 20 continue return end c----------------------------------------------------------------------- subroutine solut(neqn,t,y) integer neqn double precision t,y(neqn) c c computed at Cray C90, using Cray double precision: c Solving Andrews` squeezing mechanism using PSIDE c User input: c c give relative error tolerance: 1d-14 c give absolute error tolerance: 1d-14 c c Integration characteristics: c c number of integration steps 1112 c number of accepted steps 1022 c number of f evaluations 26000 c number of Jacobian evaluations 173 c number of LU decompositions 2056 c c CPU-time used: 56.74 sec y( 1) = 0.1581077119629904d+2 y( 2) = -0.1575637105984298d+2 y( 3) = 0.4082224013073101d-1 y( 4) = -0.5347301163226948d+0 y( 5) = 0.5244099658805304d+0 y( 6) = 0.5347301163226948d+0 y( 7) = 0.1048080741042263d+1 y( 8) = 0.1139920302151208d+4 y( 9) = -0.1424379294994111d+4 y( 10) = 0.1103291221937134d+2 y( 11) = 0.1929337464421385d+2 y( 12) = 0.5735699284790808d+0 y( 13) = -0.1929337464421385d+2 y( 14) = 0.3231791658026955d+0 y( 15) = -0.2463176316945196d+5 y( 16) = 0.5185037701610329d+5 y( 17) = 0.3241025686413781d+6 y( 18) = 0.5667493645176213d+6 y( 19) = 0.1674362929479361d+5 y( 20) = -0.5667493645176222d+6 y( 21) = 0.9826520791458422d+4 y( 22) = 0.1991753333731910d+3 y( 23) = -0.2975531228015052d+2 y( 24) = 0.2306654119098399d+2 y( 25) = 0.3145271365475927d+2 y( 26) = 0.2264249232082739d+2 y( 27) = 0.1161740700019673d+2 return end