Algoritmi

L'algoritmo di Newton-Raphson

Algoritmo di Newton-Raphson per il calcolo degli zeri di una funzione.

Esistono vari metodi per accelerare il processo di convergenza nella ricerca degli zeri. Anche in questo caso, prendiamo subito in considerazione il metodo "migliore" (se non altro quello più utilizzato). Si tratta del metodo di Newton, detto anche di Newton-Raphson, e si basa sull'assunzione che la derivata prima della funzione in esame possa essere calcolata. Questa la ricetta:

• Calcolare la derivata in un punto dell'intervallo di ricerca, il che permette di ottenere la retta tangente alla funzione in quel punto.
• Calcolare quindi il punto in cui la tangente passa per lo zero: se tale punto esce dall'intervallo di ricerca bisogna gettare la spugna.
• Calcolare di nuovo la derivata in quel punto, e ripetere dal punto (2) fino a convergenza avvenuta.
  
  

Nelle immediate vicinanze del punto di zero, il metodo funziona molto bene (la sua motivazione nasce infatti dallo sviluppo di Taylor, arrestato in questo caso al prim'ordine). Allontanandosi dal punto di zero, la funzione potrebbe avere peculiarità che rendono il metodo assai problematico, come in questo (pessimo) esempio, dove le tangenti indicate in rosso sono solo causa di guai:

Un buon approccio di programmazione ci permette però di sommare la rapidità del metodo di Newton-Raphson con la stabilità del metodo di bisezione: basta ricorrere ad un passo di bisezione quando il metodo di Newton tenderebbe ad uscire dall'intervallo oppure converge troppo lentamente.

Per comprendere meglio come implementare l'algoritmo in oggetto riporto il codice relativo in fortran 90.

 
      program main
      integer, parameter :: dp = kind (1d0)
      real (kind = dp) :: x, true
      integer :: m      
      
      interface
      function f(x)
      integer, parameter :: dp = kind(1d0)
      real (kind=dp), intent(in) :: x
      end function f
      function g(x)
      integer, parameter :: dp = kind (1d0)
      real (kind=dp), intent(in) :: x
      end function g
      end interface

      m = 6
      x  =  4.0_dp 
      call newton(f,g,x,m)
      end program main

      subroutine newton(f,g,x,m)
      interface
      function f(x)
      integer, parameter :: dp = kind(1d0)
      real (kind=dp), intent(in) :: x
      end function f
      function g(x)
      integer, parameter :: dp = kind(1d0)
      real (kind =dp), intent(in) :: x
      end function g
      end interface
      integer, parameter :: dp = kind(1d0)
      real (kind=dp), intent(in) :: x
      integer, intent(in) :: m
      
      print *, "n         x                f(x)"
      fx = f(x)   
      print *, 0, x, fx
      do n = 1,m 
           x = x - fx/g(x)     
           fx = f(x) 
           print "(i2, f20.16, f20.16)", n, x, fx      
      end do 
      end subroutine newton

      function f(x)
      integer, parameter :: dp = kind (1d0)
      real (kind = dp) , intent(in) :: x
      f = ((x - 2.0_dp)*x + 1.0_dp)*x - 3.0_dp
      end function f

      function g(x)
      integer, parameter  :: dp = kind (1d0)
      real (kind = dp), intent(in) :: x
      g = (3.0_dp*x - 4.0_dp)*x + 1.0_dp
      end function g      

• » Invia tramite EMAIL
• » Versione per la STAMPA
• » Le vostre opinioni