Nullstellensuche durch Intervallschachtelung

module mysubs

contains

  function f(x)
    use kinds
    implicit none
    real(kind=REAL8) :: f
    real(kind=REAL8), intent(in) :: x

    f = cos(x) - x
    return
  end function f

  subroutine init(xu, xo, eps)
    use kinds
    implicit none
    
    real(kind=REAL8), intent(out) :: xu, xo, eps
    
    write (*,*) "Exemplarische Darstellung einer Nullstellensuche durch "
    write (*,*) "Intervallschachtelung anhand der Funktion f=cosinus(x)-x"
    write (*,*)
    write (*,*) "Geben Sie bitte die Grenzen des Intervalls an, in dem &
                &gesucht werden soll."
    write (*,*) "(Vorsicht, (cos(xunten)-xunten)*(cos(xoben)-xoben) muss &
                &negativ sein!)"
    
    read (*,*) xu, xo
    do
       if (f(xu)*f(xo) <= 0) exit
       write (*,*) "Zwischen diesen Grenzen liegt keine Nullstelle. &
                   &Versuchen Sie es noch einmal."
       read (*,*) xu, xo
    end do
    
    write (*,*) "Wie genau soll die Nullstelle bestimmt werden?"
    read (*,*) eps
    do
       if (eps > 0) exit
       write (*,*) "Nur positive Genauigkeiten machen Sinn. Versuchen Sie es &
                   &noch einmal."
       read (*,*) eps
    end do
    return
  end subroutine init

  subroutine result(x)
    use kinds
    implicit none
    real(kind=REAL8), intent(in) :: x
    
    write (*,'(A, F10.7)') 'Die Nullstelle liegt bei', x
    return
  end subroutine result

end module mysubs


program Intervallschachtelung
  ! Nullstellensuche der Funktion f=cos(x)-x im Intervall [xu, xo] durch
  ! Intervallschachtelung

  use kinds
  use mysubs
  implicit none

  real(kind=REAL8) :: xu, xo      ! Intervallgrenzen
  real(kind=REAL8) :: eps         ! Genauigkeit
  real(kind=REAL8) :: fehler, xneu

  call init(xu, xo, eps)

  fehler = abs(xo - xu)/2.0
  do 
     if (fehler <= eps) exit
     xneu = (xo + xu)/2.0
     if (f(xneu) == 0.0) then
        call result(xneu)
        stop
     end if
     if (f(xneu)*f(xu) < 0) then
        xo = xneu
     else
        xu = xneu
        fehler = 0.5*fehler
     end if
  end do
  xneu = (xo + xu)/2.0
  call result(xneu)
  stop
end program Intervallschachtelung