Nullstellensuche mit der Regula falsi

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 mittels "
    write (*,*) "Regula Falsi 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 RegulaFalsi
  ! Nullstellensuche der Funktion f=cos(x)-x im Intervall [xu,xo]
  ! mittels Regula Falsi

  use kinds
  use mysubs
  implicit none

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

  call init(xu, xo, eps)
  fehler = abs(xo - xu)
  xneu = xu
  xdum = xo
  do
     if (fehler <= eps) exit
     xneu = xu - f(xu)*(xu - xo)/(f(xu) - f(xo))
     fehler = abs(xneu - xdum)
     if (f(xneu) == 0.0) then
        call result(xneu)
        stop
     end if
     if (f(xneu)*f(xu) < 0) then
        xo = xneu
     else
        xu = xneu
     end if
     xdum = xneu
  end do

  call result(xneu)
  stop
end program RegulaFalsi