c c This file is part of the Test Set for BVP solvers c http://www.dm.uniba.it/~bvpsolvers/ c c Problem bvpT4 c ODE of dimension 2 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/bvpT4.f cc 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(2) double precision, intent(out) :: aleft, aright, rpar(2) logical, intent(out) :: numjac, numbcjac,linear double precision :: EPS common /pars/ EPS c c For Singularly perturbed BVP (SPBVP) eps should be the c first elements of parms and the first element of rpar c fullnm = ' ' problm = ' ' type = ' ' fullnm = 'Problem bvpT4' problm = 'bvpT4' type = 'SPBVP' ms(1) = 1 ms(2) = 1 neqn = 2 nlbc = 1 aleft = -1.0d0 aright = 1.0d0 numjac = .false. numbcjac = .false. linear = .false. c c EPS is an input parameter in report.f c c rpar(1) may contains epsin, the starting value for the ontinuation codes c rpar(1) = 0 means that epsin is the default value rpar(1) = 0.5d0 c for this problem ipar is a dummy variable ipar(1) = 0 return end c 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 = .true. 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----------------------------------------------------- c Initialiser for parameter common block c used only for the R interface c c SUBROUTINE initrparms(bvpparms) EXTERNAL bvpparms c c For Singularly Perturbed BVP (SPBVP) eps should be the c first elements of parms and the first element of rpar c DOUBLE PRECISION parms(1) COMMON / pars / parms CALL bvpparms(1, parms) RETURN END c----------------------------------------------------------------------- subroutine init(neqn,aleft,aright,ms,yval0, nmsh, givmsh, + givey,xguess,yguess,dmguess,nudim) integer, intent(in) :: neqn,nudim,ms(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 = .false. givey = .false. yval0(1:neqn) = 0.0d0 nmsh = 5*6+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 ij = 1 DO i=1,neqn do j=1,ms(i) yguess(ij,1:nmsh) = 0.0d0 ij = ij+1 end do end do DO i=1,neqn dmguess(i,1:nmsh) = 0.0d0 end do return end c----------------------------------------------------------------------- subroutine feval(neqn,x,y,f,rpar,ipar) integer, intent(in) :: neqn,ipar(*) double precision, intent(in):: x,y(1:neqn) double precision, intent(inout):: rpar(*) double precision, intent(out):: f(1:neqn) double precision :: EPS common /pars/ EPS F(1)=Y(2) F(2)=((1.0D0+EPS)*Y(1)-Y(2))/EPS return end c----------------------------------------------------------------------- subroutine jeval(neqn,x,y,dfdy,rpar,ipar) integer, intent(in) :: neqn,ipar(*) double precision, intent(in):: x,y(1:neqn),rpar(*) double precision, intent(out):: dfdy(1:neqn,1:neqn) double precision :: EPS common /pars/ EPS dfdy(1,1) = 0d0 dfdy(1,2) = 1d0 dfdy(2,1) = (1d0+EPS)/eps dfdy(2,2) = -1/eps 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 double precision :: EPS common /pars/ EPS if (i .eq. 1) bc=Y(1)-1.0D0-DEXP(-2.0D0) if (i .eq. 2) bc=Y(1)-1.0D0-DEXP(-2.0D0*(1.0D0+EPS)/EPS) return end subroutine dbceval(i, ncomp,y, dbc, rpar,ipar) integer, intent(in) :: i,ncomp double precision, INTENT(IN) :: y(ncomp) double precision, INTENT(OUT) :: dbc(ncomp) double precision :: EPS common /pars/ EPS dbc(1)=1.D0 dbc(2)=0.D0 return end c c c---------------------------------------------------------------------- c Numerical Jacobian c c SUBROUTINE jeval(R, X0, Y0, JF0, RPAR,IPAR) c SUBROUTINE fnumjac(R, X0, Y0, JF0, RPAR,IPAR) c IMPLICIT NONE c INTEGER, INTENT(IN) :: R, IPAR(*) c DOUBLE PRECISION, INTENT(IN) :: X0, Y0(R), RPAR(*) c DOUBLE PRECISION, INTENT(OUT) :: JF0(1:R,1:R) c EXTERNAL feval c DOUBLE PRECISION :: DN(R), Y0N(R),YSAFE, DELT, X,FP(R),EPSPREC c INTEGER :: I c EPSPREC= epsilon(1.0d0) c CALL feval(R,X0,Y0,FP,RPAR,IPAR) c Y0N(1:R) = Y0(1:R) c DO I=1,R c YSAFE=Y0(I) c DELT=SQRT(1.D-1*EPSPREC*MAX(1.D-5,ABS(YSAFE))) c Y0N(I)=YSAFE+DELT c X = X0 c CALL feval(R,X,Y0N,DN,RPAR,IPAR) c JF0(1:R,I)=(DN(1:R)-FP(1:R))/DELT c Y0N(I)=YSAFE c END DO c RETURN c END c c---------------------------------------------------------------------- c computation of the numerical jacobian for the boundary conditions c c SUBROUTINE dbceval(I,R,y,dbc,RPAR,IPAR) c SUBROUTINE BCNUMJAC(I,R,y,dbc,RPAR,IPAR) c IMPLICIT NONE c INTEGER, INTENT(IN) :: I, R c DOUBLE PRECISION, INTENT(OUT)::dbc(R) c DOUBLE PRECISION, INTENT(IN ):: y(R) c INTEGER, INTENT(IN):: ipar(*) c DOUBLE PRECISION, INTENT(IN) :: rpar(*) c DOUBLE PRECISION :: DN, Y0N(R),YSAFE, DELT, BC, EPSPREC c INTEGER :: J c EXTERNAL bceval c EPSPREC= epsilon(1.0d0) c Y0N(1:R) = Y(1:R) c call bceval(i,R,y,bc,rpar,ipar) c DO J=1,R c YSAFE=Y(J) c DELT=SQRT(1.D-1*EPSPREC*MAX(1.D-5,ABS(YSAFE))) c Y0N(J)=YSAFE+DELT c call bceval(i,R,y0N,DN,rpar,ipar) c DBC(J)=(DN-BC)/DELT c Y0N(J)=YSAFE c END DO c RETURN c END c----------------------------------------------------------------------- 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 common /pars/ lambda DO i=1,nmsh y(1,i) = exp(x(i)-1)+exp(-(1+lambda)*(1+x(i))/lambda) y(2,i) = (exp(X(i) - 1))-((lambda + 1)/lambda)*exp(-(X(i) + 1)* + (lambda + 1)/lambda) END DO return end