CAPITOLO 1 Esercizi risolti di programmazione Fortranberga/Teaching/CalcoloNumerico/fortran.pdf ·...

CAPITOLO 1

Esercizi risolti di programmazione Fortran

In questa appendice si propongono testi e risoluzioni di alcuni temi d’esame di pro-grammazione apparsi negli ultimi 6 anni nelle prove scritte per Ingegneria IndustrialeCivile ed Ambientale. La soluzione e scritta nel vecchio linguaggio Fortan 77 ma puoessere compilata anche da piu moderni compilatori quali il Fortran90 o il Fortran95.I programmi che seguono rappresentano naturalmente una delle possibili soluzioni aitemi d’esame e quindi hanno lo scopo di suggerire una traccia allo studente.

Per ogni tema si forniscono: l’enunciato, evenutali commenti e il codice Fortran disoluzione diviso in due parti: un main che contiene il programma principale e tutti isottoprogrammi richiesti.

Prima di cominciare, e opportuno richiamare brevemente la sequenza di operazioniche e necessario compiere per potere eseguire un programma Fortran sotto il sistemaoperativo Linux, che e quello adottato per esempio nell’aula di Laboratorio per i corsidi studio in Ingegneria.

Per una trattazione approfondita e dettagliata dell’uso del Fortran per ApplicazioniNumeriche si consiglia il testo numero [?] della bibliografia.

1

2 Esercizi risolti di programmazione Fortran

1.1 Scrittura compilazione ed esecuzione di un pro-gramma Fortran

I passi da seguire per l’esecuzione di un programma Fortran possono essere schematiz-zati come segue.

1. Creare un file in cui scrivere il programma vero e proprio (programma sorgente).Tale file deve essere salvato con una estensione esplicita .f ad esempio radici.f.Usare per questa operazione un editore, cioe un programma che permette diaprire, modificare salvare un file. L’editore piu semplice e kedit che si puo farepartire semplicemente digitando kedit sulla linea di comandi.

2. Compilare il programma sorgente per la generazione del programma eseguibile.Si deve dare il seguente comando.

g77 file sorgente −o file eseguibile

Il nome del file eseguibile e arbitrario e non necessita di alcuna estensione. Adesempio

g77 radici.f −o radici

Questo comando invoca il compilatore dandogli come file di input il file ra-dici.f che abbiamo precedentemente generato e restituisce come output il fileeseguibile radici.

3. Esecuzione del programma. A questo punto per fare eseguire il programma sideve semplicemente dare l’istruzione

./file eseguibile

Nel nostro esempio

./radici

1.1#1 3

Programma sorgente per il calcolo delle soluzioni di un’equazione di se-condo grado.

radici.fprogram rad i c i

c! ca l co la le so luz ion i r e a l i di una equazione di secondo gradoc

imp l i c i t none

real a , b , c , delta , x1 , x2

write (6 ,∗ ) ’ i n s e r i s c i i tre c o e f f i c i e n t i de l la1 equazione di secondo grado ’

read (5 ,∗ ) a , b , c

delta = b∗∗2 − 4.∗a∗c

i f ( delta . l t . 0 . 0 ) thenwrite (6 ,∗ ) ’ non c i sono so luz ion i r e a l i ’

e l se i f ( delta . eq . 0 . 0 ) thenx1 = −b / ( 2 . 0∗ a )write (6 ,∗ ) ’ so luz ion i co inc ident i in x1 = ’ , x1

e lsex1 = (−b + sqrt ( delta ) ) / ( 2 . ∗ a )x2 = (−b − sqrt ( delta ) ) / ( 2 . ∗ a )write (6 ,∗ ) ’ l e due so luz ion i sono : x1 = ’ , x1 , ’ x2 = ’ , x2

end i fstopend

Esempio di esecuzione del programma.

1. Dopo aver digitato dalla linea dei comandi: ./radici

2. apparira sullo schermo la scritta inserisci i tre coefficienti della equazionedi secondo grado

3. a questo punto digitare i numeri: 1 -3 2 (separati da uno spazio o dal tasto diinvio).

4. Sullo schermo appariranno le soluzioni dell’equazione di II grado x2 − 3x+ 2 = 0

x1 =1 x2 = 2


Esempio di programma sorgente per l’implementazione del metodo delpunto fisso per la soluzione dell’equazione x = cos(x)

punto fisso.fprogram punto f issoimp l i c i t none

! dichiarazione de l l e v a r i a b i l iinteger i ter , imaxreal x0 , to l , xnew , scarto , xold

! apertura f i l e di outputopen (23 , f i l e = ’ r i s u l t a t i ’ )

! l e t tura dei dati inputwrite (6 ,∗ ) ’ i n s e r i s c i x0 t o l imax ’read (5 ,∗ ) x0 , to l , imax

! i n i z i a l i z z a z i o n e de l l e v a r i a b i l ixold = x0i t e r = 0scarto = 2. d0∗ t o lwrite (32 ,900)

c! c i c l o i t e r a t i v o per la implementazione del metodo di punto f i s s oc

do while ( scarto . gt . t o l . and . i t e r . l t . imax )i t e r = i t e r + 1xnew = cos ( xold )scarto = abs (xnew − xold )

! importante : stampa r i s u l t a t i intermedi! secondo un formato d e f i n i t o dal la riga 1000

write (32 ,1000) i ter , xnew , scartoxold = xnew

end do1000 format ( i8 , f21 .14 , d15 . 5 )

i f ( scarto . l e . t o l ) then! convergenza ! !

write (32 ,∗ ) ’ soluzione f ina l e ’ ,xnewwrite (32 ,∗ ) ’ i t e r a z i o n i ’ , i t e rwrite (32 ,∗ ) ’ scarto f ina l e ’ , scarto

e lse! non convergenza ! !

write (32 ,∗ ) ’ soluzione non trovata ’end i fc lose (32)

900 format (6x , ” i t ” ,14x , ”xk” ,12x , ” scarto ” )stopend

1.2#1 5

• compilazione:

g77 punto_fisso.f -o punto_fisso

• esecuzione:

./punto_fisso

• dati in ingresso:

1 1.e-4 100

• file in output: risultati

it xk scarto1 0.54030227661133 0.45970D+002 0.85755324363708 0.31725D+003 0.65428978204727 0.20326D+004 0.79348033666611 0.13919D+005 0.70136880874634 0.92112D-016 0.76395964622498 0.62591D-017 0.72210246324539 0.41857D-018 0.75041770935059 0.28315D-019 0.73140406608582 0.19014D-0110 0.74423736333847 0.12833D-0111 0.73560476303101 0.86326D-0212 0.74142509698868 0.58203D-0213 0.73750686645508 0.39182D-0214 0.74014735221863 0.26405D-0215 0.73836916685104 0.17782D-0216 0.73956722021103 0.11981D-0217 0.73876029253006 0.80693D-0318 0.73930388689041 0.54359D-0319 0.73893773555756 0.36615D-0320 0.73918443918228 0.24670D-0321 0.73901826143265 0.16618D-0322 0.73913019895554 0.11194D-0323 0.73905479907990 0.75400D-04

soluzione finale 0.7390548iterazioni 23scarto finale 7.5399876E-05


1.2 Temi d’esame risoltiEsercizio 1.1 Date le matrici quadrate A e B di dimensione N (N ≤ 30) si vuolecalcolare la matrice C = AB +BA, la sua traccia e la sua norma di Frobenius.

Scrivere un programma in linguaggio Fortran che:

(a) legge N,A e B;

(b) calcola la matrice C = AB +BA impiegando la subroutine PRODMAT che esegueil prodotto tra due matrici quadrate;

(c) calcola la traccia α della matrice C usando la function TRAC;

(d) calcola la norma di Frobenius β della matrice C(β =

√∑ni,j=1 c

2ij

)servendosi

della function EUCMAT;

(e) stampa con commento α e β.

Risoluzione.

Si e definita una costante nmax per dimensionare al massimo le matrici. Tale parametrodeve essere passato, unitamente alla dimensione effettiva delle matrici stesse, a tutti isottoprogrammi che coinvolgono matrici.

Si noti che e sufficiente scrivere una sola subroutine: PRODMAT per eseguire ilcalcolo della matrice C. Essa deve essere chiamata due volte con parametri invertiti (siveda il main alla pagina seguente):

call PRODMAT(...,A,B,...)call PRODMAT(...,B,A,...)

1.2#1 7

Listing 1.1: noneprogram B1impl i c i t none

! dichiarazione de l l e v a r i a b i l iinteger nmaxparameter (nmax=30)integer i , j , k , nreal∗8 a (nmax,nmax) ,b (nmax,nmax) , c (nmax,nmax)real∗8 d (nmax,nmax) , e (nmax,nmax)real∗8 alfa , beta , trac , eucmat

! le t tura dei datiopen (15 , f i l e = ’ dati ’ )read (15 ,∗ ) ndo i = 1 ,n

read (15 ,∗ ) ( a ( i , j ) , j =1 ,n )end dodo i = 1 ,n

read (15 ,∗ ) ( b ( i , j ) , j =1 ,n )end do

! chiamata subroutine PRODMAT! D = AB

c a l l PRODMAT(n ,nmax, a , b , d )! E = BA

c a l l PRODMAT(n ,nmax, b , a , e )! C = AB + BA

do i = 1 ,ndo j = 1 ,n

c ( i , j ) = d ( i , j ) + e ( i , j )end do

end doal fa = TRAC(n ,nmax,C)beta = EUCMAT(n ,nmax,C)

!write (6 ,∗ ) ’ t racc ia di C ’ , al fa , ’norma di Frobenius ’ , betac lose (15)

!stopend


sottoprogrammi–ES1! sottoprogrammi

subroutine PRODMAT(n ,nmax, a , b , c )imp l i c i t noneinteger i , j , k , n ,nmaxreal∗8 a (nmax,nmax) ,b (nmax,nmax) , c (nmax,nmax)

do i = 1 ,ndo j = 1 ,n

c ( i , j ) = 0 . d0do k = 1 ,n

c ( i , j ) = c ( i , j ) + a ( i , k ) ∗b (k , j )end do

end doend do

returnend

!real∗8 function TRAC(n ,nmax, c )imp l i c i t noneinteger i , n ,nmaxreal∗8 c (nmax,nmax)

trac = 0. d0do i = 1 ,n

trac = trac + c ( i , i )end do

returnend

!real∗8 function EUCMAT(n ,nmax, c )imp l i c i t noneinteger i , j , n ,nmaxreal∗8 c (nmax,nmax)

eucmat = 0. d0do i = 1 ,n

do j = 1 ,neucmat = eucmat + c ( i , j ) ∗∗2

end doend doeucmat = sqrt ( eucmat )

returnend

1.2#1 9

Esercizio 1.2 Si vuole risolvere il sistema lineare Ax = b di dimensione m ≤ 80mediante il metodo iterativo di Richardson:

xk+1 = (I − βA)xk + βb

con β un numero reale positivo. Scrivere un programma in linguaggio FORTRAN che:

1. legge la soluzione iniziale x0, il parametro β, la tolleranza sullo scarto TOLL e ilnumero massimo di iterazioni ITMAX;

2. ad ogni iterazione calcola il vettore z = βAxk mediante la subroutine BMATVET,la soluzione corrente xk+1 = xk − z + βb e lo scarto dk+1 = xk+1 − xk;

3. calcola la norma euclidea del vettore scarto ‖dk+1‖2 mediante la function EUC;

4. arresta le iterazioni quando ‖dk+1‖2 ≤ TOLL oppure il numero di iterazioni supe-ra ITMAX;

5. stampa la soluzione ed il numero di iterazioni effettuate.

Risoluzione.

Il metodo iterativo si implementa utilizzando due vettori xold e xnew per codificarerispettivamente xk e xk+1. Il ciclo iterativo e tradotto in Fortran utilizzando il costrut-to while: si continua ad iterare fintantoche sono contemporaneamente verificate lecondizioni: ‖dk+1‖2 > TOLL e iter < ITMAX.


B2 main.fprogram B2impl i c i t none

! dichiarazione de l l e v a r i a b i l iinteger nmaxparameter (nmax=80)real∗8 xold (nmax) ,xnew(nmax) ,A(nmax,nmax) , z (nmax) ,dk (nmax)real∗8 b (nmax)real∗8 t o l l , beta , scarto , eucinteger n , itmax , i ter , i , j

! l e t tura dati da f i l eopen (15 , f i l e = ’ dati ’ )open (16 , f i l e = ’ output ’ )read (15 ,∗ ) n , itmax , t o l l , betado i = 1 ,n

read (15 ,∗ ) ( a ( i , j ) , j =1 ,n )end doread (15 ,∗ ) ( xold ( i ) , i =1 ,n )read (15 ,∗ ) ( b ( i ) , i =1 ,n )

! i n i z i a l i z z a z i o n ei t e r = 0scarto = 1. d10

! c i c l o i t e r a t i v odo while ( i t e r . l t . itmax . and . scarto . gt . t o l l )

c a l l bmatvet (n ,nmax, beta , xold ,A, z )do i = 1 ,n

xnew( i ) = xold ( i ) − z ( i ) + beta∗b ( i )dk ( i ) = xnew( i ) − xold ( i )

end doscarto = euc (n , dk )do i = 1 ,n

xold ( i ) = xnew( i )end do

end doi f ( scarto . l e . t o l l ) then

write (16 ,∗ ) ’ soluzione raggiunta in ’ , i ter , ’ i t e r a z i o n i ’do i = 1 ,n

write (16 ,∗ ) xnew( i )end do

elsewrite (6 ,∗ ) ’ soluzione non raggiunta in ’ , itmax , ’ i t e r a z i o n i ’

end i fc lose (15)c lose (16)stopend

1.2#1 11

sottoprogrammi–ES2subroutine bmatvet (n ,nmax, beta ,A, x , z )imp l i c i t noneinteger n ,nmaxreal∗8 x (n) , z (n ) ,a (nmax,nmax) , betainteger i , j

!do i = 1 ,n

z ( i ) = 0 . d0do j = 1 ,n

z ( i ) = z ( i ) + a ( i , j ) ∗x ( j )end doz ( i ) = beta∗z ( i )

end do

returnend

!!!

real∗8 function EUC(n , x )imp l i c i t noneinteger nreal∗8 x (n)integer i

!EUC = 0. d0do i = 1 ,n

EUC = EUC + x ( i ) ∗∗2end doEUC = sqrt (EUC)

returnend


Esercizio 1.3 Date la matrice quadrata A ed il vettore x di dimensione N ≤ 45 si vuolecalcolare il vettore y = A200x secondo il seguente algoritmo:

1. y0 = x

2. yk+1 = Ayk, k = 0, · · · 199

3. y = y200

in cui il punto 2. deve essere implementato richiamando la subroutine MATVET

Scrivere un programma in linguaggio FORTRAN che:

1. legge N, A e x .

2. calcola il vettore y impiegando l’algoritmo precedente.

3. stampa con commento il vettore y.

Si fornisca inoltre la subroutine MATVET.

Risoluzione.

sottoprogrammi–ES3subroutine matvet (n ,nmax, x ,A, z )imp l i c i t noneinteger n ,nmaxreal∗8 x (n) , z (n ) ,a (nmax,nmax)integer i , j

cdo i = 1 ,n

z ( i ) = 0 . d0do j = 1 ,n

z ( i ) = z ( i ) + a ( i , j ) ∗x ( j )end do

end do

returnend

1.2#1 13


! dichiarazione de l l e v a r i a b i l iinteger nmaxparameter (nmax=45)real∗8 yold (nmax) ,ynew(nmax) ,A(nmax,nmax) , x (nmax)real∗8 t o l l , scartointeger n , i , j

! l e t tura dati da f i l eopen (15 , f i l e = ’ dati ’ )open (16 , f i l e = ’ output ’ )read (15 ,∗ ) ndo i = 1 ,n

read (15 ,∗ ) ( a ( i , j ) , j =1 ,n )end doread (15 ,∗ ) ( x ( i ) , i =1 ,n )do i = 1 ,n

yold ( i ) = x ( i )end do

cdo j = 1 ,200

c a l l matvet (n ,nmax, yold ,A, ynew)do i = 1 ,n

yold ( i ) = ynew( i )end do

end do

write (16 ,∗ ) ’ vettore y ’do i = 1 ,n

write (16 ,∗ ) ynew( i )end doc lose (15)c lose (16)

stopend


Esercizio 1.4 Si scriva un programma Fortran che calcola il quoziente di Rayleigh

q =xTAx

xTx

di una matrice quadrata di dimensione al massimo 40 × 40. Si richede in partico-lare di implementare una subroutine che esegue il prodotto di una matrice quadrataper un vettore (MATVET) e una function che esegue il prodotto scalare tra due vettori(SCAL) Il programma principale dovra leggere la matrice ed il vettore da file, chiamareopportunamente i sottoprogrammi e infine stampare con commento il risultato q.

Risoluzione.

Per quanto riguarda la subroutine MATVET si faccia riferimento al precedente eserci-zio B.3. Di seguito scriviamo l’implementazione della function SCAL.

sottoprogrammi–ES4real∗8 function scal (n , x , y )imp l i c i t none

integer nreal∗8 x (n) , y (n )integer i

sca l = 0 . d0do i = 1 ,n

sca l = sca l + x ( i ) ∗y ( i )end do

returnend

1.2#1 15


! dichiarazione de l l e v a r i a b i l iinteger nmaxparameter (nmax=45)real∗8 A(nmax,nmax) , x (nmax) , z (nmax)real∗8 q , scal , numeratore , denominatoreinteger n , i , j

! l e t tura dati da f i l eopen (15 , f i l e = ’ dati ’ )read (15 ,∗ ) ndo i = 1 ,n

read (15 ,∗ ) ( a ( i , j ) , j =1 ,n )end doread (15 ,∗ ) ( x ( i ) , i =1 ,n )

!! z = Ax!

c a l l matvet (n ,nmax, x ,A, z )!! numeratore = x ˆT z!

numeratore = scal (n , x , z )!! denominatore = x ˆT x!

denominatore = scal (n , x , x )q = numeratore / denominatore

! stampa a videowrite (6 ,∗ ) ’ i l quoziente di Rayleigh vale ’ , qc lose (15)

stopend


Esercizio 1.5 Si vuole risolvere il sistema lineare Ax = b di dimensione m ≤ 40mediante il seguente metodo iterativo:

xk+1 = Exk + q

con E matrice quadrata m ×m e q vettore noto. Scrivere un programma in linguaggioFORTRAN che:

1. legge la soluzione iniziale x0, la tolleranza sullo scarto TOLL e il numero massimodi iterazioni ITMAX; legge inoltre il vettore iniziale x0, la matrice E ed il vettoreq.

2. ad ogni iterazione calcola la soluzione corrente xk+1 utilizzando la subroutineMATVET per il prodoto matrice per vettore e lo scarto dk+1 = xk+1 − xk;

3. calcola la norma euclidea del vettore scarto ‖dk+1‖2 mediante la function EUC;

4. arresta le iterazioni quando ‖dk+1‖2 ≤ TOLL oppure il numero di iterazionisupera ITMAX;

5. stampa il numero di iterazioni effettuate e lo scarto finale.

Risoluzione.

I sottoprogrammi richiesti sono gia stati descritti in esercizi precedenti. In particolareper la function EUC si veda l’esercizio B.2 mentre per la subroutine MATVET l’esercizioB.3.

1.2#1 17


! dichiarazione de l l e v a r i a b i l iinteger nmaxparameter (nmax=40)real∗8 xold (nmax) ,xnew(nmax) ,A(nmax,nmax) ,dk (nmax) , z (nmax)real∗8 q (nmax)real∗8 t o l l , scarto , eucinteger n , itmax , i ter , i , j

! l e t tura dati da f i l eopen (15 , f i l e = ’ dati ’ )open (16 , f i l e = ’ output ’ )read (15 ,∗ ) n , itmax , t o l ldo i = 1 ,n

read (15 ,∗ ) ( a ( i , j ) , j =1 ,n )end doread (15 ,∗ ) ( xold ( i ) , i =1 ,n )read (15 ,∗ ) ( q ( i ) , i =1 ,n )

! i n i z i a l i z z a z i o n ei t e r = 0scarto = 1. d10

! c i c l o i t e r a t i v odo while ( i t e r . l t . itmax . and . scarto . gt . t o l l )

i t e r = i t e r + 1c a l l matvet (n ,nmax, xold ,A, z )do i = 1 ,n

xnew( i ) = z ( i ) + q ( i )dk ( i ) = xnew( i ) − xold ( i )

end doscarto = euc (n , dk )do i = 1 ,n

xold ( i ) = xnew( i )end do

end doi f ( scarto . l e . t o l l ) then

write (16 ,∗ ) ’ soluzione raggiunta in ’ , i ter , ’ i t e r a z i o n i ’write (16 ,∗ ) ’ norma de l l o scarto f ina l e : ’ , scarto

e lsewrite (6 ,∗ ) ’ soluzione non raggiunta in ’ , itmax , ’ i t e r a z i o n i ’

end i fc lose (15)c lose (16)stopend


Esercizio 1.6 Si scriva un programma Fortran che implementa il seguente metodoiterativo

xn+1 = xn −f(xn)

K, n = 0, 1, 2, . . .

dove x0 e il punto iniziale, f(x) = x3 − cos(x) e K e una costante assegnata. Il testdi arresto da utilizzare e basato sullo scarto e sul numero massimo di iterazioni. Ilprogramma dovra leggere da file K, x0, nmax e toll e visualizzare al termine del ciclo:la soluzione trovata, lo scarto finale ed il numero di iterazioni impiegato. Si implementicome function la funzione f .

Risoluzione.

B6 main.fprogram B6impl i c i t noneinteger j ,nmax, i t e rreal∗8 f , k , xold , xnew , x0 , scarto , t o l l

read (5 ,∗ ) x0 , t o l l ,nmax, kxold = x0scarto = 1. d0i t e r = 0do while ( abs ( scarto ) . gt . t o l l . and . i t e r . l t .nmax)

i t e r = i t e r + 1xnew = xold − f ( xold ) / kscarto = xnew − xoldxold = xnew

end doi f ( abs ( scarto ) . l e . t o l l ) then

write (6 ,∗ ) ’ soluzione = ’ ,xnewwrite (6 ,∗ ) ’ i t e r a z i o n i = ’ , i t e rwrite (6 ,∗ ) ’ scarto f ina l e= ’ , scarto

e lsewrite (6 ,∗ ) ’ soluzione non trovata in ’ ,nmax, ’ i t e r a z i o n i ’

end i fstopend

!real∗8 function f ( x )imp l i c i t nonereal∗8 xf = x∗∗3 − cos ( x )returnend

1.2#1 19

Esercizio 1.7 Si scriva un programma Fortran che approssimi l’integrale

I =

∫ 2

−1

(e−x + x3)dx

utilizzando la formula di Cavalieri Simpson composta. In particolare il programmadeve

1. Leggere da tastiera gli estremi di integrazione e il numero m di applicazioni dellaformula.

2. implementare tramite function: CAVSIM la formula semplice di CavalieriSimpson.

3. Implementare la formula composta chiamando opportunamente CAVSIM,

4. implementare tramite due function la funzione integranda f e la sua primitiva.

5. Visualizzare a video il numero di intervalli, l’approssimazione dell’integrale, ilvalore vero e l’errore vero.

Risoluzione.

sottoprogrammi–ES7real∗8 function CAVSIM( a , b )imp l i c i t nonereal∗8 a , b , fun , xmedxmed = .5 d0∗ ( a+b )cavsim = ( b−a ) ∗ ( fun ( a ) +4.d0∗fun (xmed)+fun ( b ) ) / 6 . d0returnend

creal∗8 function primitiva ( x )imp l i c i t nonereal∗8 xprimitiva = −exp(−x ) + x∗∗4∗0.25returnend

creal∗8 function fun ( x )imp l i c i t nonereal∗8 xfun = exp(−x ) + x∗∗3returnend


B7 main.fprogram B7

!imp l i c i t nonereal∗8 lung , valore vero , a , b , fun , primitivareal∗8 csim ,CAVSIM, errore , xa , xbinteger i ,m

!read (5 ,∗ ) a , b ,mvalore vero = primitiva ( b ) − primitiva ( a )write (6 , ’ ( a ,2 f8 . 2 ) ’ ) ’ i n t e rva l l o di integrazione ’ ,a , bwrite (6 , ’ ( a , f12 . 6 ) ’ ) ’ valore vero ’ , valore verocsim = 0. d0lung = ( b−a ) / d f l oa t (m)do i = 0 ,m−1

xa = a + df l oat ( i ) ∗ lungxb = xa + lungcsim = csim + CAVSIM( xa , xb )

end do!

errore = abs ( valore vero − csim )write (6 ,∗ ) ’ numero di i n t e r v a l l i ’ ,mwrite (6 ,∗ ) ’ integrale approssimato ’ , csimwrite (6 ,∗ ) ’ errore assoluto ’ , errore

!stopend

1.2#1 21

Esercizio 1.8 Data la matrice quadrataA di dimensione n ≤ 50 e i fattori triangolari Le U , basso e alto rispettivamente, tali che LU = A con A un’approssimazione sparsificatadi A, si vuole misurare la lunghezza euclidea del vettore differenza (v − x) ove x e unvettore assegnato e:

v = L−1AU−1x

Scrivere un programma in linguaggio FORTRAN che:

1. legge n, A, L, U e x;

2. calcola il vettore y = U−1x mediante sostituzioni all’indietro;

3. calcola il prodotto z = Ay con la subroutine MATVET;

4. calcola il vettore v = L−1z mediante sostituzioni in avanti;

5. calcola il vettore differenza w = v − x;

6. calcola la norma euclidea di w mediante la function EUC e la stampa.

Suggerimento: si ricordi che per le componenti i-esime dei vettore y e v valgono le re-lazioni yi =

(xi −

∑nj=i+1 uijyj

)/uii e vi =

(zi −

∑i−1j=1 lijzj

)/lii. Inoltre, le yi vanno

calcolate dall’ultima alla prima, mentre le vi dalla prima all’ultima.

Risoluzione.

Prestare attenzione alla lettura delle matrici L ed U .

Per quanto riguarda la function EUC si veda l’implementazione dell Esercizio B.1mentre la subroutine MATVET e gia presente nell implementazione dell’EsercizioB.3.



! dichiarazione de l l e v a r i a b i l iinteger nmaxparameter (nmax=50)real∗8 A(nmax,nmax) ,L(nmax,nmax) ,U(nmax,nmax)real∗8 v (nmax) , x (nmax) , z (nmax) , y (nmax) ,w(nmax)real∗8 euc ,nrminteger n , i , j

! l e t tura dati da f i l eopen (15 , f i l e = ’ dati ’ )read (15 ,∗ ) n

! le t tura matrice Ado i = 1 ,n

read (15 ,∗ ) ( a ( i , j ) , j =1 ,n )end do

! le t tura matrice Ldo i = 1 ,n

read (15 ,∗ ) ( l ( i , j ) , j =1 , i )end do

! le t tura matrice Udo i = 1 ,n

read (15 ,∗ ) (u ( i , j ) , j =i , n )end do

! le t tura vettore xread (15 ,∗ ) ( x ( i ) , i =1 ,n )

! backward subsitut iondo i = n,1,−1

y ( i ) = x ( i )do j = i +1 ,n

y ( i ) = y ( i ) − U( i , j )∗y ( j )end doy ( i ) = y ( i ) / u ( i , i )

end do! z = A y

c a l l MATVET(n ,nmax, y ,A, z )! forward subsitut ion

do i = 1 ,nv ( i ) = z ( i )do j = 1 , i−1

v ( i ) = v ( i ) − L( i , j )∗v ( j )end dov ( i ) = v ( i ) / l ( i , i )

end dodo i = 1 ,n

w( i ) = v ( i ) − x ( i )end donrm = EUC(n ,w)write (6 ,∗ ) ’ norma di w= ’ ,nrmclose (15)

stopend

1.2#1 23

Esercizio 1.9Si vuole risolvere l’equazione non lineare f(x) = 0, con f(x) =

√x2 + 1 + lnx, mediante

il seguente schema iterativo noto come schema di Steffensen:

xk+1 = xk −h2k

f(xk + hk)− hkin cui hk = f(xk). Scrivere un programma in linguaggio FORTRAN che:

1. legge la soluzione iniziale x0, la tolleranza sullo scarto τ e il numero massimo diiterazioni ITMAX;

2. implementa lo schema assegnato arrestando le iterazioni quando |xk+1 − xk| ≤ τoppure il numero di iterazioni supera ITMAX;

3. stampa la soluzione ed il numero di iterazioni effettuate.

Risoluzione.

B9 main.fprogram B9impl i c i t noneinteger j , itmax , i t e rreal∗8 f , xold , xnew , x0 , scarto , tau , hk

read (5 ,∗ ) x0 , tau , itmaxxold = x0scarto = 1. d0i t e r = 0do while ( abs ( scarto ) . gt . tau . and . i t e r . l t . itmax )

i t e r = i t e r + 1hk = f ( xold )xnew = xold − hk∗∗2 / ( f ( xold+hk)−hk)write (6 ,∗ ) i ter , xnew , scartoscarto = xnew − xoldxold = xnew

end doi f ( abs ( scarto ) . l e . tau ) then

write (6 ,∗ ) ’ soluzione = ’ ,xnewwrite (6 ,∗ ) ’ i t e r a z i o n i = ’ , i t e rwrite (6 ,∗ ) ’ scarto f ina l e= ’ , scarto

e lsewrite (6 ,∗ ) ’ soluzione non trovata in ’ , itmax , ’ i t e r a z i o n i ’

end i fstopend

!real∗8 function f ( x )imp l i c i t nonereal∗8 xf = sqrt ( x∗∗2+1.d0 ) + log ( x )returnend

CAPITOLO 1 Esercizi risolti di programmazione Fortranberga/Teaching/CalcoloNumerico/fortran.pdf ·...

Documents

Transcript of CAPITOLO 1 Esercizi risolti di programmazione Fortranberga/Teaching/CalcoloNumerico/fortran.pdf ·...