AZZERAMENTO DI UNA FUNZIONE
NON LINEARE CON METODO DICOTOMICO
Contenuto
Script per l'azzeramento di una funzione non lineare
Testo del problema
La banca dati DIPPR riporta i seguenti dati per la pressione di vapore
dell'acqua:
For WATER
The Vapor Pressure can be calculated as follows:
VP in Pa = exp(A + B/T + C*ln(T) + D*T^E)
Where: T = Temperature in Kelvin
A = 7.3649E+01 B = -7.2582E+03 C = -7.3037E+00
D = 4.1653E-06 E = 2.0000E+00
In the range: 273.16 K to 647.13 K
Range is experimental Quality code: 1
Determinare, tramite metodo dicotomico, a quale temperatura si raggiunge
una tensione di vapore pari a mezza atmosfera.
Confrontare il risultato ottenuto con quello prodotto dalla funzione
"fsolve" di Matlab dedicata alla risoluzione di sistemi di equazioni non
lineari. Esempio di sintassi:
[x,fval,exitFlag] = fsolve(@f,x0);
Dove:
f è la funzione da azzerare, x0 il valore di primo tentativo,
x la soluzione, fval il valore della funzione in corrispondenza di x ed
exitFlag l'indice di stato all'uscita di fsolve: > 0 se ha raggiunto
la convergenza, = 0 se è stato raggiunto il numero massimo di
iterazioni, < 0 se non è stata individuata la soluzione.
N.B.: il simbolo @ prima del nome della funzione da azzerare rappresenta
il puntatore a tale funzione.
Inizializzazione dell'ambiente
clear all
clc
Definizione delle grandezze globali
Le grandezze AA BB ... vengono contemporaneamente utilizzate nel
presente script e nella funzione "funz" da azzerare. Il costrutto GLOBAL
permette di NON passare via lista le costanti/variabili ad altre
porzioni di codice.
Si noti che le variabili/costanti sono separate da uno spazio.
Le costanti/variabili presenti nel costrutto GLOBAL hanno lo stesso
ordine nelle varie porzioni di codice in cui tale costrutto compare.
global AA BB CC DD EE PRESSIONE
Inizializzazione coefficienti pressione di vapore acqua secondo banca dati DIPPR
Le costanti AA BB ... verranno utilizzate dalla funzione "funz"
AA = 7.3649E+01;
BB = -7.2582E+03;
CC = -7.3037E+00;
DD = 4.1653E-06;
EE = 2.0000E+00;
PRESSIONE = 0.5 * 101325.;
Valori di primo tentativo e altre grandezze di lavoro
a = 283.15;
b = 373.15;
epsi = 1.e-5;
MAX_ITER = 30;
iConto = 0;
fa = funz(a);
fb = funz(b);
Controllo condizione necessaria applicabilità metodo dicotomico...
if sign(fa) * sign(fb) > 0
disp('Estremi del metodo dicotomico NON consistenti. Impossibile continuare...');
return
end
Inizio ciclo iterativo
while abs(b-a) > epsi & iConto < MAX_ITER
iConto = iConto + 1;
c = a + 0.5 * (b-a);
fc = funz(c);
if sign(fa) * sign(fc) < 0.
b = c;
fb = fc;
else
a = c;
fa = fc;
end
disp([num2str(iConto),') c = ',num2str(c),' f(c) = ',num2str(fc),' delta = ',num2str(b-a)]);
end
1) c = 328.15 f(c) = -34902.7672 delta = 45
2) c = 350.65 f(c) = -7886.3652 delta = 22.5
3) c = 361.9 f(c) = 16169.3852 delta = 11.25
4) c = 356.275 f(c) = 3021.643 delta = 5.625
5) c = 353.4625 f(c) = -2692.0573 delta = 2.8125
6) c = 354.8687 f(c) = 97.4146 delta = 1.4063
7) c = 354.1656 f(c) = -1313.8542 delta = 0.70313
8) c = 354.5172 f(c) = -612.3917 delta = 0.35156
9) c = 354.693 f(c) = -258.5364 delta = 0.17578
10) c = 354.7809 f(c) = -80.8235 delta = 0.087891
11) c = 354.8248 f(c) = 8.2298 delta = 0.043945
12) c = 354.8028 f(c) = -36.3133 delta = 0.021973
13) c = 354.8138 f(c) = -14.0458 delta = 0.010986
14) c = 354.8193 f(c) = -2.909 delta = 0.0054932
15) c = 354.8221 f(c) = 2.6601 delta = 0.0027466
16) c = 354.8207 f(c) = -0.12451 delta = 0.0013733
17) c = 354.8214 f(c) = 1.2678 delta = 0.00068665
18) c = 354.821 f(c) = 0.57164 delta = 0.00034332
19) c = 354.8209 f(c) = 0.22356 delta = 0.00017166
20) c = 354.8208 f(c) = 0.049523 delta = 8.5831e-005
21) c = 354.8207 f(c) = -0.037495 delta = 4.2915e-005
22) c = 354.8207 f(c) = 0.0060141 delta = 2.1458e-005
23) c = 354.8207 f(c) = -0.015741 delta = 1.0729e-005
24) c = 354.8207 f(c) = -0.0048632 delta = 5.3644e-006
Soluzione finale
disp(' ');
disp(['Soluzione finale raggiunta dopo ',num2str(iConto),' iterazioni']);
disp([' xSol = ',num2str(c),' = ',num2str(c-273.15),' ^C f(xSol) = ',num2str(fc)]);
disp(' ');
Soluzione finale raggiunta dopo 24 iterazioni
xSol = 354.8207 = 81.6707 ^C f(xSol) = -0.0048632
Verifica con il risolutore non lineare di Matlab "fsolve"
x0 = 50. + 273.15;
[x,fval,exitFlag] = fsolve(@funz,x0);
disp(['Soluzione con "fsolve" x = ',num2str(x),' f(x) = ',num2str(fval)]);
disp(['Indice di uscita da fsolve: ',num2str(exitFlag),' (se maggiore di 0 vuol dire che ha raggiunto la convergenza !!!)']);
Optimization terminated: first-order optimality is less than options.TolFun.
Soluzione con "fsolve" x = 354.8207 f(x) = 4.3656e-010
Indice di uscita da fsolve: 1
(se maggiore di 0 vuol dire che ha raggiunto la convergenza !!!)
Grafico tensione di vapore dell'acqua
Si noti la monotonia crescente di tale relazione
for i = 1: 100
T(i) = i + 273.15;
P(i) = funz(T(i)) + PRESSIONE;
end
plot(T,P,'r-o'),grid,xlabel('temperatura [K]'),ylabel('Pv(T) [Pa]'),title('Tensione di vapore dell''acqua');
Funzione da azzerare
Si noti l'utilizzo dei coefficienti necessari alla valutazione della
tensione di vapore passati tramite istruzione GLOBAL
function f = funz(T)
global AA BB CC DD EE PRESSIONE
f = exp(AA + BB / T + CC * log(T) + DD * T ^ EE) - PRESSIONE;