PROGRAM CTS
! This program solves the problems of a Test Set formulated by J. Cash. It can
! be downloaded from http://www.ma.ic.ac.uk/~jcash/BVP_software/readme.php
! There are 32 problems. Problems 1-30 are all scalar second order ODEs that
! involve a parameter EPS. They are solved as a pair of first order equations.
! Problems 31-32 are both systems of four first order ODEs. The ACDC code has
! automatic continuation. It asks for starting and ending values of the known
! parameter, so Cash specifies values of EPS for each problem. Here we use
! BVP_SOLVER to solve these problems by continuation in EPS. Continuation was
! actually useful only for problems #24,30,32, so for the other problems we solved
! only for the final value of EPS. We use default settings and in particular, use
! finite difference partial derivatives. We use a nominal error tolerance of 1D-3.
! It proved convenient to define the test set in two files. Specifically, problems
! 1-30 are defined in CTS130 and problems 31-32 in CTS3132. This program can USE
! only one of these files at a time. At run time CTS asks which problem is to be
! solved. After solving the problem, the solution is written to a data file. An
! M-file is provided that will import this data into Matlab and plot the solution.
! Import problem dependent functions and parameters. Comment out one of the
! following lines.
! USE CTS130
USE CTS3132
! Import BVP_SOLVER functions, variables, and types.
USE BVP_M
! Declare a structure for the numerical solution and associated
! information.
TYPE(BVP_SOL) :: SOL
! Working variables
INTEGER :: MXNSUB,I,OK_NPTS
DOUBLE PRECISION :: A,B,EPS0,EPSF,OK_EPS,FACTOR,GUESS(NODE)
! The arrays XGUESS,YGUESS are given their maximum size to
! simplify the program. Storage could be managed more efficiently
! by using ALLOCATABLE arrays.
DOUBLE PRECISION :: XGUESS(3000),YGUESS(NODE,3000)
DOUBLE PRECISION :: XPLOT,YPLOT(4)
! Obtain TP, the number of the test problem. It is a globally
! defined parameter.
PRINT *,'Which test problem (1-32)?'
READ *,TP
CALL PROBLEM_DEFN(A,B,EPS0,EPSF)
! We start with a value of the global parameter EPS for which the
! BVP is easy to solve. Try a default mesh of 10 equally spaced
! points and a constant guess of 1 for all solution components.
! Use a nominal accuracy tolerance of 1D-3.
EPS = EPS0
GUESS = 1D0
SOL = BVP_INIT(NODE,LEFTBC,(/A,B/),GUESS)
SOL = BVP_SOLVER(SOL,FSUB,BCSUB,TOL=1D-3)
OK_EPS = EPS
! Extract needed info from SOL fields.
OK_NPTS = SOL%NPTS
XGUESS(1:OK_NPTS) = SOL%X
YGUESS(:,1:OK_NPTS) = SOL%Y
! If continuation is going to be necessary, deallocate array fields of SOL.
IF (EPS > EPSF) THEN
CALL BVP_TERMINATE(SOL)
END IF
FACTOR = 0.1D0
DO WHILE (EPS > EPSF)
DO
EPS = FACTOR*OK_EPS
IF (EPS <= 1.1D0*EPSF) EPS = EPSF
MXNSUB = MIN(2*OK_NPTS-2,3000)
IF (EPS <= EPSF) MXNSUB = 3000
SOL = BVP_INIT(NODE,LEFTBC,XGUESS(1:OK_NPTS),YGUESS(:,1:OK_NPTS),&
MAX_NUM_SUBINTERVALS=MXNSUB)
SOL = BVP_SOLVER(SOL,FSUB,BCSUB,STOP_ON_FAIL=.FALSE.,TOL=1D-3)
IF (SOL%INFO == -1) THEN
! Computation failed. Restore the solution for OK_EPS and try
! again with a smaller change in EPS. The effect of the second
! argument in MIN is to try an EPS that is an average of the EPS
! that failed and OK_EPS, which succeeded.
FACTOR = MIN(1.2D0*FACTOR,0.5D0*(1D0 + FACTOR))
IF (MXNSUB == 3000) THEN
PRINT *,'Continuation failed.'
STOP
END IF
ELSE
! Computation succeeded. If not done, form a new guess to
! be used with a smaller value of EPS.
IF (EPS <= EPSF) EXIT
OK_EPS = EPS
! Extract needed info from SOL fields and deallocate array
! fields of SOL.
OK_NPTS = SOL%NPTS
XGUESS(1:OK_NPTS) = SOL%X
YGUESS(:,1:OK_NPTS) = SOL%Y
CALL BVP_TERMINATE(SOL)
END IF
END DO
END DO
! Write the solution to an output file. For export to Matlab
! write TP, the number of subintervals, and EPS first.
OPEN(UNIT=7,FILE="CTS.dat")
WRITE(UNIT=7,FMT="(2D25.15)") REAL(TP),REAL(SOL%NPTS+1)
WRITE(UNIT=7,FMT="(2D25.15)") EPS,EPS
DO I = 0,200
XPLOT = A + (I/200D0)*(B - A)
CALL BVP_EVAL(SOL,XPLOT,YPLOT(1:NODE))
! The component plotted depends on the problem.
IF (TP == 31) THEN
WRITE(UNIT=7,FMT="(2D25.15)") XPLOT,YPLOT(3)
ELSEIF (TP == 7 .OR. TP == 32) THEN
WRITE(UNIT=7,FMT="(2D25.15)") XPLOT,YPLOT(2)
ELSE
WRITE(UNIT=7,FMT="(2D25.15)") XPLOT,YPLOT(1)
END IF
END DO
PRINT *,' '
PRINT *,'The solution evaluated at 201 equally spaced points'
PRINT *,'in the interval can be imported into Matlab with'
PRINT *,' '
PRINT *,' >> [x,y] = CTS;'
PRINT *,' '
PRINT *,'For this the output file "CTS.dat" must be'
PRINT *,'in the same directory as the CTS.m function.'
PRINT *,' '
PRINT *,'After importing the data, CTS plots the solution.'
PRINT *,' '
PRINT *,' '
CALL BVP_TERMINATE(SOL)
END PROGRAM CTS