!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  PROGRAM RegFalsi
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  (c) 2005, Christian B. Lang, Mathematische Methoden in der Physik
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  Eine Nullstelle einer gegebenen Funktion
!  wird mit dem 'Regula Falsi' Verfahren gesucht.
!  Als Startwerte muessen die Schranken eines Intervalls
!  angegeben werden, in dem die Funktion ihr Vorzeichen
!  wechselt.
!  Geeignete Testwerte: a= -0.4, b=0.9 f. 1. Testfunktion
!                       a= -0.1, b=0.3 f. 2. Testfunktion  
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


PROGRAM RegFalsi

   implicit none

!  Genauigkeit wird vorgegeben
   real, parameter :: xeps = 0.000001,feps = 0.000001
   integer:: iter
   real::  x,a,b,fx,fa,fb,fun

10 iter = 0
   write(*,*) 'Ich suche eine Nullstelle ( abs(f) < ', feps,')'
   write(*,*) 'der gegebenen Funktion im Intervall [a,b] .'

   write(*,*) 'Untere Intervallgrenze a = '
   read(*,*) (a)
   fa = fun(a)
   write(*,*) 'f (',a,') = ', fa
   x = a
   if(Abs(fa) < feps)goto 20
   write(*,*) 'Obere  Intervallgrenze  b = '
   read(*,*) b
   fb = fun(b)
   x = b
   if(Abs(fb) < feps)goto 20
   write(*,*) 'f (',b,') = ',fb
   if(fa*fb >= 0)then
     write(*,*) 'Kein Vorzeichenwechsel im Intervall! Weiter mit RETURN.'
     read(*,*)
     goto 10
   end if

!  Erhoehe Iterationszaehler
30 iter = iter + 1
!  Ist das Intervall zu klein?
   if(abs(b - a) < xeps)goto 20
!  Berechne Schnittpunkt:
   x = (fb * a - fa * b) / (fb - fa)
   fx = fun(x)
   write(*,*) 'Zwischenergebnis:'
   write(*,*) 'Iteration : ',iter
   write(*,*) '        a : ',a, ' f(a) : ', fa
   write(*,*) '        x : ',x, ' f(x) : ', fx
   write(*,*) '        b : ',b, ' f(b) : ', fb

   write(*,*) '...weiter mit RETURN'
   read(*,*)
   if(abs(fx) < feps)goto 20
   if(fx*fb > 0)then
     b = x
     fb = fx
     goto 30
   end if
   if(fx*fa > 0)then
     a = x
     fa = fx
     goto 30
   end if

20 write(*,*) 'Iteration          = ',iter
   write(*,*) 'Bester Wert von x  = ',x
   write(*,*) '              f(x) = ',fx
   write(*,*) 'Grenzen fuer x : a = ',a
   write(*,*) '              f(a) = ',fa
   write(*,*) '                 b = ',b
   write(*,*) '              f(b) = ',fb

END PROGRAM regfalsi

REAL FUNCTION fun (x)

   implicit none

   real, intent(in)::x

!  Erste Testfunktion
   fun = 0.75 - 0.25 * x - 3*x*x + x*x*x

!  Versuchen Sie auch die andere Testfunktion:
!  fun = 1 - x - Sin(x)/Cos(x) 

!  Ein Problemfall (Nullstelle nahe 1.32, aber versuchen Sie a=1., b=1.55):
!  fun = x - Sin(x)/Cos(x) / 3

END FUNCTION fun

