c----------------------------------------------------------------------- c----------------------------------------------------------------------- c c This file is part of the Test Set for IVP solvers c http://www.dm.uniba.it/~testset/ c c Car Axis problem (in index 3 formulation) c index 3 DAE of dimension 10 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/caraxis.f c c This is revision c $Id: caraxis.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 = 'Car Axis problem' problm = 'caraxis' type = 'DAE' neqn = 10 ndisc = 0 t(0) = 0d0 t(1) = 3d0 numjac = .true. mljac = neqn mujac = neqn mlmas = 0 mumas = 0 do 10 i=1,4 ind(i) = 1 10 continue do 20 i=5,8 ind(i) = 2 20 continue do 30 i=9,10 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 i double precision M,eps,L,L0, + xl,yl,xr,yr,xla,yla,xra,yra,lam1,lam2 integer ierr,ipar double precision rpar M = 10d0 eps = 1d-2 L = 1d0 L0 = 0.5d0 xr = L xl = 0 yr = L0 yl = yr xra = -L0/L xla = xra yra = 0d0 yla = 0d0 lam1 = 0d0 lam2 = 0d0 y(1) = xl y(2) = yl y(3) = xr y(4) = yr y(5) = xla y(6) = yla y(7) = xra y(8) = yra y(9) = lam1 y(10) = lam2 call feval(neqn,0d0,y,y,yprime,ierr,rpar,ipar) do 10 i=1,4 yprime(i) = y(i+4) yprime(i+4)= 2d0/(M*eps*eps)*yprime(i+4) 10 continue do 20 i=9,10 yprime(i)=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 = 4 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,f,ierr,rpar,ipar) integer neqn,ierr,ipar(*) double precision t,y(neqn),yprime(neqn),f(neqn),rpar(*) integer i double precision M,L,L0,w,r,xb,yb,Ll,Lr,eps,g, + xl,yl,xr,yr,lam1,lam2 eps = 1d-2 M = 10d0 L = 1d0 L0 = 0.5d0 r = 0.1d0 w = 10d0 g = 1d0 yb = r*sin(w*t) xb = sqrt(L*L-yb*yb) do 10 i=1,4 f(i) = y(i+4) 10 continue xl = y(1) yl = y(2) xr = y(3) yr = y(4) lam1 = y(9) lam2 = y(10) Ll = sqrt(xl**2+yl**2) Lr = sqrt((xr-xb)**2+(yr-yb)**2) f(5) =(L0-Ll)*xl/Ll +lam1*xb+2d0*lam2*(xl-xr) f(6) =(L0-Ll)*yl/Ll +lam1*yb+2d0*lam2*(yl-yr)-M*eps*eps*g/2d0 f(7) =(L0-Lr)*(xr-xb)/Lr -2d0*lam2*(xl-xr) f(8) =(L0-Lr)*(yr-yb)/Lr -2d0*lam2*(yl-yr)-M*eps*eps*g/2d0 f(9) = xb*xl+yb*yl f(10) = (xl-xr)**2+(yl-yr)**2-L*L 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 dummy subroutine c 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 double precision M,eps M = 10d0 eps = 1d-2 do 10 i=1,4 dfddy(1,i) = 1d0 10 continue do 20 i=5,8 dfddy(1,i) = M*eps*eps/2d0 20 continue do 30 i=9,10 dfddy(1,i) = 0d0 30 continue return end c----------------------------------------------------------------------- subroutine solut(neqn,t,y) integer neqn double precision t,y(neqn) c c computed using quadruple precision GAMD on an c Alphaserver DS20E, with a 667 MHz EV67 processor. c c rtol = atol = h0 = 1.0d-24 c y( 1) = 0.493455784275402809122d-1 y( 2) = 0.496989460230171153861d0 y( 3) = 0.104174252488542151681d1 y( 4) = 0.373911027265361256927d0 y( 5) = -0.770583684040972357970d-1 y( 6) = 0.744686658723778553466d-2 y( 7) = 0.175568157537232222276d-1 y( 8) = 0.770341043779251976443d0 y( 9) = -0.473688659084893324729d-2 y( 10) = -0.110468033125734368808d-2 return end