!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  PROGRAM newton
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  (c) 2005, Christian B. Lang, Mathematische Methoden in der Physik
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  Eine Nullstelle einer gegebenen Funktion wird mit dem Newtownschen
!  Verfahren gesucht. Als Startwerte muessen die Schranken eines
!  Intervalls angegeben werden, in dem die Funktion eine Nullstelle
!  hat. Wenn das Verfahren aus dem Intervall fuehrt,
!  erfolgt eine Fehlermeldung und Neustart.
!  Geeignete Testwerte: a= 0,    b=2   f. 1.Testfunktion
!                       a= -0.1, b=0.3 f. 2.Testfunktion
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

PROGRAM newton

   implicit none
   
!  Genauigkeit wird vorgegeben
   real, parameter :: feps = 0.000001
   
   real::x,a,b,fx,fa,fb,fsx,fsa,fsb
   integer::iter
   
   
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
   call fun(a, fa, fsa)
   write(*,*) 'f (',a,') = ',fa
   x = a
   if(abs(fa) < feps) goto 20
   write(*,*) 'Obere  Intervallgrenze b = '
   read(*,*) b
   call fun(b, fb, fsb)
   write(*,*) 'f (',b,') = ',fb
   x = b
   IF(abs(fb) < feps) goto 20
   x = (a + b)/2
   write(*,*) 'Weiter mit RETURN '
   read(*,*) 

   write(*,*) 'Iteration : '
   write(*,*) '        x : '
   
   do while(iter<100)
      iter = iter + 1
      call fun(x, fx, fsx)
      write(*,*) 'Iteration          = ',iter
      write(*,*) '        x : ',x, ' f(x) : ', fx
!  Damit man sieht, was passiert kann entfernt werden:
      write(*,*) 'Weiter mit RETURN '
      read(*,*)
      if(abs(fx)<feps) goto 20
      if(fsx.ne.0) x = x - fx/fsx
!  Wenn die Ableitung verschwindet, divergiert das Verfahren}
      if((x < a).or.(x > b).or.(fsx == 0)) then
        write(*,*) 'Newton Verfahren fuehrt aus dem Intervall!'
        write(*,*) 'Weiter mit RETURN '
        goto 10
      end if
   end do
  
   write(*,*) 'Iteration konvergiert nicht: iter = ',iter

20 write(*,*) 'ERGEBNIS:'
   write(*,*) 'Iteration          = ', iter
   write(*,*) 'Bester Wert von x  = ',x
   write(*,*) '              f(x) = ',fx
   
END PROGRAM Newton

SUBROUTINE fun (x,fx,fsx)
! Berechnet Funktion fx, Ableitung fsx

   implicit none

   real, intent(in)::x
   real, intent(out)::fx,fsx

!  1. Testfunktion :
   fx = 0.75 - 0.25*x - 3*x*x + x*x*x
   fsx= -0.25 - 6*x + 3*x*x

!  Versuchen Sie auch die 2. Testfunktion:
!  fx =  1 - x - sin(x)/cos(x)
!  fsx= -1 -  /Sqr(cos(x))  

END SUBROUTINE fun



