subroutine bisection ( fatol, step_max, nprob, xatol, xa, xb, fxa, fxb ) !*****************************************************************************80 ! !! BISECTION carries out the bisection method to seek a root of F(X) = 0. ! ! Modified: ! ! 17 January 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) FATOL, an absolute error tolerance for ! the function value of the root. If an approximate root X satisfies ! ABS ( F ( X ) ) <= FATOL, then X will be accepted as the ! root and the iteration will be terminated. ! ! Input, integer STEP_MAX, the maximum number of steps allowed ! for an iteration. ! ! Input, integer NPROB, the index of the function whose root is ! to be sought. ! ! Input, real ( kind = 8 ) XATOL, an absolute error tolerance for the root. ! ! Input/output, real ( kind = 8 ) XA, XB, two points at which the ! function differs in sign. On output, these values have been adjusted ! to a smaller interval. ! ! Input/output, real ( kind = 8 ) FXA, FXB, the value of the function ! at XA and XB. ! implicit none real ( kind = 8 ) fatol real ( kind = 8 ) fxa real ( kind = 8 ) fxb real ( kind = 8 ) fxc integer nprob integer step_max integer step_num real ( kind = 8 ) xa real ( kind = 8 ) xatol real ( kind = 8 ) xb real ( kind = 8 ) xc write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Bisection' write ( *, '(a)' ) ' ' write ( *, '(a)' ) & ' Step XA XB F(XA) F(XB)' write ( *, '(a)' ) ' ' ! ! Make A the root with negative F, B the root with positive F. ! if ( 0.0D+00 < fxa ) then call r8_swap ( xa, xb ) call r8_swap ( fxa, fxb ) end if step_num = 0 ! ! Loop ! do write ( *, '(2x,i4,2x,2g16.8,2g14.6)' ) step_num, xa, xb, fxa, fxb step_num = step_num + 1 if ( step_max < step_num ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'BISECTION:' write ( *, '(a)' ) ' Maximum number of steps taken without convergence.' exit end if if ( abs ( xa - xb ) < xatol ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'BISECTION:' write ( *, '(a)' ) ' Interval small enough for convergence.' exit end if if ( abs ( fxa ) <= fatol .or. abs ( fxb ) <= fatol ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'BISECTION:' write ( *, '(a)' ) ' Function small enough for convergence.' exit end if ! ! Compute the next iterate. ! xc = 0.5D+00 * ( xa + xb ) call p00_fx ( nprob, xc, fxc ) ! ! Replace one of the old points. ! if ( fxc < 0.0D+00 ) then xa = xc fxa = fxc else xb = xc fxb = fxc end if end do return end subroutine brent ( fatol, step_max, nprob, xatol, xrtol, xa, xb, fxa, fxb ) !*****************************************************************************80 ! !! BRENT implements the Brent bisection-based zero finder. ! ! Modified: ! ! 17 January 2001 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Richard Brent, ! Algorithms for Minimization without Derivatives, ! Prentice Hall, 1973. ! ! Parameters: ! ! Input, real ( kind = 8 ) FATOL, an absolute error tolerance for the ! function value of the root. If an approximate root X satisfies ! ABS ( F ( X ) ) <= FATOL, then X will be accepted as the ! root and the iteration will be terminated. ! ! Input, integer STEP_MAX, the maximum number of steps allowed ! for an iteration. ! ! Input, integer NPROB, the index of the function whose root is ! to be sought. ! ! Input, real ( kind = 8 ) XATOL, XRTOL, absolute and relative error ! tolerances for the root. ! ! Input/output, real ( kind = 8 ) XA, XB, two points at which the ! function differs in sign. On output, these values have been adjusted ! to a smaller interval. ! ! Input/output, real ( kind = 8 ) FXA, FXB, the value of the function ! at XA and XB. ! implicit none real ( kind = 8 ) d real ( kind = 8 ) e real ( kind = 8 ) fatol real ( kind = 8 ) fxa real ( kind = 8 ) fxb real ( kind = 8 ) fxc integer step_max integer nprob integer step_num real ( kind = 8 ) p real ( kind = 8 ) q real ( kind = 8 ) r real ( kind = 8 ) s real ( kind = 8 ) xa real ( kind = 8 ) xb real ( kind = 8 ) xc real ( kind = 8 ) xm real ( kind = 8 ) xatol real ( kind = 8 ) xrtol real ( kind = 8 ) xtol ! ! Initialization. ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Brent''s Method' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Step XA XB F(XA) F(XB)' write ( *, '(a)' ) ' ' step_num = 0 call p00_fx ( nprob, xa, fxa ) call p00_fx ( nprob, xb, fxb ) ! ! Check that f(ax) and f(bx) have different signs. ! if ( ( fxa < 0.0D+00 .and. fxb < 0.0D+00 ) .or. & ( 0.0D+00 < fxa .and. 0.0D+00 < fxb ) ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'BRENT - Fatal error!' write ( *, '(a)' ) ' F(XA) and F(XB) have same sign.' return end if xc = xa fxc = fxa d = xb - xa e = d do write ( *, '(2x,i4,2x,2g16.8,2g14.6)' ) step_num, xb, xc, fxb, fxc step_num = step_num + 1 if ( step_max < step_num ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'BRENT:' write ( *, '(a)' ) ' Maximum number of steps taken.' exit end if if ( abs ( fxc ) < abs ( fxb ) ) then xa = xb xb = xc xc = xa fxa = fxb fxb = fxc fxc = fxa end if xtol = 2.0D+00 * xrtol * abs ( xb ) + 0.5D+00 * xatol ! ! XM is the halfwidth of the current change-of-sign interval. ! xm = 0.5D+00 * ( xc - xb ) if ( abs ( xm ) <= xtol ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'BRENT:' write ( *, '(a)' ) ' Interval small enough for convergence.' exit end if if ( abs ( fxb ) <= fatol ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'BRENT:' write ( *, '(a)' ) ' Function small enough for convergence.' exit end if ! ! See if a bisection is forced. ! if ( abs ( e ) < xtol .or. abs ( fxa ) <= abs ( fxb ) ) then d = xm e = d else s = fxb / fxa ! ! Linear interpolation. ! if ( xa == xc ) then p = 2.0D+00 * xm * s q = 1.0D+00 - s ! ! Inverse quadratic interpolation. ! else q = fxa / fxc r = fxb / fxc p = s * ( 2.0D+00 * xm * q * ( q - r ) - ( xb - xa ) * ( r - 1.0D+00 ) ) q = ( q - 1.0D+00 ) * ( r - 1.0D+00 ) * ( s - 1.0D+00 ) end if if ( 0.0D+00 < p ) then q = - q else p = - p end if s = e e = d if ( 3.0D+00 * xm * q - abs ( xtol * q ) <= 2.0D+00 * p .or. & abs ( 0.5D+00 * s * q ) <= p ) then d = xm e = d else d = p / q end if end if ! ! Save in XA, FXA the previous values of XB, FXB. ! xa = xb fxa = fxb ! ! Compute the new value of XB, and evaluate the function there. ! if ( xtol < abs ( d ) ) then xb = xb + d else if ( 0.0D+00 < xm ) then xb = xb + xtol else if ( xm <= 0.0D+00 ) then xb = xb - xtol end if call p00_fx ( nprob, xb, fxb ) ! ! If the new FXB has the same sign as FXC, then replace XC by XA. ! if ( ( 0.0D+00 < fxb .and. 0.0D+00 < fxc ) .or. & ( fxb < 0.0D+00 .and. fxc < 0.0D+00 ) ) then xc = xa fxc = fxa d = xb - xa e = d end if end do return end subroutine muller ( fatol, step_max, nprob, xatol, xrtol, xa, xb, xc, fxa, & fxb, fxc ) !*****************************************************************************80 ! !! MULLER carries out Muller's method for seeking a real root of a nonlinear function. ! ! Modified: ! ! 17 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) FATOL, an absolute error tolerance for the ! function value of the root. If an approximate root X satisfies ! ABS ( F ( X ) ) <= FATOL, then X will be accepted as the ! root and the iteration will be terminated. ! ! Input, integer STEP_MAX, the maximum number of steps allowed ! for an iteration. ! ! Input, integer NPROB, the index of the function whose root is ! to be sought. ! ! Input, real ( kind = 8 ) XATOL, XRTOL, absolute and relative error ! tolerances for the root. ! ! Input/output, real ( kind = 8 ) XA, XB, XC, three points. ! ! Input/output, real ( kind = 8 ) FXA, FXB, FXC, the value of the ! function at XA, XB, and XC. ! implicit none real ( kind = 8 ) a real ( kind = 8 ) b real ( kind = 8 ) c real ( kind = 8 ) fatol integer i integer step_max integer nprob real ( kind = 8 ) xa real ( kind = 8 ) xatol real ( kind = 8 ) xb real ( kind = 8 ) xc real ( kind = 8 ) xd real ( kind = 8 ) xrtol real ( kind = 8 ) fxa real ( kind = 8 ) fxb real ( kind = 8 ) fxc real ( kind = 8 ) z1 real ( kind = 8 ) z2 ! ! Initialization. ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Muller''s Method' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Step XA XB XC' write ( *, '(a)' ) ' ' i = 0 write ( *, '(2x,i4,3g14.6)' ) i, xa, xb, xc write ( *, '(2x,4x,3g14.6)' ) fxa, fxb, fxc do i = 1, step_max ! ! Determine the coefficients ! A, B, C ! of the polynomial ! Y(X) = A * (X-X2)**2 + B * (X-X2) + C ! which goes through the data: ! (X1,Y1), (X2,Y2), (X3,Y3). ! a = ( ( fxa - fxc ) * ( xb - xc ) & - ( fxb - fxc ) * ( xa - xc ) ) / & ( ( xa - xc ) * ( xb - xc ) * ( xa - xb ) ) b = ( ( fxb - fxc ) * ( xa - xc )**2 & - ( fxa - fxc ) * ( xb - xc )**2 ) / & ( ( xa - xc ) * ( xb - xc ) * ( xa - xb ) ) c = fxc ! ! Get the real roots of the polynomial, ! unless A = 0, in which case the algorithm is breaking down. ! if ( a /= 0.0D+00 ) then call r8poly2_rroot ( a, b, c, z1, z2 ) else if ( b /= 0.0D+00 ) then z2 = - c / b else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'MULLER:' write ( *, '(a)' ) ' Polynomial fitting has failed.' write ( *, '(a)' ) ' Muller''s algorithm breaks down.' return end if xd = xc + z2 ! ! Set XA, YA, based on which of XA and XB is closer to XD. ! if ( abs ( xd - xb ) < abs ( xd - xa ) ) then call r8_swap ( xa, xb ) call r8_swap ( fxa, fxb ) end if ! ! Set XB, YB, based on which of XB and XC is closer to XD. ! if ( abs ( xd - xc ) < abs ( xd - xb ) ) then call r8_swap ( xb, xc ) call r8_swap ( fxb, fxc ) end if ! ! Set XC, YC. ! xc = xd call p00_fx ( nprob, xc, fxc ) write ( *, '(2x,i4,3g14.6)' ) i, xa, xb, xc write ( *, '(2x,4x,3g14.6)' ) fxa, fxb, fxc ! ! Estimate the relative significance of the most recent correction. ! if ( abs ( z2 ) <= xrtol * abs ( xc ) + xatol ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'MULLER:' write ( *, '(a)' ) ' Stepsize small enough for convergence.' return end if if ( abs ( fxc ) < fatol ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'MULLER:' write ( *, '(a)' ) ' Function small enough for convergence.' return end if end do write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'MULLER:' write ( *, '(a)' ) ' Took maximum number of steps without convergence.' return end subroutine newton ( fatol, step_max, nprob, xatol, xmin, xmax, xa, fxa ) !*****************************************************************************80 ! !! NEWTON carries out Newton's method to seek a root of F(X) = 0. ! ! Modified: ! ! 17 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) FATOL, an absolute error tolerance for the ! function value of the root. If an approximate root X satisfies ! ABS ( F ( X ) ) <= FATOL, then X will be accepted as the ! root and the iteration will be terminated. ! ! Input, integer STEP_MAX, the maximum number of steps allowed ! for an iteration. ! ! Input, integer NPROB, the index of the function whose root is ! to be sought. ! ! Input, real ( kind = 8 ) XATOL, an absolute error tolerance for the root. ! ! Input, real ( kind = 8 ) XMIN, XMAX, the interval in which the root should ! be sought. ! ! Input/output, real ( kind = 8 ) XA. On input, the starting point for ! the iteration. On output, the current approximation to the root. ! ! Input/output, real ( kind = 8 ) FXA, the function value at XA. ! implicit none real ( kind = 8 ) fatol real ( kind = 8 ) fp real ( kind = 8 ) fxa integer step_max integer nprob integer step_num real ( kind = 8 ) step real ( kind = 8 ) xa real ( kind = 8 ) xatol real ( kind = 8 ) xmax real ( kind = 8 ) xmin step = 0.0D+00 write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Newton''s Method' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' X F(X) FP(X)' write ( *, '(a)' ) ' ' step_num = 0 call p00_fx1 ( nprob, xa, fp ) write ( *, '(2x,i4,2x,g16.8,2g14.6)' ) step_num, xa, fxa, fp do step_num = 1, step_max if ( xa < xmin .or. xmax < xa ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Newton:' write ( *, '(a,g14.6)' ) ' The iterate X = ', xa write ( *, '(a)' ) ' has left the region [XMIN,XMAX].' return end if if ( abs ( fxa ) <= fatol )then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Newton:' write ( *, '(a)' ) ' The function norm is small enough for convergence.' return end if if ( 1 < step_num .and. abs ( step ) <= xatol ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Newton:' write ( *, '(a)' ) ' The stepsize is small enough for convergence.' return end if if ( fp == 0.0D+00 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Newton:' write ( *, '(a)' ) ' F''(X)=0, the algorithm fails.' return end if step = fxa / fp xa = xa - step call p00_fx ( nprob, xa, fxa ) call p00_fx1 ( nprob, xa, fp ) write ( *, '(2x,i4,2x,g16.8,2g14.6)' ) step_num, xa, fxa, fp end do write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Newton:' write ( *, '(a)' ) ' Took maximum number of steps without convergence.' return end subroutine p00_fx ( nprob, x, fx ) !*****************************************************************************80 ! !! P00_FX evaluates a function specified by problem number. ! ! Modified: ! ! 19 February 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer NPROB, the number of the problem. ! ! Input, real ( kind = 8 ) X, the point at which F is to be evaluated. ! ! Output, real ( kind = 8 ) FX, the value of the function at X. ! implicit none real ( kind = 8 ) fx integer nprob real ( kind = 8 ) x if ( nprob == 1 ) then call p01_fx ( x, fx ) else if ( nprob == 2 ) then call p02_fx ( x, fx ) else if ( nprob == 3 ) then call p03_fx ( x, fx ) else if ( nprob == 4 ) then call p04_fx ( x, fx ) else if ( nprob == 5 ) then call p05_fx ( x, fx ) else if ( nprob == 6 ) then call p06_fx ( x, fx ) else if ( nprob == 7 ) then call p07_fx ( x, fx ) else if ( nprob == 8 ) then call p08_fx ( x, fx ) else if ( nprob == 9 ) then call p09_fx ( x, fx ) else if ( nprob == 10 ) then call p10_fx ( x, fx ) else if ( nprob == 11 ) then call p11_fx ( x, fx ) else if ( nprob == 12 ) then call p12_fx ( x, fx ) else if ( nprob == 13 ) then call p13_fx ( x, fx ) else if ( nprob == 14 ) then call p14_fx ( x, fx ) else if ( nprob == 15 ) then call p15_fx ( x, fx ) else if ( nprob == 16 ) then call p16_fx ( x, fx ) else if ( nprob == 17 ) then call p17_fx ( x, fx ) else if ( nprob == 18 ) then call p18_fx ( x, fx ) else if ( nprob == 19 ) then call p19_fx ( x, fx ) else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'P00_FX - Fatal error!' write ( *, '(a,i4)' ) ' Illegal problem number = ', nprob stop end if return end subroutine p00_fx1 ( nprob, x, fx1 ) !*****************************************************************************80 ! !! P00_FX1 evaluates the first derivative of a function specified by problem number. ! ! Modified: ! ! 19 February 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer NPROB, the number of the problem. ! ! Input, real ( kind = 8 ) X, the abscissa. ! ! Output, real ( kind = 8 ) FX1, the first derivative of the function at X. ! implicit none ! real ( kind = 8 ) fx1 integer nprob real ( kind = 8 ) x if ( nprob == 1 ) then call p01_fx1 ( x, fx1 ) else if ( nprob == 2 ) then call p02_fx1 ( x, fx1 ) else if ( nprob == 3 ) then call p03_fx1 ( x, fx1 ) else if ( nprob == 4 ) then call p04_fx1 ( x, fx1 ) else if ( nprob == 5 ) then call p05_fx1 ( x, fx1 ) else if ( nprob == 6 ) then call p06_fx1 ( x, fx1 ) else if ( nprob == 7 ) then call p07_fx1 ( x, fx1 ) else if ( nprob == 8 ) then call p08_fx1 ( x, fx1 ) else if ( nprob == 9 ) then call p09_fx1 ( x, fx1 ) else if ( nprob == 10 ) then call p10_fx1 ( x, fx1 ) else if ( nprob == 11 ) then call p11_fx1 ( x, fx1 ) else if ( nprob == 12 ) then call p12_fx1 ( x, fx1 ) else if ( nprob == 13 ) then call p13_fx1 ( x, fx1 ) else if ( nprob == 14 ) then call p14_fx1 ( x, fx1 ) else if ( nprob == 15 ) then call p15_fx1 ( x, fx1 ) else if ( nprob == 16 ) then call p16_fx1 ( x, fx1 ) else if ( nprob == 17 ) then call p17_fx1 ( x, fx1 ) else if ( nprob == 18 ) then call p18_fx1 ( x, fx1 ) else if ( nprob == 19 ) then call p19_fx1 ( x, fx1 ) else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'P00_FX1 - Fatal error!' write ( *, '(a,i4)' ) ' Illegal problem number = ', nprob stop end if return end subroutine p00_fx2 ( nprob, x, fx2 ) !*****************************************************************************80 ! !! P00_FX2 evaluates the second derivative of a function specified by problem number. ! ! Modified: ! ! 19 February 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer NPROB, the number of the problem. ! ! Input, real ( kind = 8 ) X, the abscissa. ! ! Output, real ( kind = 8 ) FX2, the second derivative of the function at X. ! implicit none real ( kind = 8 ) fx2 integer nprob real ( kind = 8 ) x if ( nprob == 1 ) then call p01_fx2 ( x, fx2 ) else if ( nprob == 2 ) then call p02_fx2 ( x, fx2 ) else if ( nprob == 3 ) then call p03_fx2 ( x, fx2 ) else if ( nprob == 4 ) then call p04_fx2 ( x, fx2 ) else if ( nprob == 5 ) then call p05_fx2 ( x, fx2 ) else if ( nprob == 6 ) then call p06_fx2 ( x, fx2 ) else if ( nprob == 7 ) then call p07_fx2 ( x, fx2 ) else if ( nprob == 8 ) then call p08_fx2 ( x, fx2 ) else if ( nprob == 9 ) then call p09_fx2 ( x, fx2 ) else if ( nprob == 10 ) then call p10_fx2 ( x, fx2 ) else if ( nprob == 11 ) then call p11_fx2 ( x, fx2 ) else if ( nprob == 12 ) then call p12_fx2 ( x, fx2 ) else if ( nprob == 13 ) then call p13_fx2 ( x, fx2 ) else if ( nprob == 14 ) then call p14_fx2 ( x, fx2 ) else if ( nprob == 15 ) then call p15_fx2 ( x, fx2 ) else if ( nprob == 16 ) then call p16_fx2 ( x, fx2 ) else if ( nprob == 17 ) then call p17_fx2 ( x, fx2 ) else if ( nprob == 18 ) then call p18_fx2 ( x, fx2 ) else if ( nprob == 19 ) then call p19_fx2 ( x, fx2 ) else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'P00_FX2 - Fatal error!' write ( *, '(a,i4)' ) ' Illegal problem number = ', nprob stop end if return end subroutine p00_interval ( nprob, xmin, xmax ) !*****************************************************************************80 ! !! P00_INTERVAL returns an interval bounding the root for any problem. ! ! Modified: ! ! 19 February 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer NPROB, the number of the problem. ! ! Output, real ( kind = 8 ) XMIN, XMAX, the minimum and maximum values of ! an interval containing the root. ! implicit none integer nprob real ( kind = 8 ) xmax real ( kind = 8 ) xmin if ( nprob == 1 ) then call p01_interval ( xmin, xmax ) else if ( nprob == 2 ) then call p02_interval ( xmin, xmax ) else if ( nprob == 3 ) then call p03_interval ( xmin, xmax ) else if ( nprob == 4 ) then call p04_interval ( xmin, xmax ) else if ( nprob == 5 ) then call p05_interval ( xmin, xmax ) else if ( nprob == 6 ) then call p06_interval ( xmin, xmax ) else if ( nprob == 7 ) then call p07_interval ( xmin, xmax ) else if ( nprob == 8 ) then call p08_interval ( xmin, xmax ) else if ( nprob == 9 ) then call p09_interval ( xmin, xmax ) else if ( nprob == 10 ) then call p10_interval ( xmin, xmax ) else if ( nprob == 11 ) then call p11_interval ( xmin, xmax ) else if ( nprob == 12 ) then call p12_interval ( xmin, xmax ) else if ( nprob == 13 ) then call p13_interval ( xmin, xmax ) else if ( nprob == 14 ) then call p14_interval ( xmin, xmax ) else if ( nprob == 15 ) then call p15_interval ( xmin, xmax ) else if ( nprob == 16 ) then call p16_interval ( xmin, xmax ) else if ( nprob == 17 ) then call p17_interval ( xmin, xmax ) else if ( nprob == 18 ) then call p18_interval ( xmin, xmax ) else if ( nprob == 19 ) then call p19_interval ( xmin, xmax ) else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'P00_INTERVAL - Fatal error!' write ( *, '(a,i4)' ) ' Illegal problem number = ', nprob stop end if return end subroutine p00_nprob ( nprob ) !*****************************************************************************80 ! !! P00_NPROB returns the number of problems available. ! ! Modified: ! ! 19 February 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer NPROB, the number of problems available. ! implicit none integer nprob nprob = 19 return end subroutine p00_root ( nprob, max_root, nroot, x ) !*****************************************************************************80 ! !! P00_ROOT returns a root for any problem. ! ! Modified: ! ! 19 February 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer NPROB, the number of the problem. ! ! Input, integer MAX_ROOT, the maximum number of roots that can be stored. ! ! Output, integer NROOT, the number of roots returned. If NROOT is zero, ! then no exact roots are returned. ! ! Output, real ( kind = 8 ) X(MAX_ROOT), NROOT roots for the given problem. ! implicit none integer max_root integer nprob integer nroot real ( kind = 8 ) x(max_root) if ( nprob == 1 ) then call p01_root ( max_root, nroot, x ) else if ( nprob == 2 ) then call p02_root ( max_root, nroot, x ) else if ( nprob == 3 ) then call p03_root ( max_root, nroot, x ) else if ( nprob == 4 ) then call p04_root ( max_root, nroot, x ) else if ( nprob == 5 ) then call p05_root ( max_root, nroot, x ) else if ( nprob == 6 ) then call p06_root ( max_root, nroot, x ) else if ( nprob == 7 ) then call p07_root ( max_root, nroot, x ) else if ( nprob == 8 ) then call p08_root ( max_root, nroot, x ) else if ( nprob == 9 ) then call p09_root ( max_root, nroot, x ) else if ( nprob == 10 ) then call p10_root ( max_root, nroot, x ) else if ( nprob == 11 ) then call p11_root ( max_root, nroot, x ) else if ( nprob == 12 ) then call p12_root ( max_root, nroot, x ) else if ( nprob == 13 ) then call p13_root ( max_root, nroot, x ) else if ( nprob == 14 ) then call p14_root ( max_root, nroot, x ) else if ( nprob == 15 ) then call p15_root ( max_root, nroot, x ) else if ( nprob == 16 ) then call p16_root ( max_root, nroot, x ) else if ( nprob == 17 ) then call p17_root ( max_root, nroot, x ) else if ( nprob == 18 ) then call p18_root ( max_root, nroot, x ) else if ( nprob == 19 ) then call p19_root ( max_root, nroot, x ) else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'P00_ROOT - Fatal error!' write ( *, '(a,i4)' ) ' Illegal problem number = ', nprob stop end if return end subroutine p00_start ( nprob, max_start, nstart, x ) !*****************************************************************************80 ! !! P00_START returns one or more starting points for any problem. ! ! Modified: ! ! 19 February 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer NPROB, the number of the problem. ! ! Input, integer MAX_START; the amount of space in X. ! ! Output, integer NSTART; the number of starting points, which will ! never be more than MAX_START. ! ! Output, real ( kind = 8 ) X(MAX_START), contains NSTART starting points. ! implicit none integer max_start integer nprob integer nstart real ( kind = 8 ) x(max_start) if ( nprob == 1 ) then call p01_start ( max_start, nstart, x ) else if ( nprob == 2 ) then call p02_start ( max_start, nstart, x ) else if ( nprob == 3 ) then call p03_start ( max_start, nstart, x ) else if ( nprob == 4 ) then call p04_start ( max_start, nstart, x ) else if ( nprob == 5 ) then call p05_start ( max_start, nstart, x ) else if ( nprob == 6 ) then call p06_start ( max_start, nstart, x ) else if ( nprob == 7 ) then call p07_start ( max_start, nstart, x ) else if ( nprob == 8 ) then call p08_start ( max_start, nstart, x ) else if ( nprob == 9 ) then call p09_start ( max_start, nstart, x ) else if ( nprob == 10 ) then call p10_start ( max_start, nstart, x ) else if ( nprob == 11 ) then call p11_start ( max_start, nstart, x ) else if ( nprob == 12 ) then call p12_start ( max_start, nstart, x ) else if ( nprob == 13 ) then call p13_start ( max_start, nstart, x ) else if ( nprob == 14 ) then call p14_start ( max_start, nstart, x ) else if ( nprob == 15 ) then call p15_start ( max_start, nstart, x ) else if ( nprob == 16 ) then call p16_start ( max_start, nstart, x ) else if ( nprob == 17 ) then call p17_start ( max_start, nstart, x ) else if ( nprob == 18 ) then call p18_start ( max_start, nstart, x ) else if ( nprob == 19 ) then call p19_start ( max_start, nstart, x ) else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'P00_START - Fatal error!' write ( *, '(a,i4)' ) ' Illegal problem number = ', nprob stop end if return end subroutine p00_title ( nprob, title ) !*****************************************************************************80 ! !! P00_TITLE returns the title for a given problem. ! ! Modified: ! ! 19 February 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer NPROB, the number of the problem. ! ! Output, character ( len = * ) TITLE, the title of the given problem. ! implicit none integer nprob character ( len = * ) title if ( nprob == 1 ) then call p01_title ( title ) else if ( nprob == 2 ) then call p02_title ( title ) else if ( nprob == 3 ) then call p03_title ( title ) else if ( nprob == 4 ) then call p04_title ( title ) else if ( nprob == 5 ) then call p05_title ( title ) else if ( nprob == 6 ) then call p06_title ( title ) else if ( nprob == 7 ) then call p07_title ( title ) else if ( nprob == 8 ) then call p08_title ( title ) else if ( nprob == 9 ) then call p09_title ( title ) else if ( nprob == 10 ) then call p10_title ( title ) else if ( nprob == 11 ) then call p11_title ( title ) else if ( nprob == 12 ) then call p12_title ( title ) else if ( nprob == 13 ) then call p13_title ( title ) else if ( nprob == 14 ) then call p14_title ( title ) else if ( nprob == 15 ) then call p15_title ( title ) else if ( nprob == 16 ) then call p16_title ( title ) else if ( nprob == 17 ) then call p17_title ( title ) else if ( nprob == 18 ) then call p18_title ( title ) else if ( nprob == 19 ) then call p19_title ( title ) else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'P00_TITLE - Fatal error!' write ( *, '(a,i4)' ) ' Illegal problem number = ', nprob stop end if return end subroutine p01_fx ( x, fx ) !*****************************************************************************80 ! !! P01_FX evaluates sin ( x ) - x / 2. ! ! Modified: ! ! 07 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the point at which F is to be evaluated. ! ! Output, real ( kind = 8 ) FX, the value of the function at X. ! implicit none real ( kind = 8 ) fx real ( kind = 8 ) x fx = sin ( x ) - 0.5D+00 * x return end subroutine p01_fx1 ( x, fx1 ) !*****************************************************************************80 ! !! P01_FX1 evaluates the derivative of the function for problem 1. ! ! Modified: ! ! 08 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the abscissa. ! ! Output, real ( kind = 8 ) FX1, the first derivative of the function at X. ! implicit none real ( kind = 8 ) fx1 real ( kind = 8 ) x fx1 = cos ( x ) - 0.5D+00 return end subroutine p01_fx2 ( x, fx2 ) !*****************************************************************************80 ! !! P01_FX2 evaluates the second derivative of the function for problem 1. ! ! Modified: ! ! 08 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the abscissa. ! ! Output, real ( kind = 8 ) FX2, the second derivative of the function at X. ! implicit none real ( kind = 8 ) fx2 real ( kind = 8 ) x fx2 = - sin ( x ) return end subroutine p01_interval ( xmin, xmax ) !*****************************************************************************80 ! !! P01_INTERVAL returns an interval bounding the root for problem 1. ! ! Modified: ! ! 07 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real ( kind = 8 ) XMIN, XMAX, the minimum and maximum values of ! an interval containing the root. ! implicit none real ( kind = 8 ) xmax real ( kind = 8 ) xmin xmin = -1000.0D+00 xmax = 1000.0D+00 return end subroutine p01_root ( max_root, nroot, x ) !*****************************************************************************80 ! !! P01_ROOT returns a root for problem 1. ! ! Modified: ! ! 07 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAX_ROOT, the maximum number of roots that can be stored. ! ! Output, integer NROOT, the number of roots returned. ! ! Output, real ( kind = 8 ) X(MAX_ROOT), NROOT roots for the given problem. ! implicit none integer max_root integer nroot real ( kind = 8 ) x(max_root) nroot = 0 if ( nroot < max_root ) then nroot = nroot + 1 x(nroot) = - 1.895494342D+00 end if if ( nroot < max_root ) then nroot = nroot + 1 x(nroot) = 0.0D+00 end if if ( nroot < max_root ) then nroot = nroot + 1 x(nroot) = 1.895494342D+00 end if return end subroutine p01_start ( max_start, nstart, x ) !*****************************************************************************80 ! !! P01_START returns one or more starting points for problem 1. ! ! Modified: ! ! 15 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAX_START; the amount of space in X. ! ! Output, integer NSTART; the number of starting points, which will ! never be more than MAX_START. ! ! Output, real ( kind = 8 ) X(MAX_START), contains NSTART starting points. ! implicit none integer max_start integer nstart real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00 real ( kind = 8 ) x(max_start) nstart = 0 if ( nstart < max_start ) then nstart = nstart + 1 x(nstart) = 0.5D+00 * pi end if if ( nstart < max_start ) then nstart = nstart + 1 x(nstart) = pi end if return end subroutine p01_title ( title ) !*****************************************************************************80 ! !! P01_TITLE returns the title of problem 1. ! ! Modified: ! ! 07 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, character ( len = * ) TITLE, the title of the problem. ! implicit none character ( len = * ) title title = 'F(X) = SIN(X) - 0.5 * X' return end subroutine p02_fx ( x, fx ) !*****************************************************************************80 ! !! P02_FX evaluates 2 * x - exp ( - x ). ! ! Modified: ! ! 07 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the point at which F is to be evaluated. ! ! Output, real ( kind = 8 ) FX, the value of the function at X. ! implicit none real ( kind = 8 ) fx real ( kind = 8 ) x fx = 2.0D+00 * x - exp ( - x ) return end subroutine p02_fx1 ( x, fx1 ) !*****************************************************************************80 ! !! P02_FX1 evaluates the derivative of the function for problem 2. ! ! Modified: ! ! 08 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the abscissa. ! ! Output, real ( kind = 8 ) FX1, the first derivative of the function at X. ! implicit none real ( kind = 8 ) fx1 real ( kind = 8 ) x fx1 = 2.0D+00 + exp ( - x ) return end subroutine p02_fx2 ( x, fx2 ) !*****************************************************************************80 ! !! P02_FX2 evaluates the second derivative of the function for problem 2. ! ! Modified: ! ! 08 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the abscissa. ! ! Output, real ( kind = 8 ) FX2, the second derivative of the function at X. ! implicit none real ( kind = 8 ) fx2 real ( kind = 8 ) x fx2 = - exp ( - x ) return end subroutine p02_interval ( xmin, xmax ) !*****************************************************************************80 ! !! P02_INTERVAL returns an interval bounding the root for problem 2. ! ! Modified: ! ! 07 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real ( kind = 8 ) XMIN, XMAX, the minimum and maximum values of ! an interval containing the root. ! implicit none real ( kind = 8 ) xmax real ( kind = 8 ) xmin xmin = -10.0D+00 xmax = 100.0D+00 return end subroutine p02_root ( max_root, nroot, x ) !*****************************************************************************80 ! !! P02_ROOT returns a root for problem 2. ! ! Modified: ! ! 07 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAX_ROOT, the maximum number of roots that can be stored. ! ! Output, integer NROOT, the number of roots returned. ! ! Output, real ( kind = 8 ) X(MAX_ROOT), NROOT roots for the given problem. ! implicit none integer max_root integer nroot real ( kind = 8 ) x(max_root) nroot = 0 if ( nroot < max_root ) then nroot = nroot + 1 x(nroot) = 0.3517337143D+00 end if return end subroutine p02_start ( max_start, nstart, x ) !*****************************************************************************80 ! !! P02_START returns one or more starting points for problem 2. ! ! Modified: ! ! 24 August 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAX_START; the amount of space in X. ! ! Output, integer NSTART; the number of starting points, which will ! never be more than MAX_START. ! ! Output, real ( kind = 8 ) X(MAX_START), contains NSTART starting points. ! implicit none integer max_start integer nstart real ( kind = 8 ) x(max_start) nstart = 0 if ( nstart < max_start ) then nstart = nstart + 1 x(nstart) = 0.0D+00 end if if ( nstart < max_start ) then nstart = nstart + 1 x(nstart) = 1.0D+00 end if if ( nstart < max_start ) then nstart = nstart + 1 x(nstart) = -5.0D+00 end if if ( nstart < max_start ) then nstart = nstart + 1 x(nstart) = 10.0D+00 end if return end subroutine p02_title ( title ) !*****************************************************************************80 ! !! P02_TITLE returns the title of problem 2. ! ! Modified: ! ! 07 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, character ( len = * ) TITLE, the title of the problem. ! implicit none character ( len = * ) title title = 'F(X) = 2 * X - EXP ( - X )' return end subroutine p03_fx ( x, fx ) !*****************************************************************************80 ! !! P03_FX evaluates x * exp ( - x ). ! ! Modified: ! ! 07 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the point at which F is to be evaluated. ! ! Output, real ( kind = 8 ) FX, the value of the function at X. ! implicit none real ( kind = 8 ) fx real ( kind = 8 ) x fx = x * exp ( - x ) return end subroutine p03_fx1 ( x, fx1 ) !*****************************************************************************80 ! !! P03_FX1 evaluates the derivative of the function for problem 3. ! ! Modified: ! ! 08 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the abscissa. ! ! Output, real ( kind = 8 ) FX1, the first derivative of the function at X. ! implicit none real ( kind = 8 ) fx1 real ( kind = 8 ) x fx1 = exp ( - x ) * ( 1.0D+00 - x ) return end subroutine p03_fx2 ( x, fx2 ) !*****************************************************************************80 ! !! P03_FX2 evaluates the second derivative of the function for problem 3. ! ! Modified: ! ! 08 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the abscissa. ! ! Output, real ( kind = 8 ) FX2, the second derivative of the function at X. ! implicit none real ( kind = 8 ) fx2 real ( kind = 8 ) x fx2 = exp ( - x ) * ( x - 2.0D+00 ) return end subroutine p03_interval ( xmin, xmax ) !*****************************************************************************80 ! !! P03_INTERVAL returns an interval bounding the root for problem 3. ! ! Modified: ! ! 07 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real ( kind = 8 ) XMIN, XMAX, the minimum and maximum values of ! an interval containing the root. ! implicit none real ( kind = 8 ) xmax real ( kind = 8 ) xmin xmin = -10.0D+00 xmax = 100.0D+00 return end subroutine p03_root ( max_root, nroot, x ) !*****************************************************************************80 ! !! P03_ROOT returns a root for problem 3. ! ! Modified: ! ! 07 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAX_ROOT, the maximum number of roots that can be stored. ! ! Output, integer NROOT, the number of roots returned. ! ! Output, real ( kind = 8 ) X(MAX_ROOT), NROOT roots for the given problem. ! implicit none integer max_root integer nroot real ( kind = 8 ) x(max_root) nroot = 0 if ( nroot < max_root ) then nroot = nroot + 1 x(nroot) = 0.0D+00 end if return end subroutine p03_start ( max_start, nstart, x ) !*****************************************************************************80 ! !! P03_START returns one or more starting points for problem 3. ! ! Modified: ! ! 15 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAX_START; the amount of space in X. ! ! Output, integer NSTART; the number of starting points, which will ! never be more than MAX_START. ! ! Output, real ( kind = 8 ) X(MAX_START), contains NSTART starting points. ! implicit none integer max_start integer nstart real ( kind = 8 ) x(max_start) nstart = 0 if ( nstart < max_start ) then nstart = nstart + 1 x(nstart) = -1.0D+00 end if if ( nstart < max_start ) then nstart = nstart + 1 x(nstart) = 0.5D+00 end if if ( nstart < max_start ) then nstart = nstart + 1 x(nstart) = 2.0D+00 end if return end subroutine p03_title ( title ) !*****************************************************************************80 ! !! P03_TITLE returns the title of problem 3. ! ! Modified: ! ! 07 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, character ( len = * ) TITLE, the title of the problem. ! implicit none character ( len = * ) title title = 'F(X) = X * EXP ( - X )' return end subroutine p04_fx ( x, fx ) !*****************************************************************************80 ! !! P04_FX evaluates exp ( x ) - 1 / ( 10 * x )^2. ! ! Modified: ! ! 07 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the point at which F is to be evaluated. ! ! Output, real ( kind = 8 ) FX, the value of the function at X. ! implicit none real ( kind = 8 ) fx real ( kind = 8 ) x fx = exp ( x ) - 1.0D+00 / ( 10.0D+00 * x )**2 return end subroutine p04_fx1 ( x, fx1 ) !*****************************************************************************80 ! !! P04_FX1 evaluates the derivative of the function for problem 4. ! ! Modified: ! ! 08 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the abscissa. ! ! Output, real ( kind = 8 ) FX1, the first derivative of the function at X. ! implicit none real ( kind = 8 ) fx1 real ( kind = 8 ) x fx1 = exp ( x ) + 2.0D+00 / ( 100.0D+00 * x**3 ) return end subroutine p04_fx2 ( x, fx2 ) !*****************************************************************************80 ! !! P04_FX2 evaluates the second derivative of the function for problem 4. ! ! Modified: ! ! 08 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the abscissa. ! ! Output, real ( kind = 8 ) FX2, the second derivative of the function at X. ! implicit none real ( kind = 8 ) fx2 real ( kind = 8 ) x fx2 = exp ( x ) - 6.0D+00 / ( 100.0D+00 * x**4 ) return end subroutine p04_interval ( xmin, xmax ) !*****************************************************************************80 ! !! P04_INTERVAL returns an interval bounding the root for problem 4. ! ! Modified: ! ! 07 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real ( kind = 8 ) XMIN, XMAX, the minimum and maximum values of ! an interval containing the root. ! implicit none real ( kind = 8 ) xmax real ( kind = 8 ) xmin xmin = 0.00001D+00 xmax = 20.0D+00 return end subroutine p04_root ( max_root, nroot, x ) !*****************************************************************************80 ! !! P04_ROOT returns a root for problem 4. ! ! Modified: ! ! 07 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAX_ROOT, the maximum number of roots that can be stored. ! ! Output, integer NROOT, the number of roots returned. ! ! Output, real ( kind = 8 ) X(MAX_ROOT), NROOT roots for the given problem. ! implicit none integer max_root integer nroot real ( kind = 8 ) x(max_root) nroot = 0 if ( nroot < max_root ) then nroot = nroot + 1 x(nroot) = 0.095344618D+00 end if return end subroutine p04_start ( max_start, nstart, x ) !*****************************************************************************80 ! !! P04_START returns one or more starting points for problem 4. ! ! Modified: ! ! 15 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAX_START; the amount of space in X. ! ! Output, integer NSTART; the number of starting points, which will ! never be more than MAX_START. ! ! Output, real ( kind = 8 ) X(MAX_START), contains NSTART starting points. ! implicit none integer max_start integer nstart real ( kind = 8 ) x(max_start) nstart = 0 if ( nstart < max_start ) then nstart = nstart + 1 x(nstart) = 0.03D+00 end if if ( nstart < max_start ) then nstart = nstart + 1 x(nstart) = 1.0D+00 end if return end subroutine p04_title ( title ) !*****************************************************************************80 ! !! P04_TITLE returns the title of problem 4. ! ! Modified: ! ! 07 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, character ( len = * ) TITLE, the title of the problem. ! implicit none character ( len = * ) title title = 'F(X) = EXP ( X ) - 1 / ( 10 * X )**2' return end subroutine p05_fx ( x, fx ) !*****************************************************************************80 ! !! P05_FX evaluates ( x + 3 ) * ( x - 1 )^2. ! ! Modified: ! ! 07 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the point at which F is to be evaluated. ! ! Output, real ( kind = 8 ) FX, the value of the function at X. ! implicit none real ( kind = 8 ) fx real ( kind = 8 ) x fx = ( x + 3.0D+00 ) * ( x - 1.0D+00 )**2 return end subroutine p05_fx1 ( x, fx1 ) !*****************************************************************************80 ! !! P05_FX1 evaluates the derivative of the function for problem 5. ! ! Modified: ! ! 08 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the abscissa. ! ! Output, real ( kind = 8 ) FX1, the first derivative of the function at X. ! implicit none real ( kind = 8 ) fx1 real ( kind = 8 ) x fx1 = ( 3.0D+00 * x + 5.0D+00 ) * ( x - 1.0D+00 ) return end subroutine p05_fx2 ( x, fx2 ) !*****************************************************************************80 ! !! P05_FX2 evaluates the second derivative of the function for problem 5. ! ! Modified: ! ! 08 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the abscissa. ! ! Output, real ( kind = 8 ) FX2, the second derivative of the function at X. ! implicit none real ( kind = 8 ) fx2 real ( kind = 8 ) x fx2 = 6.0D+00 * x + 2.0D+00 return end subroutine p05_interval ( xmin, xmax ) !*****************************************************************************80 ! !! P05_INTERVAL returns an interval bounding the root for problem 5. ! ! Modified: ! ! 07 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real ( kind = 8 ) XMIN, XMAX, the minimum and maximum values of ! an interval containing the root. ! implicit none real ( kind = 8 ) xmax real ( kind = 8 ) xmin xmin = -1000.0D+00 xmax = 1000.0D+00 return end subroutine p05_root ( max_root, nroot, x ) !*****************************************************************************80 ! !! P05_ROOT returns a root for problem 5. ! ! Modified: ! ! 07 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAX_ROOT, the maximum number of roots that can be stored. ! ! Output, integer NROOT, the number of roots returned. ! ! Output, real ( kind = 8 ) X(MAX_ROOT), NROOT roots for the given problem. ! implicit none integer max_root integer nroot real ( kind = 8 ) x(max_root) nroot = 0 if ( nroot < max_root ) then nroot = nroot + 1 x(nroot) = - 3.0D+00 end if if ( nroot < max_root ) then nroot = nroot + 1 x(nroot) = 1.0D+00 end if if ( nroot < max_root ) then nroot = nroot + 1 x(nroot) = 1.0D+00 end if return end subroutine p05_start ( max_start, nstart, x ) !*****************************************************************************80 ! !! P05_START returns one or more starting points for problem 5. ! ! Modified: ! ! 15 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAX_START; the amount of space in X. ! ! Output, integer NSTART; the number of starting points, which will ! never be more than MAX_START. ! ! Output, real ( kind = 8 ) X(MAX_START), contains NSTART starting points. ! implicit none integer max_start integer nstart real ( kind = 8 ) x(max_start) nstart = 0 if ( nstart < max_start ) then nstart = nstart + 1 x(nstart) = 2.0D+00 end if if ( nstart < max_start ) then nstart = nstart + 1 x(nstart) = - 5.0D+00 end if return end subroutine p05_title ( title ) !*****************************************************************************80 ! !! P05_TITLE returns the title of problem 5. ! ! Modified: ! ! 07 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, character ( len = * ) TITLE, the title of the problem. ! implicit none character ( len = * ) title title = 'F(X) = ( X + 3 ) * ( X - 1 )**2' return end subroutine p06_fx ( x, fx ) !*****************************************************************************80 ! !! P06_FX evaluates exp ( x ) - 2 - 1 / ( 10 * x )^2 + 2 / ( 100 * x )^3. ! ! Modified: ! ! 07 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the point at which F is to be evaluated. ! ! Output, real ( kind = 8 ) FX, the value of the function at X. ! implicit none real ( kind = 8 ) fx real ( kind = 8 ) x fx = exp ( x ) - 2.0D+00 - 1.0D+00 / ( 10.0D+00 * x )**2 & + 2.0D+00 / ( 100.0D+00 * x )**3 return end subroutine p06_fx1 ( x, fx1 ) !*****************************************************************************80 ! !! P06_FX1 evaluates the derivative of the function for problem 6. ! ! Modified: ! ! 08 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the abscissa. ! ! Output, real ( kind = 8 ) FX1, the first derivative of the function at X. ! implicit none real ( kind = 8 ) fx1 real ( kind = 8 ) x fx1 = exp ( x ) + 2.0D+00 / ( 100.0D+00 * x**3 ) & - 6.0D+00 / ( 1000000.0D+00 * x**4 ) return end subroutine p06_fx2 ( x, fx2 ) !*****************************************************************************80 ! !! P06_FX2 evaluates the second derivative of the function for problem 6. ! ! Modified: ! ! 08 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the abscissa. ! ! Output, real ( kind = 8 ) FX2, the second derivative of the function at X. ! implicit none real ( kind = 8 ) fx2 real ( kind = 8 ) x fx2 = exp ( x ) - 6.0D+00 / ( 100.0D+00 * x**4 ) & + 24.0D+00 / ( 1000000.0D+00 * x**5 ) return end subroutine p06_interval ( xmin, xmax ) !*****************************************************************************80 ! !! P06_INTERVAL returns an interval bounding the root for problem 6. ! ! Modified: ! ! 07 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real ( kind = 8 ) XMIN, XMAX, the minimum and maximum values of ! an interval containing the root. ! implicit none real ( kind = 8 ) xmax real ( kind = 8 ) xmin xmin = 0.00001D+00 xmax = 20.0D+00 return end subroutine p06_root ( max_root, nroot, x ) !*****************************************************************************80 ! !! P06_ROOT returns a root for problem 6. ! ! Modified: ! ! 07 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAX_ROOT, the maximum number of roots that can be stored. ! ! Output, integer NROOT, the number of roots returned. ! ! Output, real ( kind = 8 ) X(MAX_ROOT), NROOT roots for the given problem. ! implicit none integer max_root integer nroot real ( kind = 8 ) x(max_root) nroot = 0 if ( nroot < max_root ) then nroot = nroot + 1 x(nroot) = 0.70320484D+00 end if return end subroutine p06_start ( max_start, nstart, x ) !*****************************************************************************80 ! !! P06_START returns one or more starting points for problem 6. ! ! Modified: ! ! 15 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAX_START; the amount of space in X. ! ! Output, integer NSTART; the number of starting points, which will ! never be more than MAX_START. ! ! Output, real ( kind = 8 ) X(MAX_START), contains NSTART starting points. ! implicit none integer max_start integer nstart real ( kind = 8 ) x(max_start) nstart = 0 if ( nstart < max_start ) then nstart = nstart + 1 x(nstart) = 0.0002D+00 end if if ( nstart < max_start ) then nstart = nstart + 1 x(nstart) = 2.0D+00 end if return end subroutine p06_title ( title ) !*****************************************************************************80 ! !! P06_TITLE returns the title of problem 6. ! ! Modified: ! ! 07 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, character ( len = * ) TITLE, the title of the problem. ! implicit none character ( len = * ) title title = 'F(X) = EXP(X) - 2 - 1 / ( 10 * X )**2 - 2 / ( 100 * X )**3' return end subroutine p07_fx ( x, fx ) !*****************************************************************************80 ! !! P07_FX evaluates x^3. ! ! Modified: ! ! 07 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the point at which F is to be evaluated. ! ! Output, real ( kind = 8 ) FX, the value of the function at X. ! implicit none real ( kind = 8 ) fx real ( kind = 8 ) x fx = x**3 return end subroutine p07_fx1 ( x, fx1 ) !*****************************************************************************80 ! !! P07_FX1 evaluates the derivative of the function for problem 7. ! ! Modified: ! ! 08 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the abscissa. ! ! Output, real ( kind = 8 ) FX1, the first derivative of the function at X. ! implicit none real ( kind = 8 ) fx1 real ( kind = 8 ) x fx1 = 3.0D+00 * x**2 return end subroutine p07_fx2 ( x, fx2 ) !*****************************************************************************80 ! !! P07_FX2 evaluates the second derivative of the function for problem 7. ! ! Modified: ! ! 08 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the abscissa. ! ! Output, real ( kind = 8 ) FX2, the second derivative of the function at X. ! implicit none real ( kind = 8 ) fx2 real ( kind = 8 ) x fx2 = 6.0D+00 * x return end subroutine p07_interval ( xmin, xmax ) !*****************************************************************************80 ! !! P07_INTERVAL returns an interval bounding the root for problem 7. ! ! Modified: ! ! 07 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real ( kind = 8 ) XMIN, XMAX, the minimum and maximum values of ! an interval containing the root. ! implicit none real ( kind = 8 ) xmax real ( kind = 8 ) xmin xmin = -1000.0D+00 xmax = 1000.0D+00 return end subroutine p07_root ( max_root, nroot, x ) !*****************************************************************************80 ! !! P07_ROOT returns a root for problem 7. ! ! Modified: ! ! 07 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAX_ROOT, the maximum number of roots that can be stored. ! ! Output, integer NROOT, the number of roots returned. ! ! Output, real ( kind = 8 ) X(MAX_ROOT), NROOT roots for the given problem. ! implicit none integer max_root integer nroot real ( kind = 8 ) x(max_root) nroot = 0 if ( nroot < max_root ) then nroot = nroot + 1 x(nroot) = 0.0D+00 end if return end subroutine p07_start ( max_start, nstart, x ) !*****************************************************************************80 ! !! P07_START returns one or more starting points for problem 7. ! ! Modified: ! ! 15 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAX_START; the amount of space in X. ! ! Output, integer NSTART; the number of starting points, which will ! never be more than MAX_START. ! ! Output, real ( kind = 8 ) X(MAX_START), contains NSTART starting points. ! implicit none integer max_start integer nstart real ( kind = 8 ) x(max_start) nstart = 0 if ( nstart < max_start ) then nstart = nstart + 1 x(nstart) = 1.0D+00 end if if ( nstart < max_start ) then nstart = nstart + 1 x(nstart) = -1000.0D+00 end if return end subroutine p07_title ( title ) !*****************************************************************************80 ! !! P07_TITLE returns the title of problem 7. ! ! Modified: ! ! 07 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, character ( len = * ) TITLE, the title of the problem. ! implicit none character ( len = * ) title title = 'F(X) = X**3, only linear Newton convergence.' return end subroutine p08_fx ( x, fx ) !*****************************************************************************80 ! !! P08_FX evaluates cos ( x ) - x. ! ! Modified: ! ! 11 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the point at which F is to be evaluated. ! ! Output, real ( kind = 8 ) FX, the value of the function at X. ! implicit none real ( kind = 8 ) fx real ( kind = 8 ) x fx = cos ( x ) - x return end subroutine p08_fx1 ( x, fx1 ) !*****************************************************************************80 ! !! P08_FX1 evaluates the derivative of the function for problem 8. ! ! Modified: ! ! 11 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the abscissa. ! ! Output, real ( kind = 8 ) FX1, the first derivative of the function at X. ! implicit none real ( kind = 8 ) fx1 real ( kind = 8 ) x fx1 = - sin ( x ) - 1.0D+00 return end subroutine p08_fx2 ( x, fx2 ) !*****************************************************************************80 ! !! P08_FX2 evaluates the second derivative of the function for problem 8. ! ! Modified: ! ! 11 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the abscissa. ! ! Output, real ( kind = 8 ) FX2, the second derivative of the function at X. ! implicit none real ( kind = 8 ) fx2 real ( kind = 8 ) x fx2 = - cos ( x ) return end subroutine p08_interval ( xmin, xmax ) !*****************************************************************************80 ! !! P08_INTERVAL returns an interval bounding the root for problem 8. ! ! Modified: ! ! 11 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real ( kind = 8 ) XMIN, XMAX, the minimum and maximum values of ! an interval containing the root. ! implicit none real ( kind = 8 ) xmax real ( kind = 8 ) xmin xmin = - 10.0D+00 xmax = 10.0D+00 return end subroutine p08_root ( max_root, nroot, x ) !*****************************************************************************80 ! !! P08_ROOT returns a root for problem 8. ! ! Modified: ! ! 11 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAX_ROOT, the maximum number of roots that can be stored. ! ! Output, integer NROOT, the number of roots returned. ! ! Output, real ( kind = 8 ) X(MAX_ROOT), NROOT roots for the given problem. ! implicit none integer max_root integer nroot real ( kind = 8 ) x(max_root) nroot = 0 if ( nroot < max_root ) then nroot = nroot + 1 x(nroot) = 0.73908514D+00 end if return end subroutine p08_start ( max_start, nstart, x ) !*****************************************************************************80 ! !! P08_START returns one or more starting points for problem 8. ! ! Modified: ! ! 15 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAX_START; the amount of space in X. ! ! Output, integer NSTART; the number of starting points, which will ! never be more than MAX_START. ! ! Output, real ( kind = 8 ) X(MAX_START), contains NSTART starting points. ! implicit none integer max_start integer nstart real ( kind = 8 ) x(max_start) nstart = 0 if ( nstart < max_start ) then nstart = nstart + 1 x(nstart) = 1.0D+00 end if if ( nstart < max_start ) then nstart = nstart + 1 x(nstart) = 0.5D+00 end if if ( nstart < max_start ) then nstart = nstart + 1 x(nstart) = - 1.6D+00 end if return end subroutine p08_title ( title ) !*****************************************************************************80 ! !! P08_TITLE returns the title of problem 8. ! ! Modified: ! ! 02 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, character ( len = * ) TITLE, the title of the problem. ! implicit none character ( len = * ) title title = 'F(X) = COS(X) - X' return end subroutine p09_fx ( x, fx ) !*****************************************************************************80 ! !! P09_FX evaluates the Newton Baffler. ! ! Modified: ! ! 02 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the point at which F is to be evaluated. ! ! Output, real ( kind = 8 ) FX, the value of the function at X. ! implicit none integer, parameter :: max_root = 1 real ( kind = 8 ) fx integer nroot real ( kind = 8 ) root(max_root) real ( kind = 8 ) x real ( kind = 8 ) x2 call p09_root ( max_root, nroot, root ) x2 = x - root(1) if ( x2 < - 0.25D+00 ) then fx = 0.75D+00 * x2 - 0.3125D+00 else if ( x2 < 0.25D+00 ) then fx = 2.0D+00 * x2 else fx = 0.75D+00 * x2 + 0.3125D+00 end if return end subroutine p09_fx1 ( x, fx1 ) !*****************************************************************************80 ! !! P09_FX1 evaluates the derivative of the function for problem 9. ! ! Modified: ! ! 02 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the abscissa. ! ! Output, real ( kind = 8 ) FX1, the first derivative of the function at X. ! implicit none integer, parameter :: max_root = 1 real ( kind = 8 ) fx1 integer nroot real ( kind = 8 ) root(max_root) real ( kind = 8 ) x real ( kind = 8 ) x2 call p09_root ( max_root, nroot, root ) x2 = x - root(1) if ( x2 < - 0.25D+00 ) then fx1 = 0.75D+00 else if ( x2 < 0.25D+00 ) then fx1 = 2.0D+00 else fx1 = 0.75D+00 end if return end subroutine p09_fx2 ( x, fx2 ) !*****************************************************************************80 ! !! P09_FX2 evaluates the second derivative of the function for problem 9. ! ! Modified: ! ! 02 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the abscissa. ! ! Output, real ( kind = 8 ) FX2, the second derivative of the function at X. ! implicit none real ( kind = 8 ) fx2 real ( kind = 8 ) x fx2 = 0.0D+00 return end subroutine p09_interval ( xmin, xmax ) !*****************************************************************************80 ! !! P09_INTERVAL returns an interval bounding the root for problem 9. ! ! Modified: ! ! 02 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real ( kind = 8 ) XMIN, XMAX, the minimum and maximum values of ! an interval containing the root. ! implicit none integer, parameter :: max_root = 1 integer nroot real ( kind = 8 ) root(max_root) real ( kind = 8 ) xmax real ( kind = 8 ) xmin call p09_root ( max_root, nroot, root ) xmin = anint ( root(1) - 10.0D+00 ) xmax = anint ( root(1) + 10.0D+00 ) return end subroutine p09_root ( max_root, nroot, x ) !*****************************************************************************80 ! !! P09_ROOT returns a root for problem 9. ! ! Modified: ! ! 02 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAX_ROOT, the maximum number of roots that can be stored. ! ! Output, integer NROOT, the number of roots returned. ! ! Output, real ( kind = 8 ) X(MAX_ROOT), NROOT roots for the given problem. ! implicit none integer max_root integer nroot real ( kind = 8 ) x(max_root) nroot = 0 if ( nroot < max_root ) then nroot = nroot + 1 x(nroot) = 6.25D+00 end if return end subroutine p09_start ( max_start, nstart, x ) !*****************************************************************************80 ! !! P09_START returns one or more starting points for problem 9. ! ! Modified: ! ! 15 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAX_START; the amount of space in X. ! ! Output, integer NSTART; the number of starting points, which will ! never be more than MAX_START. ! ! Output, real ( kind = 8 ) X(MAX_START), contains NSTART starting points. ! implicit none integer max_start integer, parameter :: max_root = 1 integer nroot integer nstart real ( kind = 8 ) root(max_root) real ( kind = 8 ) x(max_start) call p09_root ( max_root, nroot, root ) nstart = 0 if ( nstart < max_start ) then nstart = nstart + 1 x(nstart) = root(1) + 5.0D+00 end if if ( nstart < max_start ) then nstart = nstart + 1 x(nstart) = root(1) - 1.0D+00 end if if ( nstart < max_start ) then nstart = nstart + 1 x(nstart) = root(1) + 0.1D+00 end if return end subroutine p09_title ( title ) !*****************************************************************************80 ! !! P09_TITLE returns the title of problem 9. ! ! Modified: ! ! 02 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, character ( len = * ) TITLE, the title of the problem. ! implicit none character ( len = * ) title title = 'The Newton Baffler' return end subroutine p10_fx ( x, fx ) !*****************************************************************************80 ! !! P10_FX evaluates the Repeller. ! ! Modified: ! ! 02 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the point at which F is to be evaluated. ! ! Output, real ( kind = 8 ) FX, the value of the function at X. ! implicit none real ( kind = 8 ) fx real ( kind = 8 ) x fx = 20.0D+00 * x / ( 100.0D+00 * x * x + 1.0D+00 ) return end subroutine p10_fx1 ( x, fx1 ) !*****************************************************************************80 ! !! P10_FX1 evaluates the derivative of the function for problem 10. ! ! Modified: ! ! 02 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the abscissa. ! ! Output, real ( kind = 8 ) FX1, the first derivative of the function at X. ! implicit none real ( kind = 8 ) fx1 real ( kind = 8 ) x fx1 = ( 1.0D+00 - 10.0D+00 * x ) * ( 1.0D+00 + 10.0D+00 * x ) & / ( 100.0D+00 * x * x + 1.0D+00 )**2 return end subroutine p10_fx2 ( x, fx2 ) !*****************************************************************************80 ! !! P10_FX2 evaluates the second derivative of the function for problem 10. ! ! Modified: ! ! 02 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the abscissa. ! ! Output, real ( kind = 8 ) FX2, the second derivative of the function at X. ! implicit none real ( kind = 8 ) fx2 real ( kind = 8 ) x fx2 = - 200.0D+00 * x * ( 3.0D+00 - 100.0D+00 * x**2 ) & / ( 100.0D+00 * x * x + 1.0D+00 )**3 return end subroutine p10_interval ( xmin, xmax ) !*****************************************************************************80 ! !! P10_INTERVAL returns an interval bounding the root for problem 10. ! ! Modified: ! ! 02 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real ( kind = 8 ) XMIN, XMAX, the minimum and maximum values of ! an interval containing the root. ! implicit none real ( kind = 8 ) xmax real ( kind = 8 ) xmin xmin = - 10.0D+00 xmax = + 10.0D+00 return end subroutine p10_root ( max_root, nroot, x ) !*****************************************************************************80 ! !! P10_ROOT returns a root for problem 10. ! ! Modified: ! ! 02 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAX_ROOT, the maximum number of roots that can be stored. ! ! Output, integer NROOT, the number of roots returned. ! ! Output, real ( kind = 8 ) X(MAX_ROOT), NROOT roots for the given problem. ! implicit none integer max_root integer nroot real ( kind = 8 ) x(max_root) nroot = 0 if ( nroot < max_root ) then nroot = nroot + 1 x(nroot) = 0.0D+00 end if return end subroutine p10_start ( max_start, nstart, x ) !*****************************************************************************80 ! !! P10_START returns one or more starting points for problem 10. ! ! Modified: ! ! 15 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAX_START; the amount of space in X. ! ! Output, integer NSTART; the number of starting points, which will ! never be more than MAX_START. ! ! Output, real ( kind = 8 ) X(MAX_START), contains NSTART starting points. ! implicit none integer max_start integer nstart real ( kind = 8 ) x(max_start) nstart = 0 if ( nstart < max_start ) then nstart = nstart + 1 x(nstart) = 1.0D+00 end if if ( nstart < max_start ) then nstart = nstart + 1 x(nstart) = - 0.14D+00 end if if ( nstart < max_start ) then nstart = nstart + 1 x(nstart) = 0.041D+00 end if return end subroutine p10_title ( title ) !*****************************************************************************80 ! !! P10_TITLE returns the title of problem 10. ! ! Modified: ! ! 02 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, character ( len = * ) TITLE, the title of the problem. ! implicit none character ( len = * ) title title = 'The Repeller' return end subroutine p11_fx ( x, fx ) !*****************************************************************************80 ! !! P11_FX evaluates the Pinhead. ! ! Modified: ! ! 09 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the point at which F is to be evaluated. ! ! Output, real ( kind = 8 ) FX, the value of the function at X. ! implicit none real ( kind = 8 ) epsilon real ( kind = 8 ) fx real ( kind = 8 ) x epsilon = 0.00001D+00 if ( epsilon == 0.0D+00 ) then fx = ( 16.0D+00 - x**4 ) / ( 16.0D+00 * x**4 ) else fx = ( 16.0D+00 - x**4 ) / ( 16.0D+00 * x**4 + epsilon ) end if return end subroutine p11_fx1 ( x, fx1 ) !*****************************************************************************80 ! !! P11_FX1 evaluates the derivative of the function for problem 11. ! ! Modified: ! ! 09 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the abscissa. ! ! Output, real ( kind = 8 ) FX1, the first derivative of the function at X. ! implicit none real ( kind = 8 ) epsilon real ( kind = 8 ) fx1 real ( kind = 8 ) x epsilon = 0.00001D+00 if ( epsilon == 0.0D+00 ) then fx1 = - 4.0D+00 / x**5 else fx1 = - 4.0D+00 * x**3 * ( epsilon + 256.0D+00 ) & / ( 16.0D+00 * x**4 + epsilon )**2 end if return end subroutine p11_fx2 ( x, fx2 ) !*****************************************************************************80 ! !! P11_FX2 evaluates the second derivative of the function for problem 11. ! ! Modified: ! ! 09 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the abscissa. ! ! Output, real ( kind = 8 ) FX2, the second derivative of the function at X. ! implicit none real ( kind = 8 ) epsilon real ( kind = 8 ) fx2 real ( kind = 8 ) x epsilon = 0.00001D+00 if ( epsilon == 0.0D+00 ) then fx2 = 20.0D+00 / x**6 else fx2 = - 4.0D+00 * ( epsilon + 256.0D+00 ) & * ( 3.0D+00 * epsilon - 80.0D+00 * x**4 ) * x**2 & / ( 16.0D+00 * x**4 + epsilon )**3 end if return end subroutine p11_interval ( xmin, xmax ) !*****************************************************************************80 ! !! P11_INTERVAL returns an interval bounding the root for problem 11. ! ! Modified: ! ! 02 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real ( kind = 8 ) XMIN, XMAX, the minimum and maximum values of ! an interval containing the root. ! implicit none real ( kind = 8 ) xmax real ( kind = 8 ) xmin xmin = 0.0D+00 xmax = + 10.0D+00 return end subroutine p11_root ( max_root, nroot, x ) !*****************************************************************************80 ! !! P11_ROOT returns a root for problem 11. ! ! Modified: ! ! 02 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAX_ROOT, the maximum number of roots that can be stored. ! ! Output, integer NROOT, the number of roots returned. ! ! Output, real ( kind = 8 ) X(MAX_ROOT), NROOT roots for the given problem. ! implicit none integer max_root integer nroot real ( kind = 8 ) x(max_root) nroot = 0 if ( nroot < max_root ) then nroot = nroot + 1 x(nroot) = - 2.0D+00 end if if ( nroot < max_root ) then nroot = nroot + 1 x(nroot) = 2.0D+00 end if return end subroutine p11_start ( max_start, nstart, x ) !*****************************************************************************80 ! !! P11_START returns one or more starting points for problem 11. ! ! Modified: ! ! 15 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAX_START; the amount of space in X. ! ! Output, integer NSTART; the number of starting points, which will ! never be more than MAX_START. ! ! Output, real ( kind = 8 ) X(MAX_START), contains NSTART starting points. ! implicit none integer max_start integer nstart real ( kind = 8 ) x(max_start) nstart = 0 if ( nstart < max_start ) then nstart = nstart + 1 x(nstart) = 0.25D+00 end if if ( nstart < max_start ) then nstart = nstart + 1 x(nstart) = 5.0D+00 end if if ( nstart < max_start ) then nstart = nstart + 1 x(nstart) = 1.1D+00 end if return end subroutine p11_title ( title ) !*****************************************************************************80 ! !! P11_TITLE returns the title of problem 11. ! ! Modified: ! ! 02 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, character ( len = * ) TITLE, the title of the problem. ! implicit none character ( len = * ) title title = 'The Pinhead' return end subroutine p12_fx ( x, fx ) !*****************************************************************************80 ! !! P12_FX evaluates Flat Stanley. ! ! Modified: ! ! 09 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the point at which F is to be evaluated. ! ! Output, real ( kind = 8 ) FX, the value of the function at X. ! implicit none real ( kind = 8 ) factor real ( kind = 8 ) fx real ( kind = 8 ) s real ( kind = 8 ) x real ( kind = 8 ) y factor = 1000.0D+00 if ( x == 1.0D+00 ) then fx = 0.0D+00 else y = x - 1.0D+00 s = sign ( 1.0D+00, y ) fx = s * exp ( log ( factor ) + log ( abs ( y ) ) - 1.0D+00 / y**2 ) end if return end subroutine p12_fx1 ( x, fx1 ) !*****************************************************************************80 ! !! P12_FX1 evaluates the derivative of the function for problem 12. ! ! Modified: ! ! 09 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the abscissa. ! ! Output, real ( kind = 8 ) FX1, the first derivative of the function at X. ! implicit none real ( kind = 8 ) factor real ( kind = 8 ) fx1 real ( kind = 8 ) x real ( kind = 8 ) y factor = 1000.0D+00 if ( x == 1.0D+00 ) then fx1 = 0.0D+00 else y = x - 1.0D+00 fx1 = factor * exp ( - 1.0D+00 / y**2 ) * ( y**2 + 2.0D+00 ) / y**2 end if return end subroutine p12_fx2 ( x, fx2 ) !*****************************************************************************80 ! !! P12_FX2 evaluates the second derivative of the function for problem 12. ! ! Modified: ! ! 09 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the abscissa. ! ! Output, real ( kind = 8 ) FX2, the second derivative of the function at X. ! implicit none real ( kind = 8 ) factor real ( kind = 8 ) fx2 real ( kind = 8 ) x real ( kind = 8 ) y factor = 1000.0D+00 if ( x == 1.0D+00 ) then fx2 = 0.0D+00 else y = x - 1.0D+00 fx2 = - 2.0D+00 * factor * exp ( - 1.0D+00 / y**2 ) & * ( y**2 - 2.0D+00 ) / y**5 end if return end subroutine p12_interval ( xmin, xmax ) !*****************************************************************************80 ! !! P12_INTERVAL returns an interval bounding the root for problem 12. ! ! Modified: ! ! 06 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real ( kind = 8 ) XMIN, XMAX, the minimum and maximum values of ! an interval containing the root. ! implicit none real ( kind = 8 ) xmax real ( kind = 8 ) xmin xmin = - 4.0D+00 xmax = 4.0D+00 return end subroutine p12_root ( max_root, nroot, x ) ! !*****************************************************************************80 ! !! P12_ROOT returns a root for problem 12. ! ! ! Modified: ! ! 06 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAX_ROOT, the maximum number of roots that can be stored. ! ! Output, integer NROOT, the number of roots returned. ! ! Output, real ( kind = 8 ) X(MAX_ROOT), NROOT roots for the given problem. ! implicit none integer max_root integer nroot real ( kind = 8 ) x(max_root) nroot = 0 if ( nroot < max_root ) then nroot = nroot + 1 x(nroot) = 1.0D+00 end if return end subroutine p12_start ( max_start, nstart, x ) !*****************************************************************************80 ! !! P12_START returns one or more starting points for problem 12. ! ! Modified: ! ! 15 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAX_START; the amount of space in X. ! ! Output, integer NSTART; the number of starting points, which will ! never be more than MAX_START. ! ! Output, real ( kind = 8 ) X(MAX_START), contains NSTART starting points. ! implicit none integer max_start integer nstart real ( kind = 8 ) x(nstart) nstart = 0 if ( nstart < max_start ) then nstart = nstart + 1 x(nstart) = 2.0D+00 end if if ( nstart < max_start ) then nstart = nstart + 1 x(nstart) = 0.50D+00 end if if ( nstart < max_start ) then nstart = nstart + 1 x(nstart) = 4.0D+00 end if return end subroutine p12_title ( title ) !*****************************************************************************80 ! !! P12_TITLE returns the title of problem 12. ! ! Modified: ! ! 09 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, character ( len = * ) TITLE, the title of the problem. ! implicit none character ( len = * ) title title = 'Flat Stanley (ALL derivatives are zero at the root.)' return end subroutine p13_fx ( x, fx ) !*****************************************************************************80 ! !! P13_FX evaluates Lazy Boy. ! ! Modified: ! ! 09 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the point at which F is to be evaluated. ! ! Output, real ( kind = 8 ) FX, the value of the function at X. ! implicit none integer, parameter :: max_root = 1 real ( kind = 8 ) fx integer nroot real ( kind = 8 ) root(max_root) real ( kind = 8 ) slope real ( kind = 8 ) x slope = 0.00000000001D+00 call p13_root ( max_root, nroot, root ) fx = slope * ( x - root(1) ) return end subroutine p13_fx1 ( x, fx1 ) !*****************************************************************************80 ! !! P13_FX1 evaluates the derivative of the function for problem 13. ! ! Modified: ! ! 09 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the abscissa. ! ! Output, real ( kind = 8 ) FX1, the first derivative of the function at X. ! implicit none real ( kind = 8 ) fx1 real ( kind = 8 ) slope real ( kind = 8 ) x slope = 0.00000000001D+00 fx1 = slope return end subroutine p13_fx2 ( x, fx2 ) !*****************************************************************************80 ! !! P13_FX2 evaluates the second derivative of the function for problem 13. ! ! Modified: ! ! 09 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the abscissa. ! ! Output, real ( kind = 8 ) FX2, the second derivative of the function at X. ! implicit none real ( kind = 8 ) fx2 real ( kind = 8 ) x fx2 = 0.0D+00 return end subroutine p13_interval ( xmin, xmax ) !*****************************************************************************80 ! !! P13_INTERVAL returns an interval bounding the root for problem 13. ! ! Modified: ! ! 09 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real ( kind = 8 ) XMIN, XMAX, the minimum and maximum values of ! an interval containing the root. ! implicit none real ( kind = 8 ) xmax real ( kind = 8 ) xmin xmin = - 10000000000000.0D+00 xmax = 10000000000000.0D+00 return end subroutine p13_root ( max_root, nroot, x ) !*****************************************************************************80 ! !! P13_ROOT returns a root for problem 13. ! ! Modified: ! ! 09 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAX_ROOT, the maximum number of roots that can be stored. ! ! Output, integer NROOT, the number of roots returned. ! ! Output, real ( kind = 8 ) X(MAX_ROOT), NROOT roots for the given problem. ! implicit none integer max_root integer nroot real ( kind = 8 ) x(max_root) nroot = 0 if ( nroot < max_root ) then nroot = nroot + 1 x(nroot) = 100.0D+00 end if return end subroutine p13_start ( max_start, nstart, x ) !*****************************************************************************80 ! !! P13_START returns one or more starting points for problem 13. ! ! Modified: ! ! 15 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAX_START; the amount of space in X. ! ! Output, integer NSTART; the number of starting points, which will ! never be more than MAX_START. ! ! Output, real ( kind = 8 ) X(MAX_START), contains NSTART starting points. ! implicit none integer max_start integer nstart real ( kind = 8 ) x(max_start) nstart = 0 if ( nstart < max_start ) then nstart = nstart + 1 x(nstart) = 100000000.0D+00 end if if ( nstart < max_start ) then nstart = nstart + 1 x(nstart) = 100000013.0D+00 end if if ( nstart < max_start ) then nstart = nstart + 1 x(nstart) = - 100000000000.0D+00 end if return end subroutine p13_title ( title ) !*****************************************************************************80 ! !! P13_TITLE returns the title of problem 13. ! ! Modified: ! ! 09 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, character ( len = * ) TITLE, the title of the problem. ! implicit none character ( len = * ) title title = 'Lazy Boy (Linear function, almost flat.)' return end subroutine p14_fx ( x, fx ) !*****************************************************************************80 ! !! P14_FX evaluates the Camel. ! ! Modified: ! ! 15 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the point at which F is to be evaluated. ! ! Output, real ( kind = 8 ) FX, the value of the function at X. ! implicit none real ( kind = 8 ) fx real ( kind = 8 ) x fx = 1.0D+00 / ( ( x - 0.3D+00 )**2 + 0.01D+00 ) & + 1.0D+00 / ( ( x - 0.9D+00 )**2 + 0.04D+00 ) + 2.0D+00 * x - 5.2D+00 return end subroutine p14_fx1 ( x, fx1 ) !*****************************************************************************80 ! !! P14_FX1 evaluates the derivative of the function for problem 14. ! ! Modified: ! ! 15 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the abscissa. ! ! Output, real ( kind = 8 ) FX1, the first derivative of the function at X. ! implicit none real ( kind = 8 ) fx1 real ( kind = 8 ) x fx1 = - 2.0D+00 * ( x - 0.3D+00 ) / ( ( x - 0.3D+00 )**2 + 0.01D+00 )**2 & - 2.0D+00 * ( x - 0.9D+00 ) / ( ( x - 0.9D+00 )**2 + 0.04D+00 )**2 & + 2.0D+00 return end subroutine p14_fx2 ( x, fx2 ) !*****************************************************************************80 ! !! P14_FX2 evaluates the second derivative of the function for problem 14. ! ! Modified: ! ! 15 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the abscissa. ! ! Output, real ( kind = 8 ) FX2, the second derivative of the function at X. ! implicit none real ( kind = 8 ) fx2 real ( kind = 8 ) x fx2 = 0.0D+00 return end subroutine p14_interval ( xmin, xmax ) !*****************************************************************************80 ! !! P14_INTERVAL returns an interval bounding the root for problem 14. ! ! Modified: ! ! 15 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real ( kind = 8 ) XMIN, XMAX, the minimum and maximum values of ! an interval containing the root. ! implicit none real ( kind = 8 ) xmax real ( kind = 8 ) xmin xmin = - 10.0D+00 xmax = 10.0D+00 return end subroutine p14_root ( max_root, nroot, x ) !*****************************************************************************80 ! !! P14_ROOT returns a root for problem 14. ! ! Modified: ! ! 15 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAX_ROOT, the maximum number of roots that can be stored. ! ! Output, integer NROOT, the number of roots returned. ! ! Output, real ( kind = 8 ) X(MAX_ROOT), NROOT roots for the given problem. ! implicit none integer max_root integer nroot real ( kind = 8 ) x(max_root) nroot = 0 if ( nroot < max_root ) then nroot = nroot + 1 x(nroot) = - 0.15348049D+00 end if if ( nroot < max_root ) then nroot = nroot + 1 x(nroot) = 1.8190323D+00 end if if ( nroot < max_root ) then nroot = nroot + 1 x(nroot) = 2.1274323D+00 end if return end subroutine p14_start ( max_start, nstart, x ) !*****************************************************************************80 ! !! P14_START returns one or more starting points for problem 14. ! ! Modified: ! ! 15 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAX_START; the amount of space in X. ! ! Output, integer NSTART; the number of starting points, which will ! never be more than MAX_START. ! ! Output, real ( kind = 8 ) X(MAX_START), contains NSTART starting points. ! implicit none integer max_start integer nstart real ( kind = 8 ) x(max_start) nstart = 0 if ( nstart < max_start ) then nstart = nstart + 1 x(nstart) = 3.0D+00 end if if ( nstart < max_start ) then nstart = nstart + 1 x(nstart) = - 0.5D+00 end if if ( nstart < max_start ) then nstart = nstart + 1 x(nstart) = 0.0D+00 end if if ( nstart < max_start ) then nstart = nstart + 1 x(nstart) = 2.12742D+00 end if return end subroutine p14_title ( title ) !*****************************************************************************80 ! !! P14_TITLE returns the title of problem 14. ! ! Modified: ! ! 15 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, character ( len = * ) TITLE, the title of the problem. ! implicit none character ( len = * ) title title = 'The Camel (double hump and some shallow roots.)' return end subroutine p15_fx ( x, fx ) !*****************************************************************************80 ! !! P15_FX evaluates a pathological function for Newton's method. ! ! Modified: ! ! 10 October 2000 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! George Donovan, Arnold Miller, Timothy Moreland, ! Pathological Functions for Newton's Method, ! American Mathematical Monthly, January 1993, pages 53-58. ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the point at which F is to be evaluated. ! ! Output, real ( kind = 8 ) FX, the value of the function at X. ! implicit none real ( kind = 8 ) fx real ( kind = 8 ) r8_cube_root real ( kind = 8 ) x fx = r8_cube_root ( x ) * exp ( - x**2 ) return end subroutine p15_fx1 ( x, fx1 ) !*****************************************************************************80 ! !! P15_FX1 evaluates the derivative of the function for problem 15. ! ! Modified: ! ! 10 October 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the abscissa. ! ! Output, real ( kind = 8 ) FX1, the first derivative of the function at X. ! implicit none real ( kind = 8 ) fx1 real ( kind = 8 ) r8_cube_root real ( kind = 8 ) x fx1 = ( 1.0D+00 - 6.0D+00 * x**2 ) * r8_cube_root ( x ) & * exp ( - x**2 ) / ( 3.0D+00 * x ) return end subroutine p15_fx2 ( x, fx2 ) !*****************************************************************************80 ! !! P15_FX2 evaluates the second derivative of the function for problem 15. ! ! Modified: ! ! 10 October 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the abscissa. ! ! Output, real ( kind = 8 ) FX2, the second derivative of the function at X. ! implicit none real ( kind = 8 ) fx2 real ( kind = 8 ) r8_cube_root real ( kind = 8 ) x fx2 = ( - 2.0D+00 - 30.0D+00 * x**2 + 36.0D+00 * x**4 ) * r8_cube_root ( x ) & * exp ( - x**2 ) / ( 9.0D+00 * x**2 ) return end subroutine p15_interval ( xmin, xmax ) !*****************************************************************************80 ! !! P15_INTERVAL returns an interval bounding the root for problem 15. ! ! Modified: ! ! 10 October 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real ( kind = 8 ) XMIN, XMAX, the minimum and maximum values of ! an interval containing the root. ! implicit none real ( kind = 8 ) xmax real ( kind = 8 ) xmin xmin = - 10.0D+00 xmax = 10.0D+00 return end subroutine p15_root ( max_root, nroot, x ) !*****************************************************************************80 ! !! P15_ROOT returns a root for problem 15. ! ! Modified: ! ! 10 October 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAX_ROOT, the maximum number of roots that can be stored. ! ! Output, integer NROOT, the number of roots returned. ! ! Output, real ( kind = 8 ) X(MAX_ROOT), NROOT roots for the given problem. ! implicit none integer max_root integer nroot real ( kind = 8 ) x(max_root) nroot = 0 if ( 1 <= max_root ) then nroot = 1 x(nroot) = 0.0D+00 end if return end subroutine p15_start ( max_start, nstart, x ) !*****************************************************************************80 ! !! P15_START returns one or more starting points for problem 15. ! ! Modified: ! ! 10 October 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAX_START; the amount of space in X. ! ! Output, integer NSTART; the number of starting points, which will ! never be more than MAX_START. ! ! Output, real ( kind = 8 ) X(MAX_START), contains NSTART starting points. ! implicit none integer max_start integer nstart real ( kind = 8 ) x(max_start) nstart = 0 if ( 1 <= max_start ) then nstart = nstart + 1 x(nstart) = 0.01D+00 end if if ( 2 <= max_start ) then nstart = nstart + 1 x(nstart) = - 0.25D+00 end if return end subroutine p15_title ( title ) !*****************************************************************************80 ! !! P15_TITLE returns the title of problem 15. ! ! Modified: ! ! 10 October 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, character ( len = * ) TITLE, the title of the problem. ! implicit none character ( len = * ) title title = 'Donovan/Miller/Moreland Pathological Function' return end subroutine p16_fx ( x, fx ) !*****************************************************************************80 ! !! P16_FX evaluates Kepler's Equation. ! ! Discussion: ! ! This is Kepler's equation. The equation has the form: ! ! X = M + E * sin ( X ) ! ! X represents the eccentric anomaly of a planet, the angle between the ! perihelion (the point on the orbit nearest to the sun) ! through the sun to the center of the ellipse, and the ! line from the center of the ellipse to the planet. ! ! There are two parameters: ! ! E is the eccentricity of the orbit, which should be between 0 and 1.0; ! ! M is the angle from the perihelion made by a fictitious planet traveling ! on a circular orbit centered at the sun, and traveling at a constant ! angular velocity equal to the average angular velocity of the true planet. ! M is usually between 0 and 180 degrees, but can have any value. ! ! For convenience, X and M are measured in degrees. ! ! Sample results: ! ! E M X ! ----- --- ---------- ! 0.100 5 5.554589 ! 0.200 5 6.246908 ! 0.300 5 7.134960 ! 0.400 5 8.313903 ! 0.500 5 9.950063 ! 0.600 5 12.356653 ! 0.700 5 16.167990 ! 0.800 5 22.656579 ! 0.900 5 33.344447 ! 0.990 5 45.361023 ! 0.990 1 24.725822 ! 0.990 33 89.722155 ! 0.750 70 110.302 ! 0.990 2 32.361007 ! ! Modified: ! ! 11 February 2001 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Peter Colwell, ! Solving Kepler's Equation Over Three Centuries, ! Willmann-Bell, 1993 ! ! Jean Meeus, ! Astronomical Algorithms, ! Willman-Bell, Inc, 1991, ! QB51.3.E43M42 ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the point at which F is to be evaluated. ! ! Output, real ( kind = 8 ) FX, the value of the function at X. ! implicit none real ( kind = 8 ) e real ( kind = 8 ) fx real ( kind = 8 ) m real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00 real ( kind = 8 ) x call p16_p_get ( e, m ) fx = ( pi * ( x - m ) / 180.0D+00 ) - e * sin ( pi * x / 180.0D+00 ) return end subroutine p16_fx1 ( x, fx1 ) !*****************************************************************************80 ! !! P16_FX1 evaluates the derivative of the function for problem 16. ! ! Modified: ! ! 11 February 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the abscissa. ! ! Output, real ( kind = 8 ) FX1, the first derivative of the function at X. ! implicit none real ( kind = 8 ) e real ( kind = 8 ) fx1 real ( kind = 8 ) m real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00 real ( kind = 8 ) x call p16_p_get ( e, m ) fx1 = ( pi / 180.0D+00 ) & - e * pi * cos ( pi * x / 180.0D+00 ) / 180.0D+00 return end subroutine p16_fx2 ( x, fx2 ) !*****************************************************************************80 ! !! P16_FX2 evaluates the second derivative of the function for problem 16. ! ! Modified: ! ! 11 February 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the abscissa. ! ! Output, real ( kind = 8 ) FX2, the second derivative of the function at X. ! implicit none real ( kind = 8 ) e real ( kind = 8 ) fx2 real ( kind = 8 ) m real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00 real ( kind = 8 ) x call p16_p_get ( e, m ) fx2 = e * pi**2 * sin ( pi * x / 180.0D+00 ) / 180.0D+00**2 return end subroutine p16_interval ( xmin, xmax ) !*****************************************************************************80 ! !! P16_INTERVAL returns an interval bounding the root for problem 16. ! ! Modified: ! ! 11 February 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real ( kind = 8 ) XMIN, XMAX, the minimum and maximum values of ! an interval containing the root. ! implicit none real ( kind = 8 ) e real ( kind = 8 ) m real ( kind = 8 ) xmax real ( kind = 8 ) xmin call p16_p_get ( e, m ) xmin = m - 180.0D+00 xmax = m + 180.0D+00 return end subroutine p16_p ( action, e, m ) !*****************************************************************************80 ! !! P16_P sets, gets, or prints the parameter values for problem 16. ! ! Modified: ! ! 12 February 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) ACTION, is 'GET', 'SET', or PRINT. ! ! Input/output, real ( kind = 8 ) E, M, the current values of E and M. ! implicit none character ( len = * ) action real ( kind = 8 ) e real ( kind = 8 ), save :: e_save = 0.8D+00 real ( kind = 8 ) m real ( kind = 8 ), save :: m_save = 5.0D+00 if ( action(1:1) == 'G' .or. action(1:1) == 'g' ) then e = e_save m = m_save else if ( action(1:1) == 'S' .or. action(1:1) == 's' ) then e_save = e m_save = m else if ( action(1:1) == 'P' .or. action(1:1) == 'p' ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'P16_P - Current parameter values:' write ( *, '(a,g14.6)' ) ' E, orbital eccentricity, = ', e_save write ( *, '(a,g14.6)' ) ' M, some angle in degrees, = ', m_save else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'P16_P - Fatal error!' write ( *, '(a)' ) ' Illegal input value of ACTION = ' // trim ( action ) stop end if return end subroutine p16_p_get ( e, m ) !*****************************************************************************80 ! !! P16_P_GET gets the parameter values for problem 16. ! ! Modified: ! ! 11 February 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real ( kind = 8 ) E, M, the current values of E and M. ! implicit none real ( kind = 8 ) e real ( kind = 8 ) m call p16_p ( 'GET', e, m ) return end subroutine p16_p_print ( ) !*****************************************************************************80 ! !! P16_P_PRINT prints the parameter values for problem 16. ! ! Modified: ! ! 12 February 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! implicit none real ( kind = 8 ) e real ( kind = 8 ) m call p16_p ( 'PRINT', e, m ) return end subroutine p16_p_set ( e, m ) !*****************************************************************************80 ! !! P16_P_SET sets the parameter values for problem 16. ! ! Modified: ! ! 11 February 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) E, M, the new values of E and M. ! implicit none real ( kind = 8 ) e real ( kind = 8 ) m call p16_p ( 'SET', e, m ) return end subroutine p16_root ( max_root, nroot, x ) !*****************************************************************************80 ! !! P16_ROOT returns a root for problem 16. ! ! Modified: ! ! 11 February 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAX_ROOT, the maximum number of roots that can be stored. ! ! Output, integer NROOT, the number of roots returned. ! ! Output, real ( kind = 8 ) X(MAX_ROOT), NROOT roots for the given problem. ! implicit none integer max_root integer nroot real ( kind = 8 ) x(max_root) nroot = 0 return end subroutine p16_start ( max_start, nstart, x ) !*****************************************************************************80 ! !! P16_START returns one or more starting points for problem 16. ! ! Modified: ! ! 11 February 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAX_START; the amount of space in X. ! ! Output, integer NSTART; the number of starting points, which will ! never be more than MAX_START. ! ! Output, real ( kind = 8 ) X(MAX_START), contains NSTART starting points. ! implicit none integer max_start real ( kind = 8 ) e real ( kind = 8 ) m integer nstart real ( kind = 8 ) x(max_start) call p16_p_get ( e, m ) nstart = 0 if ( 1 <= max_start ) then nstart = nstart + 1 x(nstart) = 0.0D+00 end if if ( 2 < max_start ) then nstart = nstart + 1 x(nstart) = m end if if ( 3 < max_start ) then nstart = nstart + 1 x(nstart) = m + 180.0D+00 end if return end subroutine p16_title ( title ) !*****************************************************************************80 ! !! P16_TITLE returns the title of problem 16. ! ! Modified: ! ! 11 February 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, character ( len = * ) TITLE, the title of the problem. ! implicit none character ( len = * ) title title = 'Kepler''s Eccentric Anomaly Equation, in degrees' return end subroutine p17_fx ( x, fx ) !*****************************************************************************80 ! !! P17_FX evaluates the bogus function (returns random numbers). ! ! Discussion: ! ! This is a bogus function, which simply returns random numbers ! for all requested quantities. ! ! An intelligent zero finder ought to realize, at some point, ! that it is being lied to. Or at the very least, that the ! function behavior is beyond its ability to handle. ! ! Modified: ! ! 17 December 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the point at which F is to be evaluated. ! ! Output, real ( kind = 8 ) FX, the value of the function at X. ! implicit none real ( kind = 8 ) fx real ( kind = 8 ) x call random_number ( harvest = fx ) return end subroutine p17_fx1 ( x, fx1 ) !*****************************************************************************80 ! !! P17_FX1 evaluates the derivative of the function for problem 17. ! ! Modified: ! ! 17 December 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the abscissa. ! ! Output, real ( kind = 8 ) FX1, the first derivative of the function at X. ! implicit none real ( kind = 8 ) fx1 real ( kind = 8 ) x call random_number ( harvest = fx1 ) return end subroutine p17_fx2 ( x, fx2 ) !*****************************************************************************80 ! !! P17_FX2 evaluates the second derivative of the function for problem 17. ! ! Modified: ! ! 17 December 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the abscissa. ! ! Output, real ( kind = 8 ) FX2, the second derivative of the function at X. ! implicit none real ( kind = 8 ) fx2 real ( kind = 8 ) x call random_number ( harvest = fx2 ) return end subroutine p17_interval ( xmin, xmax ) !*****************************************************************************80 ! !! P17_INTERVAL returns an interval bounding the root for problem 17. ! ! Modified: ! ! 17 December 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real ( kind = 8 ) XMIN, XMAX, the minimum and maximum values of ! an interval containing the root. ! implicit none real ( kind = 8 ) xmax real ( kind = 8 ) xmin xmin = -huge ( xmin ) xmax = huge ( xmax ) return end subroutine p17_root ( max_root, nroot, x ) !*****************************************************************************80 ! !! P17_ROOT returns a root for problem 17. ! ! Modified: ! ! 17 December 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAX_ROOT, the maximum number of roots that can be stored. ! ! Output, integer NROOT, the number of roots returned. ! ! Output, real ( kind = 8 ) X(MAX_ROOT), NROOT roots for the given problem. ! implicit none integer max_root integer nroot real ( kind = 8 ) x(max_root) if ( 1 <= max_root )then nroot = 1 call random_number ( harvest = x(1:nroot) ) end if return end subroutine p17_start ( max_start, nstart, x ) !*****************************************************************************80 ! !! P17_START returns one or more starting points for problem 17. ! ! Modified: ! ! 17 December 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAX_START; the amount of space in X. ! ! Output, integer NSTART; the number of starting points, which will ! never be more than MAX_START. ! ! Output, real ( kind = 8 ) X(MAX_START), contains NSTART starting points. ! implicit none integer max_start integer nstart real ( kind = 8 ) x(max_start) nstart = max_start call random_number ( harvest = x(1:max_start) ) return end subroutine p17_title ( title ) !*****************************************************************************80 ! !! P17_TITLE returns the title of problem 17. ! ! Modified: ! ! 17 December 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, character ( len = * ) TITLE, the title of the problem. ! implicit none character ( len = * ) title title = 'The Bogus Function, (Random Values for All Quantities)' return end subroutine p18_fx ( x, fx ) !*****************************************************************************80 ! !! P18_FX evaluates Wallis's function, f(x) = x^3 - 2*x - 5. ! ! Discussion: ! ! This simple example is of historical interest, since it was used ! by Wallis to illustrate the use of Newton's method, and has been ! a common example ever since. ! ! Modified: ! ! 18 December 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the point at which F is to be evaluated. ! ! Output, real ( kind = 8 ) FX, the value of the function at X. ! implicit none real ( kind = 8 ) fx real ( kind = 8 ) x fx = x**3 - 2.0D+00 * x - 5.0D+00 return end subroutine p18_fx1 ( x, fx1 ) !*****************************************************************************80 ! !! P18_FX1 evaluates the derivative of the function for problem 18. ! ! Modified: ! ! 18 December 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the abscissa. ! ! Output, real ( kind = 8 ) FX1, the first derivative of the function at X. ! implicit none real ( kind = 8 ) fx1 real ( kind = 8 ) x fx1 = 3.0D+00 * x**2 - 2.0D+00 return end subroutine p18_fx2 ( x, fx2 ) !*****************************************************************************80 ! !! P18_FX2 evaluates the second derivative of the function for problem 18. ! ! Modified: ! ! 18 December 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the abscissa. ! ! Output, real ( kind = 8 ) FX2, the second derivative of the function at X. ! implicit none real ( kind = 8 ) fx2 real ( kind = 8 ) x fx2 = 6.0D+00 * x return end subroutine p18_interval ( xmin, xmax ) !*****************************************************************************80 ! !! P18_INTERVAL returns an interval bounding the root for problem 18. ! ! Modified: ! ! 18 December 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real ( kind = 8 ) XMIN, XMAX, the minimum and maximum values of ! an interval containing the root. ! implicit none real ( kind = 8 ) xmax real ( kind = 8 ) xmin xmin = 2.0D+00 xmax = 3.0D+00 return end subroutine p18_root ( max_root, nroot, x ) !*****************************************************************************80 ! !! P18_ROOT returns a root for problem 18. ! ! Modified: ! ! 18 December 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAX_ROOT, the maximum number of roots that can be stored. ! ! Output, integer NROOT, the number of roots returned. ! ! Output, real ( kind = 8 ) X(MAX_ROOT), NROOT roots for the given problem. ! implicit none integer max_root integer nroot real ( kind = 8 ) x(max_root) nroot = 0 if ( 1 <= max_root ) then nroot = nroot + 1 x(nroot) = 2.0945516D+00 end if return end subroutine p18_start ( max_start, nstart, x ) !*****************************************************************************80 ! !! P18_START returns one or more starting points for problem 18. ! ! Modified: ! ! 18 December 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAX_START; the amount of space in X. ! ! Output, integer NSTART; the number of starting points, which will ! never be more than MAX_START. ! ! Output, real ( kind = 8 ) X(MAX_START), contains NSTART starting points. ! implicit none integer max_start integer nstart real ( kind = 8 ) x(max_start) nstart = 0 if ( 1 <= max_start ) then nstart = nstart + 1 x(nstart) = 2.0D+00 end if if ( 2 <= max_start ) then nstart = nstart + 1 x(nstart) = 3.0D+00 end if return end subroutine p18_title ( title ) !*****************************************************************************80 ! !! P18_TITLE returns the title of problem 18. ! ! Modified: ! ! 18 December 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, character ( len = * ) TITLE, the title of the problem. ! implicit none character ( len = * ) title title = 'The Wallis example, x^3-2x-5=0' return end subroutine p19_fx ( x, fx ) !*****************************************************************************80 ! !! P19_FX evaluates the function with a very thin pole. ! ! Discussion: ! ! This function looks positive, but has a pole at x = pi, ! and two zeroes nearby. None of this will show up computationally. ! ! Modified: ! ! 19 February 2003 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Arnold Krommer, Christoph Ueberhuber, ! Numerical Integration on Advanced Systems, ! Springer, 1994, pages 185-186. ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the point at which F is to be evaluated. ! ! Output, real ( kind = 8 ) FX, the value of the function at X. ! implicit none real ( kind = 8 ) fx real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00 real ( kind = 8 ) x fx = 3.0D+00 * x**2 + 1.0D+00 + ( log ( ( x - pi )**2 ) ) / pi**4 return end subroutine p19_fx1 ( x, fx1 ) !*****************************************************************************80 ! !! P19_FX1 evaluates the derivative of the function for problem 19. ! ! Modified: ! ! 19 February 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the abscissa. ! ! Output, real ( kind = 8 ) FX1, the first derivative of the function at X. ! implicit none real ( kind = 8 ) fx1 real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00 real ( kind = 8 ) x fx1 = 6.0D+00 * x + ( 2.0D+00 / ( x - pi ) ) / pi**4 return end subroutine p19_fx2 ( x, fx2 ) !*****************************************************************************80 ! !! P19_FX2 evaluates the second derivative of the function for problem 19. ! ! Modified: ! ! 19 February 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the abscissa. ! ! Output, real ( kind = 8 ) FX2, the second derivative of the function at X. ! implicit none real ( kind = 8 ) fx2 real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00 real ( kind = 8 ) x fx2 = 6.0D+00 + ( - 2.0D+00 / ( x - pi )**2 ) / pi**4 return end subroutine p19_interval ( xmin, xmax ) !*****************************************************************************80 ! !! P19_INTERVAL returns an interval bounding the root for problem 19. ! ! Modified: ! ! 19 February 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real ( kind = 8 ) XMIN, XMAX, the minimum and maximum values of ! an interval containing the root. ! implicit none real ( kind = 8 ) xmax real ( kind = 8 ) xmin xmin = 2.0D+00 xmax = 4.0D+00 return end subroutine p19_root ( max_root, nroot, x ) !*****************************************************************************80 ! !! P19_ROOT returns a root for problem 19. ! ! Modified: ! ! 19 February 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAX_ROOT, the maximum number of roots that can be stored. ! ! Output, integer NROOT, the number of roots returned. ! ! Output, real ( kind = 8 ) X(MAX_ROOT), NROOT roots for the given problem. ! implicit none integer max_root integer nroot real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00 real ( kind = 8 ) x(max_root) nroot = 0 if ( 1 <= max_root ) then nroot = nroot + 1 x(nroot) = ( 1.0D+00 - epsilon ( x(1) ) ) * pi end if if ( 2 <= max_root ) then nroot = nroot + 1 x(nroot) = ( 1.0D+00 + epsilon ( x(1) ) ) * pi end if return end subroutine p19_start ( max_start, nstart, x ) !*****************************************************************************80 ! !! P19_START returns one or more starting points for problem 19. ! ! Modified: ! ! 19 February 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAX_START; the amount of space in X. ! ! Output, integer NSTART; the number of starting points, which will ! never be more than MAX_START. ! ! Output, real ( kind = 8 ) X(MAX_START), contains NSTART starting points. ! implicit none integer max_start integer nstart real ( kind = 8 ) x(max_start) nstart = 0 if ( 1 <= max_start ) then nstart = nstart + 1 x(nstart) = 2.0D+00 end if return end subroutine p19_title ( title ) !*****************************************************************************80 ! !! P19_TITLE returns the title of problem 19. ! ! Modified: ! ! 19 February 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, character ( len = * ) TITLE, the title of the problem. ! implicit none character ( len = * ) title title = 'The "Thin Pole", x^2+1+log((pi-x)**2)/pi^4' return end function r8_cube_root ( x ) !*****************************************************************************80 ! !! R8_CUBE_ROOT returns the cube root of a real number. ! ! Discussion: ! ! This routine is designed to avoid the possible problems that can occur ! when formulas like 0.0**(1/3) or (-1.0)**(1/3) are to be evaluated. ! ! Modified: ! ! 01 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the number whose cube root is desired. ! ! Output, real ( kind = 8 ) R8_CUBE_ROOT, the cube root of X. ! implicit none real ( kind = 8 ) r8_cube_root real ( kind = 8 ) x if ( 0.0D+00 < x ) then r8_cube_root = x**(1.0D+00/3.0D+00) else if ( x == 0.0D+00 ) then r8_cube_root = 0.0D+00 else r8_cube_root = - ( abs ( x ) )**(1.0D+00/3.0D+00) end if return end subroutine r8_swap ( x, y ) !*****************************************************************************80 ! !! R8_SWAP switches two R8's. ! ! Modified: ! ! 30 November 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, real ( kind = 8 ) X, Y. On output, the values of X and ! Y have been interchanged. ! implicit none real ( kind = 8 ) x real ( kind = 8 ) y real ( kind = 8 ) z z = x x = y y = z return end subroutine r8poly2_rroot ( a, b, c, r1, r2 ) !*****************************************************************************80 ! !! R8POLY2_RROOT returns the real parts of the roots of a quadratic polynomial. ! ! Example: ! ! A B C roots R1 R2 ! -- -- -- ------------------ -- -- ! 1 -4 3 1 3 1 3 ! 1 0 4 2*i - 2*i 0 0 ! 2 -6 5 3 + i 3 - i 3 3 ! ! Modified: ! ! 16 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) A, B, C, the coefficients of the quadratic ! polynomial A * X**2 + B * X + C = 0 whose roots are desired. ! A must not be zero. ! ! Output, real ( kind = 8 ) R1, R2, the real parts of the two roots ! of the polynomial. ! implicit none real ( kind = 8 ) a real ( kind = 8 ) b real ( kind = 8 ) c real ( kind = 8 ) disc real ( kind = 8 ) q real ( kind = 8 ) r1 real ( kind = 8 ) r2 if ( a == 0.0D+00 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8POLY2_RROOT - Fatal error!' write ( *, '(a)' ) ' The coefficient A is zero.' stop end if disc = b * b - 4.0D+00 * a * c disc = max ( disc, 0.0D+00 ) q = ( b + sign ( 1.0D+00, b ) * sqrt ( disc ) ) r1 = - 0.5D+00 * q / a r2 = - 2.0D+00 * c / q return end subroutine regula_falsi ( fatol, step_max, nprob, xatol, xa, xb, fxa, fxb ) !*****************************************************************************80 ! !! REGULA_FALSI carries out the Regula Falsi method to seek a root of F(X) = 0. ! ! Modified: ! ! 10 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) FATOL, an absolute error tolerance for the ! function value of the root. If an approximate root X satisfies ! ABS ( F ( X ) ) <= FATOL, then X will be accepted as the ! root and the iteration will be terminated. ! ! Input, integer STEP_MAX, the maximum number of steps allowed ! for an iteration. ! ! Input, integer NPROB, the index of the function whose root is ! to be sought. ! ! Input, integer NSTART, the number of starting points supplied. ! ! Input, real ( kind = 8 ) XATOL, an absolute error tolerance for the root. ! ! Input, real ( kind = 8 ) XSTART(NSTART), several starting points for the ! iteration. ! implicit none real ( kind = 8 ) fatol real ( kind = 8 ) fxa real ( kind = 8 ) fxb real ( kind = 8 ) fxc integer step_max integer nprob integer step_num real ( kind = 8 ) xa real ( kind = 8 ) xb real ( kind = 8 ) xatol real ( kind = 8 ) xc ! ! The method requires a change-of-sign interval. ! if ( sign ( 1.0D+00, fxa ) == sign ( 1.0D+00, fxb ) ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'REGULA_FALSI - Fatal error!' write ( *, '(a)' ) ' Function does not change sign at endpoints.' stop end if ! ! Make A the root with negative F, B the root with positive F. ! if ( 0.0D+00 < fxa ) then call r8_swap ( xa, xb ) call r8_swap ( fxa, fxb ) end if write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Regula Falsi' write ( *, '(a)' ) ' ' write ( *, '(a)' ) & ' Step XA XB F(XA) F(XB)' write ( *, '(a)' ) ' ' step_num = 0 write ( *, '(2x,i4,2x,2g16.8,2g14.6)' ) step_num, xa, xb, fxa, fxb do step_num = 1, step_max if ( abs ( xa - xb ) < xatol ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'REGULA_FALSI' write ( *, '(a)' ) ' Interval small enough for convergence.' return end if if ( abs ( fxa ) <= fatol .or. abs ( fxb ) <= fatol ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'REGULA_FALSI' write ( *, '(a)' ) ' Function small enough for convergence.' return end if xc = ( fxa * xb - fxb * xa ) / ( fxa - fxb ) call p00_fx ( nprob, xc, fxc ) if ( fxc < 0.0D+00 ) then xa = xc fxa = fxc else xb = xc fxb = fxc end if write ( *, '(2x,i4,2x,2g16.8,2g14.6)' ) step_num, xa, xb, fxa, fxb end do write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'REGULA_FALSI' write ( *, '(a)' ) ' Took maximum number of steps without convergence.' return end subroutine secant ( fatol, step_max, nprob, xatol, xmin, xmax, xa, xb, fxa, & fxb ) !*****************************************************************************80 ! !! SECANT carries out the secant method to seek a root of F(X) = 0. ! ! Modified: ! ! 17 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) FATOL, an absolute error tolerance for the ! function value of the root. If an approximate root X satisfies ! ABS ( F ( X ) ) <= FATOL, then X will be accepted as the ! root and the iteration will be terminated. ! ! Input, integer STEP_MAX, the maximum number of steps allowed ! for an iteration. ! ! Input, integer NPROB, the index of the function whose root is ! to be sought. ! ! Input, integer NSTART, the number of starting points supplied. ! ! Input, real ( kind = 8 ) XATOL, an absolute error tolerance for the root. ! ! Input, real ( kind = 8 ) XMAX, XMIN, the interval in which the root should ! be sought. ! ! Input, real ( kind = 8 ) XSTART(NSTART), several starting points for the ! iteration. ! implicit none real ( kind = 8 ) fatol real ( kind = 8 ) fxa real ( kind = 8 ) fxb real ( kind = 8 ) fxc integer step_max integer nprob integer step_num real ( kind = 8 ) xa real ( kind = 8 ) xatol real ( kind = 8 ) xb real ( kind = 8 ) xc real ( kind = 8 ) xmax real ( kind = 8 ) xmin write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' SECANT METHOD' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' X F(X)' write ( *, '(a)' ) ' ' step_num = -1 write ( *, '(2x,i4,2x,g16.8,g14.6)' ) step_num, xa, fxa if ( abs ( fxa ) <= fatol ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SECANT:' write ( *, '(a)' ) ' Function small enough for convergence.' return end if step_num = 0 write ( *, '(2x,i4,2x,g16.8,g14.6)' ) step_num, xb, fxb do step_num = 1, step_max if ( abs ( fxb ) <= fatol ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SECANT:' write ( *, '(a)' ) ' Function small enough for convergence.' return end if if ( abs ( xa - xb ) < xatol ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SECANT:' write ( *, '(a)' ) ' Interval small enough for convergence.' return end if if ( xb < xmin .or. xmax < xb ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SECANT:' write ( *, '(a)' ) ' Iterate has left the region [XMIN,XMAX].' return end if if ( fxa == fxb ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SECANT:' write ( *, '(a)' ) ' F(A) = F(B), algorithm fails.' return end if xc = ( fxa * xb - fxb * xa ) / ( fxa - fxb ) call p00_fx ( nprob, xc, fxc ) xa = xb fxa = fxb xb = xc fxb = fxc write ( *, '(2x,i4,2x,g16.8,g14.6)' ) step_num, xb, fxb end do write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SECANT:' write ( *, '(a)' ) ' Took maximum number of steps.' return end subroutine timestamp ( ) !*****************************************************************************80 ! !! TIMESTAMP prints the current YMDHMS date as a time stamp. ! ! Example: ! ! 31 May 2001 9:45:54.872 AM ! ! Modified: ! ! 06 August 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! None ! implicit none character ( len = 8 ) ampm integer d integer h integer m integer mm character ( len = 9 ), parameter, dimension(12) :: month = (/ & 'January ', 'February ', 'March ', 'April ', & 'May ', 'June ', 'July ', 'August ', & 'September', 'October ', 'November ', 'December ' /) integer n integer s integer values(8) integer y call date_and_time ( values = values ) y = values(1) m = values(2) d = values(3) h = values(5) n = values(6) s = values(7) mm = values(8) if ( h < 12 ) then ampm = 'AM' else if ( h == 12 ) then if ( n == 0 .and. s == 0 ) then ampm = 'Noon' else ampm = 'PM' end if else h = h - 12 if ( h < 12 ) then ampm = 'PM' else if ( h == 12 ) then if ( n == 0 .and. s == 0 ) then ampm = 'Midnight' else ampm = 'AM' end if end if end if write ( *, '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) & d, trim ( month(m) ), y, h, ':', n, ':', s, '.', mm, trim ( ampm ) return end