Il problema è quello di fornire una stima approssimata dell'integrale definito di una funzione reale di una variabile reale, e una stima dell'errore commesso nel passaggio dal continuo al discreto.

Descrizione dell'algoritmo

L'algoritmo implementa un procedimento di quadratura automatica adattativa, che consta di un modulo di quadratura locale (LQM), in questo caso il modulo Cavalieri-Simpson, di una strategia per selezionare il sottointervallo da suddividere ad ogni passo, e di un criterio di arresto.

La strategia adottata per la scelta dell'intervallo è di tipo locale, suddivide l'intervallo di lavoro solo nel caso in cui l'errore ad esso relativo derivante dall'applicazione dell'LQM è maggiore della tolleranza accettata per quell'intervallo (in base alla sua ampiezza); per contrastare poi la tendenza naturale dell'algoritmo ad infittire il calcolo nel caso in cui la tolleranza inserita sia molto piccola, si preferisce richiedere in ingresso anche una misura dell'ampiezza minima alla quale si può ridurre un intervallo di integrazione nonché un numero massimo di valutazioni eseguibili.

La gestione dei sottointervalli di integrazione è eseguita per mezzo di una pila che ha nel top l'intervallo di ampiezza minore.

Ogni qualvolta un intervallo è scaricato dalla pila, cioè quando l'errore dell'LQM è minore della tolleranza relativa al sottointervallo di calcolo, è aggiunto al valore dell'integrale Q ed a quello dell'errore E il contributo locale calcolato.

L'algoritmo termina o per lo svuotamento della pila o per il raggiungimento del numero massimo di valutazioni consentite.

Raccomandazioni sull'uso

La tolleranza inserita deve essere non negativa. Il numero massimo di valutazioni eseguibili deve essere maggiore di cinque. Il passo minimo di integrazione deve essere positivo e minore dell'ampiezza dell'intervallo di integrazione

f - funzione integranda
a - reale, estremo inferiore d'integrazione
b - reale, estremo superiore d'integrazione
tol - reale, tolleranza assoluta
hmin - reale, peso minimo d'integrazione
nmax - intero, massimo numero di valutazioni

output

q - reale, valore approssimato dell'integrale
e - reale, stima dell'errore
nval - intero, numero di valutazioni occorse
ifail - intero, indicatore d'errore

      Program quadd
	   
c     Programma dimostrativo della subroutine quad
      
      external quad
      
      real f1,f2,f3,f4,f5
      external f1,f2,f3,f4,f5 
      
      integer dim
      parameter (dim=25)

      real a,b,tol,hmin,q,e,x(dim),y(dim)
      integer nmax,nval,ifail
      character risp
      character*12 fin,fout
      logical warn
      
      print*,'*** PROGRAMMA DIMOSTRATIVO DELLA SUBROUTINE QUAD ***'
      print*
   10 continue
      print*,'Funzione integranda :'
      print*
      print*,'1) sin(x)'
      print*
      print*,'2) (100/x**2)*sin(10/x)'
      print*
      print*,'3) 10*(x-1)*(x-0.5)*(x-0.25)*(x-0.125)*(x-0.0625)',
    %               '*(x-.03125)'
      print*
      print*,'4) 1/sqrt(x)'
      print*
      print*,'5) e^(-x^2)'
      print*
      print*,'Numero della funzione= '
      read(*,*)n
      print*,'Integrale definito/indefinito [d/i]= '
      read(*,'(A)')risp                                        
      if (risp .eq. 'd') then
         print*,'Estremi d''integrazione= '
         read(*,*)a,b         
         num=1
      endif
      if (risp .eq. 'i') then
         print*,'Punto iniziale= '
         read(*,*)a
         print*,'Valore iniziale= '
         read(*,*)Fa
         print*,'Numero di punti di valutazione= '
         read(*,*)num
         print*,'Nome del file di input = '
         read(*,'(A)')fin
     
         open(11,FILE=fin,STATUS='OLD')
         do 15 i=1,num
            read(11,*)x(i)
   15    continue
         close(11)
      endif
     
      print*,'Tolleranza richiesta= '
      read(*,*)tol
      print*,'Passo minimo d''integrazione= '
      read(*,*)hmin
      print*,'Numero massimo di valutazioni (minimo 5)= '
      read(*,*)nmax
      warn=.false.
      do 20 i=1,num
         if (risp .eq. 'i') b=x(i)
         if (n .eq. 1)
    %      call quad(f1,a,b,tol,hmin,nmax,q,e,nval,ifail)
         if (n .eq. 2) 
    %      call quad(f2,a,b,tol,hmin,nmax,q,e,nval,ifail)
         if (n .eq. 3)
    %      call quad(f3,a,b,tol,hmin,nmax,q,e,nval,ifail)
         if (n .eq. 4)
    %      call quad(f4,a,b,tol,hmin,nmax,q,e,nval,ifail)
         if (n .eq. 5)
    %      call quad(f5,a,b,tol,hmin,nmax,q,e,nval,ifail)
            
         if (ifail .eq. 2)
    %      print*,'Errore! tolleranza assoluta non valida'
         if (ifail .eq. 16)
    %      print*,'Errore! numero di valutazioni non valido'
         if (ifail .eq. 17)
    %      print*,'Errore! passo minimo d''integrazione non valido'
         if (ifail .eq. 18)
    %      print*,'Errore! stack overflow'
         if (ifail .ge. 2 .and. ifail .le. 18) goto 30
         if (ifail .eq. 19) warn=.true.
            
         if (risp .eq. 'i') y(i)=Fa+q
   20 continue
         
      if (warn) then
        print*,'Attenzione! raggiunto il numero massimo ',
    %         'di valutazioni'
        print*
      endif
      if (risp .eq. 'd') then
         write(*,40)'Stima dell''integrale= ',q
         print*
         write(*,40)'Stima dell''errore= ',e
         print*
         write(*,50)'Valutazioni effettuate= ',nval
      endif
      if (risp .eq. 'i') then
         print*,'Nome del file di output='
         read(*,'(A)')fout
         open(12,FILE=fout)
         write(12,60)(y(i),i=1,num)
         close(12)     
         print*,'Valutazione approssimata dell''integrale ',
    %          'indefinito nei punti assegnati:'
         write(*,60)(y(i),i=1,num)
      endif
   30 print*
      print*,'c per continuare, un altro carattere per terminare:'
      read(*,'(A)')risp
      if (risp .eq. 'c') goto 10
      
   40 format(1X,A,E14.7)
   50 format(1X,A,I6)                
   60 format(1X,E14.7)
   
      End