Calcolo numerico

Programmazione in Fortran 77.

Integrazione numerica

 In questa rubrica presento una completa serie di algoritmi implementati in Digital Visual Fortran 90 ( totalmente compatibile con il FORTRAN 77 ).

Nel 2005 questi programmi sono serviti agli studenti dei corsi serali per studenti lavoratori per superare l'esame universitario "Calcolo numerico e programmazione" presso la facoltą di ingegneria di Padova ( corso della professoressa Michela Redivo Zaglia per allievi aereospaziali).

I programmi sono testati e funzionanti ma essi non assicurano il superamento della prova d'esame senza un approfondito studio teorico nel testo "Calcolo numerico"  dell'autrice Redivo Zaglia.

Attenzione: Nell'anno in cui questi programmi sono stati creati sono stati prodotti in sede d'esame (quindi vanno ben memorizzati nella loro essenza e non solo a memoria) usando il compilatore in linea del sistema operativo Linux installato in aula Taliercio della facoltą di ingegneria in via Belzoni  7 (PD). E' indispensabile al fine di superare la prova d'esame, sapere compilare in ambiente Linux ed effettuare la visualizzazione degli eventuali grafici con i tools di default presenti in quel sistema operativo.

I programmi  name.f  posso essere tagliati ed incollati nell'editor di programmazione per un corretto funzionamento.   Ove sia necessario bisognerą far trovare nel direttorio di compilazione i file di testo contenenti l'insieme di input "formattato" ovvero con righe e colonne al posto in cui l'algoritmo le aspetta.  Nella stessa directory troveremo il file testuale di output.

Marco Gottardo.

Metodo di Cavalieri Simpson
cavsimpson.f cav.out
C====================================================================
C Programma per il calcolo integrale (Metodo CAVALIERI SIMPSN)
C====================================================================
program CavSimpson
implicit none
C====================
C Parte dichiarativa
C====================
real*8 funz
real*8 a,b,int
integer m
real*8 x(0:100)
integer n
character*20 nomefile
C====================
C Parte esecutiva
C====================

C====================
C Ingresso dati
C====================
print *, 'Dammi gli estremi a e b dell''intervallo iniziale'
read*, a, b
b=dacos(b)

print *, 'Dammi il numero di suddivisioni per integrare m'
read*, m
do while (mod(m,2).ne.0)
print *, 'Reinserire m pari '
read*, m
enddo

print *, 'Dammi il nome del file in uscita'
read*, nomefile

open(unit=10,file=nomefile)

C============================================
C Chiamata della formula CAVALIERI SIMPSON
C============================================

Call cavsimp(m,b,a,int)

C====================
C Uscita dati
C====================

print*, 'Valore approssimato integrale int =', int
write(10,*)'**********************************'
write(10,*)'* Gottardo Marco matr 347349/IT             *'
write(10,*)'* esame di calcolo numerico                        *'
write(10,*)'* 03/07/2003                                              *'
write(10,*)'* Metodo di integrazione Cavalieri Simpson *'
write(10,*)'* Comparazione dei risultati                          *'
write(10,*)'* funz=funz=exp(x)*dcos(x)                        *'
write(10,*)'* vedi tabella risultati pag 208                      *'
write(10,*)'* tra 0 e pi cav(512) -12.0703463160        *'
write(10,*)'*                                                                 *'
write(10,*)'***********************************'
write(10,*)
write(10,*)'Estremo iniziale a = ',a
write(10,*)'Estremo finale b = ',b
write(10,*)'Suddivisioni m = ',m
write(10,*)'Valore approssimato dell''integrale = ',int
close(10)

stop
end

C====================================================================
C Algoritmo del calcolo integrale con metodo CAVALIERI SIMPSON
C Applicazione dell'integrazione nell'intervallo assegnato
C====================================================================
Subroutine cavsimp(m,b,a,int)
implicit none
real*8 h,a,b,x,funz,int,fx
integer m,i
h=(b-a)/m
int=funz(a)+funz(b)
do i= 1, m-1, 2
x=a+(i*h)
Int=Int+(4.0d0*funz(x))
enddo
do i=2, m-2, 2
x=a+(i*h)
Int= Int+(2.0d0*funz(x))
enddo
Int= h*(Int/3.0d0)
return
end

C=======================
C Funzione assegnata
C=======================
real*8 function funz(x)
implicit none
real*8 x
funz= exp(x)*dcos(x)
return
end		  ********************************
* Gottardo Marco matr 347349/IT          *
* esame di calcolo numerico                     *
* 03/07/2003                                           *
* Metodi di integrazione numerica             *
* Comparazione dei risultati                      *
* funz=funz=exp(x)*dcos(x)                     *
* vedi tabella risultati pag 208                   *
* tra 0 e pi cav(512) -12.0703463160     *
* tra 0 e pi trap(512) -12.0704220570     *
*                                                              *
*********************************

Estremo iniziale a = 0.000000000000000E+000
Estremo finale b = 3.14159265358979

numero di divisioni sull'intervallo = 512

Valore approssimato dell'integrale
Metodo dei trapezi -> I= -12.2128471815589
Metodo Cavalieri Simpson -> I= -12.0703463163890
 

 

Programma di calcolo dell'integrale di una funzione predefinita, in un intervallo assegnato con il metodo dei trapezi.

Metodo dei trapezi.
trapezi.f trapezi.out
C====================================================================
C Programma per il calcolo integrale (Metodo Trapezi)
C====================================================================
program trapezi
implicit none
C====================
C Parte dichiarativa
C====================
real*8 funz
real*8 a,b,int
integer m
real*8 x(0:100)
integer n
character*20 nomefile
C====================
C Parte esecutiva
C====================

C====================
C Ingresso dati
C====================
print *, 'Dammi gli estremi a e b dell''intervallo iniziale'
read*, a, b
b=dacos(b)

print *, 'Dammi il numero di suddivisioni per integrare m'
read*, m
do while (mod(m,2).ne.0)
print *, 'Reinserire m pari '
read*, m
enddo

print *, 'Dammi il nome del file in uscita'
read*, nomefile

open(unit=10,file=nomefile)

C============================================
C Chiamata della formula Trapezi
C============================================

Call trapez(m,b,a,int)

C====================
C Uscita dati
C====================

print*, 'Valore approssimato integrale int =', int
write(10,*)'*********************************'
write(10,*)'* Gottardo Marco matr 347349/IT            *'
write(10,*)'* esame di calcolo numerico                       *'
write(10,*)'* 03/07/2003                                             *'
write(10,*)'* Metodo di integrazione numerica Trapezi *'
write(10,*)'* Comparazione dei risultati                       *'
write(10,*)'* funz=funz=exp(x)*dcos(x)                      *'
write(10,*)'* vedi tabella risultati pag 208                   *'
write(10,*)'* tra 0 e pi trap(512) -12.0704220570    *'
write(10,*)'*                                                              *'
write(10,*)'*********************************'
write(10,*)
write(10,*)'Estremo iniziale a = ',a
write(10,*)'Estremo finale b = ',b
write(10,*)'Suddivisioni m = ',m
write(10,*)'Valore approssimato dell''integrale = ',int
close(10)

stop
end

C====================================================================
C Algoritmo del calcolo integrale con metodo Trapezi
C Applicazione dell'integrazione nell'intervallo assegnato
C====================================================================
Subroutine trapez(m,b,a,int)
implicit none
real*8 h,a,b,x,funz,int
integer m,i
h=(b-a)/m
int=funz(a)+funz(b)
do i=1, m-1
x=a+(i*h)
Int= Int+(2.0d0*funz(x))
enddo
Int= h*(Int/2.0d0)
return
end

C=======================
C Funzione assegnata
C=======================
real*8 function funz(x)
implicit none
real*8 x
funz= exp(x)*dcos(x)
return
end		  ********************************
* Gottardo Marco matr 347349/IT           *
* esame di calcolo numerico                      *
* 03/07/2003                                            *
* Metodo di integrazione numerica Trapezi *
* Comparazione dei risultati                        *
* funz=funz=exp(x)*dcos(x)                       *
* vedi tabella risultati pag 208                     *
* tra 0 e pi trap(512) -12.0704220570       *
*                                                                 *
***********************************

Estremo iniziale a = 0.000000000000000E+000
Estremo finale b = 3.14159265358979
Suddivisioni m = 512
Valore approssimato dell'integrale = -12.0704220570084
 

Programma comparativo dei due metodi di integrazione.

Programma comparativo dei trapezi/Cavalieri S.
integrali.f integrali.out
C====================================================================
C Gottardo Marco matr 347349/IT
C Esame di calcolo numerico
C Sessione estiva 3/07/2003
C Programma per il calcolo integrale (comparazione tra i metodi)
C====================================================================
program integrale
implicit none
C====================
C Parte dichiarativa
C====================
real*8 funz,trapez,cav_sim,C_S
real*8 a,b,integr,integ,integS,CS,caval
integer divis
real*8 x(0:100)
integer n
character*20 nomefile
C====================
C Parte esecutiva
C====================
C Ingresso ed uscita dati
C
print *, 'Dammi gli estremi a e b dell''intervallo iniziale'
read*, a, b
print *, 'Dammi il numero di suddivisioni per integrare'
read*, divis
c
print *, 'Dammi il nome del file in uscita'
read*, nomefile
c
integr = trapez(divis,b,a,integ)
caval = C_S(a,b,divis,integS)

open(unit=10,file=nomefile)
write(10,*)'************************************'
write(10,*)'* Gottardo Marco matr 347349/IT                 *'
write(10,*)'* esame di calcolo numerico                           *'
write(10,*)'* 03/07/2003                                                 *'
write(10,*)'* Metodi di integrazione numerica                   *'
write(10,*)'* Comparazione dei risultati                             *'
write(10,*)'*funz=(1.0d0/(1.0d0+(x-(acos(-1.0d0)))**2))*'
write(10,*)'*                                                                     *'
write(10,*)'* tra 0 e 5 I=2.33976628367 circa                 *'
write(10,*)'************************************'
write(10,*)
write(10,*)'Estremo iniziale a = ',a
write(10,*)'Estremo finale b = ',b
write(10,*)
write(10,*)'numero di divisioni sull''intervallo = ',divis/4
write(10,*)
write(10,*)'Valore approssimato dell''integrale '
write(10,*)'Metodo dei trapezi -> I=',integ
write(10,*)'Metodo Cavalieri Simpson -> I=',integS
close(10)
stop
end

C====================================================================
C Funzione assegnata
C====================================================================
real*8 function funz(x)
implicit none
real*8 x
funz=(1.0d0/(1.0d0+(x-(acos(-1.0d0)))**2))
return
end
C====================================================================
C Algoritmo del calcolo integrale con metodo dei Trapezi
C Applicazione dell'integrazione nell'intervallo assegnato
C====================================================================
real*8 function trapez(divis,b,a,integ)
implicit none
real*8 h,a,b,sum,xi,funz,integ
integer divis,index
h=(b-a)/divis
print *,'Ampiezza base trapezi h ',h
do while (index.le.divis)
xi=a+(index*h)
sum=(((funz(xi)+funz(xi+h))*h)/2)
integ=sum+integ
index=index+1
end do
write(6,*)
print *, ' Integrale calcolato'
write(*,'(a,f20.15)') 'Soluzione appross. metodo trapezi',integ
return
end
C====================================================================
C Algoritmo del calcolo integrale con metodo di Cavalieri Simpson
C Applicazione dell'integrazione nell'intervallo assegnato
C====================================================================
real*8 function C_S(a,b,n,integS)
implicit none
integer i,n,iter,n1,n2
real*8 a,b,x,In,I2n,h,S1,Xo,X1,X2,INT(2),funz,integS
iter =0
30 iter=iter+1
h=(b-a)/(2.0d0*n)
S1=0.0d0
X2=a
do 40 i=1,n
Xo=X2
X1=Xo+h
X2=X1+h
S1=S1+h*(funz(xo)+4.0d0*funz(x1)+funz(X2))/3.0d0
40 continue
INT(iter)=S1
n=2*n
if(iter.lt.2) goto 30
IN=INT(1)
I2n=INT(2)
n1=n/4
n2=n/2
c write(6,*)n1,In,n2,I2n
integS=I2n
write(6,*)'Soluzione approssimata cavalieri Simpson:',I2n
return
end		  ************************************
* Gottardo Marco matr 347349/IT                 *
* esame di calcolo numerico                            *
* 03/07/2003                                                  *
* Metodi di integrazione numerica                    *
* Comparazione dei risultati                             *
*funz=(1.0d0/(1.0d0+(x-(acos(-1.0d0)))**2))*
*                                                                     *
* tra 0 e 5 I=2.33976628367 circa                 *
************************************

Estremo iniziale a = 0.000000000000000E+000
Estremo finale b = 5.00000000000000

numero di divisioni sull'intervallo = 500

Valore approssimato dell'integrale
Metodo dei trapezi -> I= 2.34200029777255
Metodo Cavalieri Simpson -> I= 2.33976628366840