Concorso Energia in gioco IT

Calcolo numerico

Programmazione in Fortran 77.

Equazioni differenziali.

Valgono le seguenti regole:

I file nome.f contengono l'implementazione dell'argoritmo di calcolo.
I file nome.out contengono l'output dell'argoritmo di calcolo e sono creati dal programma.
I file con estensione .dsp .dsw sono creati dal compilatore e non sono qui riportati perchè inutili per la teoria.
il caratter "C" davanti una riga indica Commento ovvero codice che non viene compilato.

Metodo di Eulero
eulero.f	eulero.out
C==================================================================== C Programma per il Metodo di EULERO esplicito C Versione con controlli in ingresso C==================================================================== program Eul implicit none C==================== C Parte dichiarativa C==================== real8 funz real8 funzs, sol external funz real8 b , h integer maxdim parameter (maxdim=500) real8 x(0:maxdim), y(0:maxdim) integer m, i character20 nomefile C==================== C Parte esecutiva C==================== C Ingresso dati C C inserire 0 1 C b = 0.2 C h =0.02 C la tabella è a pagina 235 per verifica dei risultati print , 'Dammi i valori x0 ed y0 ' read, x(0), y(0) print , 'Dammi l''estremo b dell''intervallo' read, b print , 'Dammi il passo h ' read, h do while ((b.ge.x(0) .and. h.le.0.0d0) .or. (b.le.x(0) .and. h.ge.0.0d0)) print , 'Valori non corretti Reinserire x0, y0, b e h' read, x(0), y(0), b, h end do print , 'Dammi il nome del file in uscita' read, nomefile open(unit=10,file=nomefile) C C Esecuzione metodo di Eulero esplicito C CALL euleroe (funz, b, h, maxdim, x, y, m) C C Uscita finale risultati C write(10,) 'Valori x0 = ',x(0),' y0 = ',y(0), ' b = ', b write(10,) 'Numero di intervalli ', m write(10,) 'Passo ', h write(10,) write(10,) 'Vettori delle discretizzazioni ed approssimazioni' do i = 0, m write(10,'(i5,f10.5,g25.15)') i, x(i), y(i) end do C C Calcolo e stampa dell'errore relativo C sol = funzs(x(m)) write(10,) write(10,'(T1,A,T13,A,T33,A,T50,A)') * 'Estremo b', 'Valore finale', 'Valore esatto', 'Errore relativo' write(10,'(f8.4,2g20.10,e12.5)') x(m), y(m), sol, * dabs((y(m)-sol)/sol) close(10) stop end C==================================================================== C Subroutine per metodo di Eulero Esplicito C==================================================================== SUBROUTINE euleroe (f, b, h, maxdim, x, y, m) implicit none C C Variabili INGRESSO: C f funzione y'=f(x,y) C b estremo finale dell'intervallo C h passo di integrazione C maxdim numero massimo di intervalli (da non superare) C C Variabili INGRESSO/USCITA: C x vettore contenente in ingresso in x(0) il valore x0 C in uscita contiene i valori della discretizzazione C (lunghezza 0:m) C y vettore contenente in ingresso in y(0) il valore y0 C in uscita contiene i valori della approssimazioni C (lunghezza 0:m) C C Variabili USCITA: C m numero degli intervalli (m <= maxdim) C C==================== C Parte dichiarativa C==================== real8 f external f real8 b, h integer maxdim, m, n real8 x(0:), y(0:) C==================== C Parte esecutiva C==================== m = idnint((b - x(0))/h) if (m .gt. maxdim) then print , * ' m >= maxdim. Viene preso m=maxdim e ricalcolato il passo h' m = maxdim h = (b - x(0))/m print , 'Nuovo passo h = ', h end if do n = 1, m y(n) = y(n-1) + h f(x(n-1),y(n-1)) x(n) = x(0) + n * h end do return end C==================================================================== C Funzione assegnata C==================================================================== real8 function funz(x,y) implicit none real8 x, y funz = 1.0d0+x+y*2 return end C==================================================================== C Funzione soluzione C==================================================================== real8 function funzs(x) implicit none real*8 x funzs = 1.53268118d0 return end	Valori x0 = 0.000000000000000E+000 y0 = 1.00000000000000 b = 0.200000000000000 Numero di intervalli 10 Passo 2.000000000000000E-002 Vettori delle discretizzazioni ed approssimazioni 0 0.00000 1.00000000000000 1 0.02000 1.04000000000000 2 0.04000 1.08203200000000 3 0.06000 1.12624786498048 4 0.08000 1.17281655004794 5 0.10000 1.22192652324927 6 0.12000 1.27378861181367 7 0.14000 1.32863936036539 8 0.16000 1.38674501136364 9 0.18000 1.44840624589448 10 0.20000 1.51396385895740 Estremo b Valore finale Valore esatto Errore relativo 0.2000 1.513963859 1.532681180 0.12212E-01

Un algoritmo più efficente di quello di Eulero è il metodo di Eulero Modificato di cui riporto sotto l'implementazione.

Metodo di Eulero modificato
euleromod.f	euleromod.out
C==================================================================== C Programma per il Metodo di Eulero modificato C Versione con controlli in ingresso C==================================================================== program Eulero modificato implicit none C==================== C Parte dichiarativa C==================== real8 funz real8 funzs, sol external funz real8 b , h integer maxdim parameter (maxdim=500) real8 x(0:maxdim), y(0:maxdim) integer m, i character20 nomefile C==================== C Parte esecutiva C==================== C Ingresso dati C print , 'Dammi i valori x0 ed y0 ' read, x(0), y(0) print , 'Dammi l''estremo b dell''intervallo' read, b print , 'Dammi il passo h ' read, h do while ((b.ge.x(0) .and. h.le.0.0d0) .or. (b.le.x(0) .and. h.ge.0.0d0)) print , 'Valori non corretti Reinserire x0, y0, b e h' read, x(0), y(0), b, h end do print , 'Dammi il nome del file in uscita' read, nomefile open(unit=10,file=nomefile) C C Esecuzione metodo di Eulero modificato C CALL Eulmod(funz, b, h, maxdim, x, y, m) C C Uscita finale risultati C write(10,) 'Valori x0 = ',x(0),' y0 = ',y(0), ' b = ', b write(10,) 'Numero di intervalli ', m write(10,) 'Passo ', h write(10,) write(10,) 'Vettori delle discretizzazioni ed approssimazioni' do i = 0, m write(10,'(i5,f10.5,g25.15)') i, x(i), y(i) end do C C Calcolo e stampa dell'errore relativo C sol = funzs(x(m)) write(10,) write(10,'(T1,A,T13,A,T33,A,T50,A)') * 'Estremo b', 'Valore finale', 'Valore esatto', 'Errore relativo' write(10,'(f8.4,2g20.10,e12.5)') x(m), y(m), sol, * dabs((y(m)-sol)/sol) close(10) stop end C==================================================================== C Subroutine per metodo di Eulero modificato C==================================================================== SUBROUTINE Eulmod (f, b, h, maxdim, x, y, m) implicit none C C Variabili INGRESSO: C f funzione y'=f(x,y) C b estremo finale dell'intervallo C h passo di integrazione C maxdim numero massimo di intervalli (da non superare) C C Variabili INGRESSO/USCITA: C x vettore contenente in ingresso in x(0) il valore x0 C in uscita contiene i valori della discretizzazione C (lunghezza 0:m) C y vettore contenente in ingresso in y(0) il valore y0 C in uscita contiene i valori della approssimazioni C (lunghezza 0:m) C C Variabili USCITA: C m numero degli intervalli (m <= maxdim) C C==================== C Parte dichiarativa C==================== real8 f external f real8 b, h,k1,k2 integer maxdim, m, n real8 x(0:), y(0:) C==================== C Parte esecutiva C==================== c calcolo dell'intero più vicino con la funzione predefinita m = idnint((b - x(0))/h) if (m .gt. maxdim) then print , * ' m >= maxdim. Viene preso m=maxdim e ricalcolato il passo h' m = maxdim h = (b - x(0))/m print , 'Nuovo passo h = ', h end if do n = 1, m k1=f(x(n-1),y(n-1)) k2=f(x(n-1)+(h/2.0d0),y(n-1)+((h/2.0d0)k1)) y(n) = y(n-1) + (h * k2) x(n) = x(0) + n * h end do return end C==================================================================== C Funzione assegnata C==================================================================== real8 function funz(x,y) implicit none real8 x, y funz = 1.0d0+x+y*2 return end C==================================================================== C Funzione soluzione C==================================================================== real8 function funzs(x) implicit none real*8 x funzs = 1.53268118d0 return end	Valori x0 = 0.000000000000000E+000 y0 = 1.00000000000000 b = 0.200000000000000 Numero di intervalli 10 Passo 2.000000000000000E-002 Vettori delle discretizzazioni ed approssimazioni 0 0.00000 1.00000000000000 1 0.02000 1.04100800000000 2 0.04000 1.08416679064475 3 0.06000 1.12964571228276 4 0.08000 1.17763420159369 5 0.10000 1.22834482400350 6 0.12000 1.28201687376056 7 0.14000 1.33892066960675 8 0.16000 1.39936270799904 9 0.18000 1.46369188037965 10 0.20000 1.53230701979822 Estremo b Valore finale Valore esatto Errore relativo 0.2000 1.532307020 1.532681180 0.24412E-03

Algoritmo di risoluzione delle equazioni differenziali con il metodo di Heun.

Metodo di Heun
heun.f	heun.out
C==================================================================== program Heun implicit none C==================== C Parte dichiarativa C==================== real8 funz real8 funzs, sol external funz real8 b , h integer maxdim parameter (maxdim=500) real8 x(0:maxdim), y(0:maxdim) integer m, i character20 nomefile C==================== C Parte esecutiva C==================== C Ingresso dati C print , 'Dammi i valori x0 ed y0 ' read, x(0), y(0) print , 'Dammi l''estremo b dell''intervallo' read, b print , 'Dammi il passo h ' read, h do while ((b.ge.x(0) .and. h.le.0.0d0) .or. (b.le.x(0) .and. h.ge.0.0d0)) print , 'Valori non corretti Reinserire x0, y0, b e h' read, x(0), y(0), b, h end do print , 'Dammi il nome del file in uscita' read, nomefile open(unit=10,file=nomefile) C C Esecuzione metodo di Heun C CALL eun (funz, b, h, maxdim, x, y, m) C C Uscita finale risultati C write(10,) 'Valori x0 = ',x(0),' y0 = ',y(0), ' b = ', b write(10,) 'Numero di intervalli ', m write(10,) 'Passo ', h write(10,) write(10,) 'Vettori delle discretizzazioni ed approssimazioni' do i = 0, m write(10,'(i5,f10.5,g25.15)') i, x(i), y(i) end do C C Calcolo e stampa dell'errore relativo C sol = funzs(x(m)) write(10,) write(10,'(T1,A,T13,A,T33,A,T50,A)') * 'Estremo b', 'Valore finale', 'Valore esatto', 'Errore relativo' write(10,'(f8.4,2g20.10,e12.5)') x(m), y(m), sol, * dabs((y(m)-sol)/sol) close(10) stop end C==================================================================== C Subroutine per metodo di Heun C==================================================================== SUBROUTINE eun (f, b, h, maxdim, x, y, m) implicit none C C Variabili INGRESSO: C f funzione y'=f(x,y) C b estremo finale dell'intervallo C h passo di integrazione C maxdim numero massimo di intervalli (da non superare) C C Variabili INGRESSO/USCITA: C x vettore contenente in ingresso in x(0) il valore x0 C in uscita contiene i valori della discretizzazione C (lunghezza 0:m) C y vettore contenente in ingresso in y(0) il valore y0 C in uscita contiene i valori della approssimazioni C (lunghezza 0:m) C C Variabili USCITA: C m numero degli intervalli (m <= maxdim) C C==================== C Parte dichiarativa C==================== real8 f external f real8 b, h,k1,k2 integer maxdim, m, n real8 x(0:), y(0:) C==================== C Parte esecutiva C==================== m = idnint((b - x(0))/h) if (m .gt. maxdim) then print , * ' m >= maxdim. Viene preso m=maxdim e ricalcolato il passo h' m = maxdim h = (b - x(0))/m print , 'Nuovo passo h = ', h end if do n = 1, m k1=f(x(n-1),y(n-1)) x(n) = x(0) + n h k2=f(x(n),y(n-1)+(hk1)) y(n) = y(n-1) + h ((k1+k2)/2) end do return end C==================================================================== C Funzione assegnata C==================================================================== real8 function funz(x,y) implicit none real8 x, y funz = 1.0d0+x+y*2 return end C==================================================================== C Funzione soluzione C==================================================================== real8 function funzs(x) implicit none real*8 x funzs = 1.53268118d0 return end	Valori x0 = 0.000000000000000E+000 y0 = 1.00000000000000 b = 0.200000000000000 Numero di intervalli 10 Passo 2.000000000000000E-002 Vettori delle discretizzazioni ed approssimazioni 0 0.00000 1.00000000000000 1 0.02000 1.04101600000000 2 0.04000 1.08418398880415 3 0.06000 1.12967350449933 4 0.08000 1.17767421987840 5 0.10000 1.22839898356911 6 0.12000 1.28208743112089 7 0.14000 1.33901029467996 8 0.16000 1.39947457414558 9 0.18000 1.46382977757830 10 0.20000 1.53247549790791 Estremo b Valore finale Valore esatto Errore relativo 0.2000 1.532475498 1.532681180 0.13420E-03

Metodo di Runnge Kutta.

Metodo di Runge Kutta Rango 4
rungek4.f	rk4.out
C==================================================================== C Programma per il Metodo di Runge Kutta rango 4 C Versione con controlli in ingresso C==================================================================== program Runge Kutta rango 4 implicit none C==================== C Parte dichiarativa C==================== real8 funz real8 funzs, sol external funz real8 b , h integer maxdim parameter (maxdim=500) real8 x(0:maxdim), y(0:maxdim) integer m, i character20 nomefile C==================== C Parte esecutiva C==================== C Ingresso dati C print , 'Dammi i valori x0 ed y0 ' read, x(0), y(0) print , 'Dammi l''estremo b dell''intervallo' read, b print , 'Dammi il passo h ' read, h do while ((b.ge.x(0) .and. h.le.0.0d0) .or. (b.le.x(0) .and. h.ge.0.0d0)) print , 'Valori non corretti Reinserire x0, y0, b e h' read, x(0), y(0), b, h end do print , 'Dammi il nome del file in uscita' read, nomefile open(unit=10,file=nomefile) C C Esecuzione metodo di Runge Kutta rango 4 C CALL rungeK4(funz, b, h, maxdim, x, y, m) C C Uscita finale risultati C write(10,) 'Valori x0 = ',x(0),' y0 = ',y(0), ' b = ', b write(10,) 'Numero di intervalli ', m write(10,) 'Passo ', h write(10,) write(10,) 'Vettori delle discretizzazioni ed approssimazioni' do i = 0, m write(10,'(i5,f10.5,g25.15)') i, x(i), y(i) end do C C Calcolo e stampa dell'errore relativo C sol = funzs(x(m)) write(10,) write(10,'(T1,A,T13,A,T33,A,T50,A)') * 'Estremo b', 'Valore finale', 'Valore esatto', 'Errore relativo' write(10,'(f8.4,2g20.10,e12.5)') x(m), y(m), sol, * dabs((y(m)-sol)/sol) close(10) stop end C==================================================================== C Subroutine per metodo di Runge Kutta rango 4 C==================================================================== SUBROUTINE rungeK4 (f, b, h, maxdim, x, y, m) implicit none C C Variabili INGRESSO: C f funzione y'=f(x,y) C b estremo finale dell'intervallo C h passo di integrazione C maxdim numero massimo di intervalli (da non superare) C C Variabili INGRESSO/USCITA: C x vettore contenente in ingresso in x(0) il valore x0 C in uscita contiene i valori della discretizzazione C (lunghezza 0:m) C y vettore contenente in ingresso in y(0) il valore y0 C in uscita contiene i valori della approssimazioni C (lunghezza 0:m) C C Variabili USCITA: C m numero degli intervalli (m <= maxdim) C C==================== C Parte dichiarativa C==================== real8 f external f real8 b, h,k1,k2,k3,k4 integer maxdim, m, n real8 x(0:), y(0:) C==================== C Parte esecutiva C==================== c calcolo dell'intero più vicino con la funzione predefinita m = idnint((b - x(0))/h) if (m .gt. maxdim) then print , * ' m >= maxdim. Viene preso m=maxdim e ricalcolato il passo h' m = maxdim h = (b - x(0))/m print , 'Nuovo passo h = ', h end if do n = 1, m k1=f(x(n-1),y(n-1)) k2=f(x(n-1)+(h/2.0d0),y(n-1)+((h/2.0d0)k1)) k3=f(x(n-1)+(h/2.0d0),y(n-1)+((h/2.0d0)k2)) x(n) = x(0) + (n h) k4=f(x(n),y(n-1)+(hk3)) y(n) = y(n-1) + ((h /6.0d0)(k1+2.0d0k2+2.0d0k3+k4)) end do return end C==================================================================== C Funzione assegnata C==================================================================== real8 function funz(x,y) implicit none real8 x, y funz = 1.0d0+x+y*2 return end C==================================================================== C Funzione soluzione C==================================================================== real8 function funzs(x) implicit none real*8 x funzs = 1.53268118d0 return end	Valori x0 = 0.000000000000000E+000 y0 = 1.00000000000000 b = 0.200000000000000 Numero di intervalli 10 Passo 2.000000000000000E-002 Vettori delle discretizzazioni ed approssimazioni 0 0.00000 1.00000000000000 1 0.02000 1.04102465769336 2 0.04000 1.08420280842809 3 0.06000 1.12970427426891 4 0.08000 1.17771907215813 5 0.10000 1.22846046963387 6 0.12000 1.28216861455636 7 0.14000 1.33911486875439 8 0.16000 1.39960701022635 9 0.18000 1.46399551405448 10 0.20000 1.53268118238470 Estremo b Valore finale Valore esatto Errore relativo 0.2000 1.532681182 1.532681180 0.15559E-08