!************************************************************************ !* Solve a boundary value problem for a first order DE system of the * !* form: * !* * !* Y1' = F1(X,Y1,Y2,...,YN) * !* Y2' = F2(X,Y1,Y2,...,YN) * !* ... * !* YN' = FN(X,Y1,Y2,...,YN) * !* * !* via the shooting method by determining an approximation for the * !* initial value Y(A). * !* -------------------------------------------------------------------- * !* SAMPLE RUNS; * !* * !* TEST EXAMPLE (WITH METHOD #1): * !* Method: Runge-Kutta embedding formula of 4/5th order. * !* ============= * !* ANALYZED SET OF DIFFERENTIAL EQUATIONS: * !* --------------------------------------- * !* Y'(1) = Y(2) * !* Y'(2) = -Y(1)**3 * !* WITH BOUNDARY CONDITIONS : VALUE OF Y(1) * !* AT POINT A = 0.000000000000000E+000: 0.000000000000000E+000 * !* AT POINT B = 1.000000000000000 : 0.000000000000000E+000 * !* * !* AT LOCATION A = 0.000000000000000E+000 * !* THE FOLLOWING VALUES ARE PROVIDED: * !* * !* Y(1) = 0.000000E+00 * !* Y(2) = 0.120000E+02 * !* * !* REQUIRED PARAMETER: * !* ------------------: * !* -PRECISION FOR THE IVP = 1.000000000000000E-007 * !* -PRECISION FOR THE BVP = 1.000000000000000E-008 * !* -STEP WIDTH FOR THE IVP= 1.000000000000000E-002 * !* -MAX. NUMBER OF F-EVALUATIONS FOR SOLUTION OF THE IVP'S = 20000 * !* -MAX. NUMBER OF NEWTON-ITERATION STEPS = 1000 * !* * !* SOLUTION OF THE PROBLEM: * !* ------------------------ * !* THE START VALUE OF A SOLUTION OF THE BVP IS APPROXIMATED * !* IN POINT A = 0.000000000000000E+000 AS FOLLOWS: * !* * !* Y(1) = 0.0000000000E+00 * !* Y(2) = 0.9722981024E+01 * !* * !* THIS REQUIRES 4 NEWTON-ITERATIONS! * !* * !* DETERMINATION COMPLETED AS PLANNED! * !* * !* * !* TEST EXAMPLE (WITH METHOD #2): * !* Method: Predictor-corrector method of order 4 by Adams-Bashforth- * !* Moulton. * !* ============= * !* ANALYZED SET OF DIFFERENTIAL EQUATIONS: * !* --------------------------------------- * !* Y'(1) = Y(2) * !* Y'(2) = -Y(1)**3 * !* WITH BOUNDARY CONDITIONS : VALUE OF Y(1) * !* AT POINT A = 0.000000000000000E+000: 0.000000000000000E+000 * !* AT POINT B = 1.000000000000000 : 0.000000000000000E+000 * !* * !* AT LOCATION A = 0.000000000000000E+000 * !* THE FOLLOWING VALUES ARE PROVIDED: * !* * !* Y(1) = 0.000000E+00 * !* Y(2) = 0.120000E+02 * !* * !* REQUIRED PARAMETER: * !* ------------------: * !* -PRECISION FOR THE IVP = 1.000000000000000E-007 * !* -PRECISION FOR THE BVP = 1.000000000000000E-008 * !* -STEP WIDTH FOR THE IVP= 1.000000000000000E-002 * !* -MAX. NUMBER OF F-EVALUATIONS FOR SOLUTION OF THE IVP'S = 20000 * !* -MAX. NUMBER OF NEWTON-ITERATION STEPS = 1000 * !* * !* SOLUTION OF THE PROBLEM: * !* ------------------------ * !* THE START VALUE OF A SOLUTION OF THE BVP IS APPROXIMATED * !* IN POINT A = 0.000000000000000E+000 AS FOLLOWS: * !* * !* Y(1) = 0.0000000000E+00 * !* Y(2) = 0.9722980952E+01 * !* * !* THIS REQUIRES 4 NEWTON-ITERATIONS! * !* * !* DETERMINATION COMPLETED AS PLANNED! * !* * !* * !* TEST EXAMPLE (WITH METHOD #3): * !* Method: Extrapolation method of Bulirsch-Stoer. * !* ============= * !* ANALYZED SET OF DIFFERENTIAL EQUATIONS: * !* --------------------------------------- * !* Y'(1) = Y(2) * !* Y'(2) = -Y(1)**3 * !* WITH BOUNDARY CONDITIONS : VALUE OF Y(1) * !* WITH BOUNDARY CONDITIONS : VALUE OF Y(1) * !* AT POINT A = 0.000000000000000E+000: 0.000000000000000E+000 * !* AT POINT B = 1.000000000000000 : 0.000000000000000E+000 * !* * !* AT LOCATION A = 0.000000000000000E+000 * !* THE FOLLOWING VALUES ARE PROVIDED: * !* * !* Y(1) = 0.000000E+00 * !* Y(2) = 0.120000E+02 * !* * !* REQUIRED PARAMETER: * !* ------------------: * !* -PRECISION FOR THE IVP = 1.000000000000000E-007 * !* -PRECISION FOR THE BVP = 1.000000000000000E-008 * !* -STEP WIDTH FOR THE IVP= 1.000000000000000E-002 * !* -MAX. NUMBER OF F-EVALUATIONS FOR SOLUTION OF THE IVP'S = 20000 * !* -MAX. NUMBER OF NEWTON-ITERATION STEPS = 1000 * !* * !* SOLUTION OF THE PROBLEM: * !* ------------------------ * !* THE START VALUE OF A SOLUTION OF THE BVP IS APPROXIMATED * !* IN POINT A = 0.000000000000000E+000 AS FOLLOWS: * !* * !* Y(1) = 0.0000000000E+00 * !* Y(2) = 0.9722981026E+01 * !* * !* THIS REQUIRES 4 NEWTON-ITERATIONS! * !* * !* DETERMINATION COMPLETED AS PLANNED! * !* * !* -------------------------------------------------------------------- * !* REF.: "Numerical Algorithms with C, By Gisela Engeln-Muellges * !* and Frank Uhlig, Springer-Verlag, 1996" [BIBLI 11]. * !* * !* F90 Release By J-P Moreau, Paris. * !* (www.jpmoreau.fr) * !************************************************************************ ! Initialize data; test examples can be exchanged by adjusting ! the subroutine dgl below, randbed in module rwp and ! the variable N (and NMAX=N-1) in main program. Program Test_Rwp integer, parameter :: NMAX=1 character*65 TMethod(1:3) character*12 ystri(0:NMAX) real*8 yanf(0:NMAX), yanf0(0:NMAX) integer N real*8 a, b, h, wa, wb integer ifmax, itmax real*8 epsawp, epsrb integer fehler, iter, method, j, rwp call Init_TMethod(TMethod) N=2 !Number of equations yanf0(0) = 0.d0 yanf0(1) = 12.d0 ystri(0) = 'Y(2)' !alpha-numeric form} ystri(1) = '-Y(1)**3' !of the DE} a = 0.d0 b = 1.d0 h = 0.01d0 wa = 0.d0 wb = 0.d0 ifmax = 20000; itmax = 1000; epsawp = 1d-7 epsrb = 1d-8 !output of the test examples for three varying methods pause 'Adjust size of screen...' do method = 1, 3 call Copy_vector(yanf, yanf0, N) !yanf = yanf0 write(*,10) method write(*,*) '* Method: ', TMethod(method) write(*,*) '=============' write(*,*) 'ANALYZED SET OF DIFFERENTIAL EQUATIONS:' write(*,*) '-------------------------------------- ' do j = 0, N-1 write(*,20) j+1, ystri(j) end do write(*,*) 'WITH BOUNDARY CONDITIONS : VALUE OF Y(1)' write(*,*) '* AT POINT A = ',a,': ',wa write(*,*) '* AT POINT B = ',b,': ',wb write(*,*) '* AT LOCATION A = ',a,', THE FOLLOWING VALUES ARE PROVIDED:' do j = 0, N-1 write(*,30) j+1, yanf(j) end do write(*,*) 'REQUIRED PARAMETER:' write(*,*) '------------------ ' write(*,*) '* -PRECISION FOR THE IVP = ', epsawp write(*,*) '* -PRECISION FOR THE BVP = ', epsrb write(*,*) '* -STEP WIDTH FOR THE IVP= ', h write(*,*) '* -MAX. NUMBER OF F-EVALUATIONS FOR SOLUTION OF THE IVP''S = ', ifmax write(*,*) '* -MAX. NUMBER OF NEWTON-ITERATION STEPS = ', itmax fehler = rwp(a, b, h, yanf, N, method, epsawp, epsrb, ifmax, itmax, iter) !print results if (fehler.ne.0) then write(*,*) '*' Select Case(fehler) case (1) write(*,*) 'ERRORBOUND(S) TOO SMALL' case (2) write(*,*) 'ERROR: B <= A' case (3) write(*,*) 'ERROR: STEP SIZE H <= 0' case (4) write(*,*) 'ERROR: NUMBER N OF DEQ''S NOT CORRECT' case (5) write(*,*) 'ERROR: WRONG IVP-PROGRAM CHOICE' case (6) write(*,*) 'IFMAX TO SMALL FOR SOLUTION OF THE IVP-PROGRAM' case (7) write(*,*) 'AFTER ITMAX STEPS PRECISION WAS NOT REACHED!' case (8) write(*,*) 'ERROR: JACOBI-MATRIX IS SINGULAR!' case (9) write(*,*) 'LACK of MEMORY!' End Select else !no error ? !Output of solution write(*,*) '*' write(*,*) 'SOLUTION OF THE PROBLEM:' write(*,*) '----------------------- ' write(*,*) 'THE START VALUE OF A SOLUTION OF THE BVP IS APPROXIMATED' write(*,*) '* IN POINT A = ', a,' AS FOLLOWS:' do j = 0, N-1 write(*,40) j+1, yanf(j) end do write(*,*) '*' write(*,*) '* THIS REQUIRES ',iter,' NEWTON-ITERATIONS!' write(*,*) 'DETERMINATION COMPLETED AS PLANNED! ' end if pause end do stop 10 format(' * TEST EXAMPLE (WITH METHOD #',I1,')') 20 format(' * Y''(',I1,') = ', A12) 30 format(' * Y''(',I1,') = ', E12.6) 40 format(' * Y(',I1,') = ', E16.10) END !Diff. system to integrate Subroutine dgl(num, n, x, y, f) real*8 y(0:n-1), f(0:n-1) x = x f(0) = y(1) f(1) = -y(0)**3 End; !Init table of method names Subroutine Init_TMethod(TMethod) character*65 TMethod(1:3) TMethod(1)='Runge-Kutta embedding formula of 4/5th order.' TMethod(2)='Predictor-corrector method of order 4 by Adams-Bashforth-Moulton.' TMethod(3)='Extrapolation method of Bulirsch-Stoer.' return End ! --------------------------- END m_rwp.f90 ----------------------------