L'algoritmo della bisezione

Fortran di Giuseppe Ciaburro guida dal 04-06-2002

Popolarità:

economia news e media viaggi informatica internet salute e benessere int rattenimento e spettacolo sport tempo libero istruzio ne e formazione arte cultura scienza

Le categorie della guida

I link preferiti della guida

Il Sondaggio

Quale versione del fortran utilizzi?

Guarda i risultati

Algoritmi

A cura di Giuseppe Ciaburro

Pubblicato il 24/08/2001

» Invia tramite EMAIL

» Versione per la STAMPA

» Le vostre opinioni

Algoritmo della bisezione per il calcolo degli zeri di una funzione.

Il criterio di ricerca degli zeri più semplice consiste nel continuare a dividere in due l'intervallo di ricerca iniziato, e continuare a scegliere l'intervallo agli estremi del quale la funzione cambia segno.
Non è un criterio molto efficiente, ma sicuramente converge a qualcosa, anche se la funzione è discontinua o ha punti di singolarità. Se nell'intervallo ci sono più zeri, convergerà ad uno di questi. Se c'è una discontinuità a gradino, convergerà alla discontinuità, e se c'è un punto singolare convergerà al punto singolare.

Supponiamo che

sia una funzione continua nell'intervallo

tale che

cioè la funzione assume valori di segno opposto nei due estremi. Allora per il teorema dell'esistenza degli zeri $\exists c : a < c < b, f(c)=0$ . Definiamo $x_{1}=\frac{a+b}{2}$ allora si può avere uno dei seguenti tre casi:

$f(x_{1})=0$ allora abbiamo avuto fortuna: $x_{1}$ è una radice e l'algoritmo termina.
$f(a)f(x_{1})<0$ allora possiamo reiterare il procedimento sul nuovo intervallo $[a,x_{1}]$ .
$f(b)f(x_{1})<0$ allora possiamo reiterare il procedimento sul nuovo intervallo $[x_{1},b]$ .

In un passo abbiamo in questo modo dimezzato l'ampiezza dell'intervallo di confidenza, quindi reiterando questo metodo è possibile trovare un intervallo di confidenza di ampiezza sempre più piccola. Poiché $x_{1}=\frac{a+b}{2}$ si ha $\vert\overline{x}-x_{1}\vert<\frac{b-a}{2}$ . In generale $\vert\overline{x}-x_{i}\vert<\frac{b_{i}-a_{i}}{2}$ . Sviluppando questa relazione si ottiene

$\displaystyle \vert\overline{x}-x_{i}\vert<\frac{b-a}{2^{i}} \leq \epsilon$
$\displaystyle b-a \leq \epsilon 2^{i}$
$\displaystyle \frac{b-a}{\epsilon} \leq 2^{i}$
$\displaystyle \log_{2}{(b-a)}-\log_{2}{(\epsilon)} \leq i$
$\displaystyle i \geq \lceil \log_{2}{(b-a)}-\log_{2}{(\epsilon)} \rceil = \nu$			(1.1)

In questo modo, fissata la accuratezza $\epsilon$ , viene fissato anche il numero massimo di iterazioni, $\nu$ . Nella nostra implementazione, epsilon viene dato dall'utente e $\nu$ viene calcolato secondo la formula (1.1) e restituito in nu. Un altro problema che si presenta nell'implementazione dell'algoritmo è la condizione $f(x_{i})=0$ che è difficilmente verificabile in aritmetica finita e dev'essere più realisticamente sostituita con una condizione del tipo $\vert f(x_{i})\vert<\delta$ dove $\delta$ soddisfa

$\displaystyle \delta=\delta(\epsilon) \textrm{ tale che } \vert f(x_i)\vert<\delta \Longrightarrow \vert\overline{x}-x_{i}\vert<\epsilon$

(1.2)

Tale $\delta$ risulta essere

$\displaystyle \delta \approx \epsilon \vert f'(x_{i})\vert$

(1.3)

Tuttavia calcolare esplicitamente $f'(x_{i})$ sarebbe più costoso del metodo di bisezione stesso. Quello che possiamo fare è approssimare $f'(x_{i})$ con il rapporto incrementale che ha come estremi l'intervallo corrente $[a_{i},b_{i}]$ :

$\displaystyle \delta \approx \epsilon \Big\vert\frac{f(b_{i})-f(a_{i})}{b_{i}-a_{i}}\Big\vert$

(1.4)

Di seguito è riportato il codice in Fortran 90 che impelmenta l'algoritmo della bisezione per il calcolo degli zeri di una funzione.


      program main
      !programma in fortran 90 per il calcolo
      !degli zeri di una funzione con il metodo
      !di bisezione
      integer, parameter ::m=20
      real :: fa = 0.0, fb =1.0, f 
      real :: ga = 0.5, gb =2.0, g 

      interface
      function f(x)
      real, intent(in) :: x
      end function 
      function g(x)
      real, intent(in) :: x
      end function g 
      end interface

      print *, "first function f"
      print *, "a =",fa, "b =",fb
      call bisect(f,fa,fb,m)
      print *, "second function g"
      print *, "a =", ga, "b =",gb
      call bisect(g,ga,gb,m)
      end program main

      function f(x)  
      real, intent (in) :: x
      f = ((x)*x - 3.0)*x + 1.0     
      end function f 

      function g(x)
      real, intent (in) :: x
      g = x**3 - 2.0*sin(x)
      end function
      
      subroutine bisect(f,a,b,m)
      integer, intent (in) :: m 
      real :: a, b, c, fa, fb, fc, error
      interface  
      function f(x) 
      real, intent(in) :: x  
      end function f          
      end interface

      fa = f(a)             
      fb = f(b)
      if ( sign(1.0,fa) == sign(1.0,fb) ) then
        print *,"function has same sign at",a,b  
      else       
        print *,"n    c    f(c)     error"  
        error = b - a
        do n = 0,m        
            error = error/2.0
            c  = a + error
            fc = f(c)
            if ( sign(1.0,fa) /= sign(1.0,fc) ) then     
               b = c                  
               fb = fc
            else  
               a = c    
               fa = fc  
            endif  
            print *, n, c, fc, error   
         end do    
      end if       
      end subroutine bisect

Feed