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