c c This file is part of the Test Set for BVP solvers c http://www.dm.uniba.it/~bvpsolvers/ c c Problem MEASLES c BVP of dimension 6 c c DISCLAIMER: see c http://www.dm.uniba.it/~bvpsolvers/testsetbvpsolvers/?page_id=34 c c The most recent version of this source file can be found at c http://www.dm.uniba.it/~bvpsolvers/bvpTestSet/fortransrc/problems/Measles.f c c c c----------------------------------------------------------------------- c integer function pidate() pidate = 20130531 return end c----------------------------------------------------------------------- subroutine prob(fullnm,problm,type, + neqn,aleft,aright,nlbc,ms, + numjac,numbcjac,linear,rpar,ipar) implicit none character*(25), intent(out) :: fullnm character*(15), intent(out) ::problm character*(5), intent(out) :: type integer, intent(out) :: neqn,ipar(1), nlbc, ms(1:6) double precision, intent(out) :: aleft, aright, rpar(4) logical, intent(out) :: numjac, numbcjac,linear double precision :: pi fullnm = ' ' problm = ' ' type = ' ' fullnm = 'Problem Measles' problm = 'Measles' type = 'BVP' ms(1:6)= 1 neqn = 6 nlbc = 3 aleft = 0.0d0 aright = 1.0d0 numjac = .false. numbcjac = .false. linear = .false. rpar(1) = 0.5 PI=4.0D0*DATAN(1.0D0) rpar(2) = pi ipar(1) = 0 return end c c subroutine settolerances(neqn,rtol,atol,tolvec,ntol,ltol) integer, intent(in) :: neqn integer, intent(out) :: ntol, ltol(*) logical, intent(out) :: tolvec double precision, intent(in out) :: rtol(neqn), atol(neqn) tolvec = .true. ntol = neqn ltol(1)=1 DO i=2,ntol ltol(i) = i rtol(i) = rtol(1) atol(i) = atol(1) ENDDO return end c--------------------------------------------------------------------- subroutine setoutput(neqn,solref,printsolout, + nindsol,indsol) logical, intent(out) :: solref, printsolout integer, intent(in) :: neqn integer, intent(out) :: nindsol integer, intent(out) :: indsol(*) c the reference solution is available solref = .false. c output file is required printsolout = .true. c default values if printsolout is .true. nindsol = neqn c only nindsol component of indsol are referenced do i=1,nindsol indsol(i) = i end do return end c----------------------------------------------------------------------- subroutine init(neqn,aleft,aright,ms,yval0, nmsh, givmsh,givey, + xguess,yguess,dmguess,nudim) integer, intent(in) :: neqn, ms(1:neqn) double precision, intent(in) :: aleft,aright integer, intent(out) :: nmsh double precision, intent(out) :: xguess(*),yguess(nudim,*) double precision, intent(out) :: dmguess(nudim,*) double precision, intent(out) :: yval0(1:neqn) logical, intent(out) :: givmsh, givey double precision :: h givmsh = .true. givey = .true. yval0(1:neqn) = 0.5d0 nmsh = 3*5+1 c---- dummy values, not used if givmsh=FALSE, givey=false h = (aright-aleft)/(nmsh-1) xguess(1) = aleft do i=2,nmsh xguess(i)=xguess(i-1)+h end do xguess(nmsh)= aright DO i=1,neqn yguess(i,1:nmsh) = yval0(i) end do DO i=1,neqn dmguess(i,1:nmsh) = 0.0d0 end do return end c----------------------------------------------------------------------- subroutine feval(neqn,t,y,f,rpar,ipar) integer, intent(in) :: neqn,ipar(*) double precision, intent(in):: t,y(1:neqn),rpar(*) double precision, intent(out):: f(1:neqn) double precision :: mu, lambda, eta, B0, bt,pigreco mu=0.02d0 lambda=0.0279d0 eta=0.01d0 B0=1575d0 pigreco = rpar(2) bt=B0*(1.0d0+(COS(2*pigreco*t))) F(1)=mu-(bt*Y(1)*Y(3)) F(2)=(bt*Y(1)*Y(3))-(Y(2)/lambda) F(3)=(Y(2)/lambda)-(Y(3)/eta) F(4)=0 F(5)=0 F(6)=0 return end c----------------------------------------------------------------------- subroutine bceval(i,neqn,y,bc,rpar,ipar) integer, INTENT(IN) :: i, neqn, ipar(*) double precision, INTENT(IN) :: y(neqn), rpar(*) double precision, INTENT(OUT) :: bc if (i .eq. 1) bc=Y(1)-Y(4) if (i .eq. 2) bc=Y(2)-Y(5) if (i .eq. 3) bc=Y(3)-Y(6) if (i .eq. 4) bc=Y(1)-Y(4) if (i .eq. 5) bc=Y(2)-Y(5) if (i .eq. 6) bc=Y(3)-Y(6) return end SUBROUTINE jeval(NCOMP,X,Z,DF,RPAR,IPAR) IMPLICIT NONE INTEGER NCOMP, IPAR(*), I, J DOUBLE PRECISION X, Z, DF, RPAR(*) DIMENSION Z(*),DF(NCOMP,*) double precision :: mu, lambda, eta, B0, bt,pigreco mu=0.02d0 lambda=0.0279d0 eta=0.01d0 B0=1575d0 pigreco = rpar(2) bt=B0*(1.0d0+(COS(2*pigreco*X))) DF(1:ncomp,1:ncomp)=0.0d0 DF(1,1)=-bt*Z(3) DF(1,3)=-bt*Z(1) DF(2,1)= bt*Z(3) DF(2,2)=-1.0d0/lambda DF(2,3)= bt*Z(1) DF(3,2)= 1.0d0/lambda DF(3,3)= -1.0d0/eta RETURN END c---------------------------------------------------------------------- c computation of the jacobian for the boundary conditions c c SUBROUTINE dbceval(I,NCOMP,Z,DG,RPAR,IPAR) IMPLICIT NONE INTEGER I, NCOMP, IPAR(*) DOUBLE PRECISION Z, DG, RPAR(*) DIMENSION Z(*),DG(*) DG(1)=0.D0 DG(2)=0.D0 DG(3)=0.D0 DG(4)=0.D0 DG(5)=0.D0 DG(6)=0.D0 IF (I.EQ.1 .OR. I .EQ. 4) THEN DG(1)= 1.0d0 DG(4)= -1.0d0 ENDIF IF (I.EQ.2.OR.I.EQ.5) THEN DG(2)= 1.0d0 DG(5)=-1.0d0 ENDIF IF (I.EQ.3.OR.I.EQ.6) THEN DG(3)= 1.0d0 DG(6)=-1.0d0 ENDIF RETURN END c----------------------------------------------------------------------- subroutine solut(neqn,x,y,nmsh,rpar,ipar) integer, intent(in) :: neqn, ipar(*) double precision, intent(in) :: x(nmsh),rpar(*) double precision, intent(out) :: y(1:neqn,1:nmsh) double precision lambda, lambda2 return end