SUBROUTINE Newton(init_X, root_eps, f_eps) INTERVAL(8), INTENT(IN) :: init_X REAL(8), INTENT(IN) :: f_eps, root_eps REAL(8) mid_of_X INTEGER, PARAMETER :: STACK_SIZE = 1000 INTERVAL(8) domain_stack(STACK_SIZE) INTEGER stack_pointer INTERVAL(8) cur_X, new_X, I_mid_of_X, I_deriv, I_val_f stack_pointer = 1 domain_stack(stack_pointer) = init_X MAIN_LOOP: DO WHILE (stack_pointer .NE. 0) cur_X = domain_stack(stack_pointer) stack_pointer = stack_pointer - 1 DO I_val_f = f_interval(cur_X) IF (.NOT. (0.0 .IN. I_val_f)) CYCLE MAIN_LOOP IF (wid(cur_X) < root_eps .AND. mag(I_val_f) < f_eps) THEN CALL print_root(cur_X) CYCLE MAIN_LOOP ENDIF ! (wid(cur_X) < root_eps ) mid_of_X = mid(cur_X) I_deriv = f_deriv(cur_X) IF (0.0 .IN. I_deriv) THEN IF (mid_of_X == inf(cur_X) .OR. mid_of_X == sup(cur_X)) & CYCLE MAIN_LOOP stack_pointer = stack_pointer + 1 IF ( stack_pointer > STACK_SIZE ) THEN PRINT*, "More than ", STACK_SIZE, " intervals required." PRINT*, "Increase PARAMETER STACK_SIZE in SUBROUTINE Newton." RETURN ENDIF domain_stack(stack_pointer) = INTERVAL(mid_of_X, sup(cur_X)) cur_X = INTERVAL(inf(cur_X), mid_of_X) CYCLE ELSE ! f' is bounded away from 0 => a Newton step I_mid_of_X = INTERVAL(mid_of_X) new_X = cur_X .IX. I_mid_of_X - f_interval(I_mid_of_X)/I_deriv IF (ISEMPTY(new_X)) THEN CYCLE MAIN_LOOP ELSEIF (cur_X == new_X) THEN CALL print_root(cur_X) CYCLE MAIN_LOOP ELSE cur_X = new_X CYCLE ENDIF ENDIF ! (0.0 .IN. I_deriv) ENDDO ENDDO MAIN_LOOP CONTAINS INTERVAL(8) FUNCTION f_interval(x) ! interval function INTERVAL(8), INTENT(IN) :: x f_interval = exp(x) - x - 5.0 END FUNCTION f_interval ! INTERVAL(8) FUNCTION f_deriv(x) ! interval first derivative INTERVAL(8), INTENT(IN) :: x f_deriv = exp(x) - 1.0 END FUNCTION f_deriv ! END SUBROUTINE Newton ! PROGRAM TEST_NEWTON INTERVAL(8) X_0 REAL(8) f_eps, root_eps X_0 = [-100.0_8, 100.0_8] f_eps = 1.0E-5_8 root_eps = 1.0E-5_8 PRINT*, "Over the interval ", X_0, "use the interval " PRINT*, "Newton algorithm to find all roots of the nonlinear " PRINT*, "function f(X) = EXP(X) - X - 5.0." PRINT*, "Roots are printed using the [inf, sup] and single-number" PRINT*, "format in which root + [-1, +1]_uld contains the root." PRINT* CALL Newton(X_0, root_eps, f_eps) END PROGRAM TEST_NEWTON ! SUBROUTINE print_root(X) INTERVAL(8), INTENT(IN) :: X WRITE(UNIT=6, FMT="(' root = ', VF40.15, ' or ', Y25.15)") X, X END SUBROUTINE print_root