MODULE DEFINE_FCN ! Define global variables for the number of differential ! equations, NEQNS, and the number of boundary conditions ! at the left end of the interval, LEFTBC. The number of ! boundary conditions at the right end is NODE-LEFTBC. INTEGER, PARAMETER :: NEQNS=2,LEFTBC=1,RIGHTBC=NEQNS-LEFTBC CONTAINS SUBROUTINE FSUB(X,Y,F) DOUBLE PRECISION :: X,Y(NEQNS),F(NEQNS) F = (/ Y(2), -(X*Y(1) - 1D0)*Y(1) /) END SUBROUTINE FSUB SUBROUTINE DF(X,Y,PD) DOUBLE PRECISION :: X,Y(NEQNS),PD(NEQNS,NEQNS) PD = 0D0 PD(1,2) = 1D0 PD(2,1) = 1D0 - 2D0*X*Y(1) END SUBROUTINE DF SUBROUTINE BCSUB(YA,YB,BCA,BCB) DOUBLE PRECISION :: YA(NEQNS),YB(NEQNS),BCA(LEFTBC),BCB(RIGHTBC) BCA(1) = YA(2) BCB(1) = YB(1) + YB(2) END SUBROUTINE BCSUB SUBROUTINE GUESS_Y(X,Y) DOUBLE PRECISION :: X,Y(NEQNS) IF (X <= 1.5D0) THEN Y(1) = 2D0 Y(2) = 0D0 ELSE Y(1) = 2D0*EXP(1.5D0 - X) Y(2) = - Y(1) END IF END SUBROUTINE GUESS_Y END MODULE DEFINE_FCN !************************************************************ PROGRAM ESI ! U.M. Ascher and R.D. Russell, Reformulation of boundary ! value problems into `standard' form, SIAM Review 23 (1981) ! 238-254 use a BVP from electromagnetic self-interaction ! theory to discuss problems with singular coefficients set ! on an infinite interval. After some preparation they solve ! ! u'' + 4*u'/t + (t*u - 1)*u = 0 ! ! with u'(0) = 0 and u(L) + u'(L) = 0. The problem is set on ! an infinite interval, so some experimentation is necessary ! to verify that a sufficiently large L has been specified. ! They present results for u(0) when L = 5,8,10,20. They ! use the initial guess u(t) = 2 for 0 <= t <= 1.5 and ! u(t) = 2*exp(1.5 - t) for t > 1.5. We use this guess for ! the first L and thereafter use the solution for one L as ! the guess for the next, extending to the right with the ! value from their guess (which has the right asymptotic ! behavior). USE DEFINE_FCN USE BVP_M TYPE(BVP_SOL) :: SOL DOUBLE PRECISION :: S(NEQNS,NEQNS),L(4),AR(4),YA(NEQNS),YBNEW(NEQNS) INTEGER :: I ! Define the singular term. S = 0D0 S(2,2) = -4D0 L = (/ 5D0, 8D0, 10D0, 20D0 /) ! Results obtained by Ascher/Russell. AR = (/ 2.03464, 2.11986, 2.11997, 2.11997 /) PRINT *,' L BVP_SOLVER ASCHER/RUSSELL' PRINT *,' ' DO I = 1,4 IF (I == 1) THEN SOL = BVP_INIT(NEQNS,LEFTBC,(/ 0D0,L(I) /),GUESS_Y) ELSE ! Use current approximate solution info to extend initial guess. ! Extend interval to the right. Left end is fixed at 0. CALL BVP_EVAL(SOL,0D0,YA) ! Get solution at 0.0D0. CALL GUESS_Y(L(I),YBNEW) ! Get solution at L(I). SOL = BVP_EXTEND(SOL,0D0,L(I),YA,YBNEW) END IF ! The code can solve this problem with finite difference Jacobians, ! but it is easy to provide the analytical Jacobian and doing so ! helps confirm that the code handles this option properly when ! the BVP is singular. SOL = BVP_SOLVER(SOL,FSUB,BCSUB,SINGULARTERM=S,DFDY=DF) CALL BVP_EVAL(SOL,0.0D0,YA) WRITE(*,FMT="(F5.0,F14.5,F17.5)") L(I),YA(1),AR(I) END DO CALL BVP_TERMINATE(SOL) END PROGRAM ESI