How to fix infinity limit in the integral equation
28 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
Kumaresh Kumaresh
il 16 Lug 2022
Commentato: Kumaresh Kumaresh
il 11 Ago 2022
Hello everyone, Hope everyone is doing well.
Here is my code attached below.
The code is solved using trapezoidal integration. The variable x in the code need to be fixed with the limits ranging from E0 to Infinity. For the sake of testing the case, the variable x is fixed with the limits ranging from E0 to EA+210000, and the case works good but as expected the outcomes are not satisfied.
My query:
How to fix the upper limit of the integral with infinity ? By fixing infinity, the code returns Nan :-(.
Thank you
%diary verse_15thJuly
clear; clc; close all;
% INITIAL CONDITIONS
dtemp=1; finaltemp=1050; dt=(dtemp*60/3);
H2O=0.023; CO2=0.00555;
fprintf('T\t Totvol\n');
for T = 350:dtemp:finaltemp
TK=T+273; E0 = 183000;
EA = ((189.95*log(T)) - 1013.5)*1000;
RT=8.3144*TK;
%x = linspace(E0,EA+210000);
x = linspace(E0,Inf);
FE_H2O = exp(-((((x)-E0)/48000).^8));
FE_CO2 = exp(-((((x)-E0)/78000).^4));
A_PVM = 1-exp(-1.3E13.*dt.*exp(-(x)./RT));
H2Ox = -H2O*(-8/(48000^8)).*(((x)-E0).^7).*FE_H2O.*A_PVM;
CO2x = -CO2*(-4/(78000^4)).*(((x)-E0).^3).*FE_CO2.*A_PVM;
H2Oy = H2O - trapz(x,H2Ox);
%H2Oy = H2O - expint(H2Ox);
CO2y = CO2 - trapz(x,CO2x);
Y=H2Oy+CO2y+0.97145;
VM = 1- Y;
fprintf('%0.0f\t %0.8f\n',[T,VM]);
end
%diary off
1 Commento
Bruno Luong
il 16 Lug 2022
Modificato: Bruno Luong
il 16 Lug 2022
x = linspace(E0,Inf)
Nice try. You should think more about that.
Risposta accettata
Bruno Luong
il 16 Lug 2022
Use integral function rather trapz
% INITIAL CONDITIONS
dtemp=1; finaltemp=1050; dt=(dtemp*60/3);
H2O=0.023; CO2=0.00555;
fprintf('T\t Totvol\n');
for T = 350:dtemp:finaltemp
TK=T+273; E0 = 183000;
EA = ((189.95*log(T)) - 1013.5)*1000;
RT=8.3144*TK;
FE_H2O_fun = @(x) exp(-((((x)-E0)/48000).^8));
FE_CO2_fun = @(x) exp(-((((x)-E0)/78000).^4));
A_PVM_fun = @(x) 1-exp(-1.3E13.*dt.*exp(-(x)./RT));
H2Ox_fun = @(x) -H2O*(-8/(48000^8)).*(((x)-E0).^7).*FE_H2O_fun(x).*A_PVM_fun(x);
CO2x_fun = @(x) -CO2*(-4/(78000^4)).*(((x)-E0).^3).*FE_CO2_fun(x).*A_PVM_fun(x);
H2Oy = H2O - integral(H2Ox_fun, E0, Inf);
CO2y = CO2 - integral(CO2x_fun, E0, Inf);
Y=H2Oy+CO2y+0.97145;
VM = 1- Y;
fprintf('%0.0f\t %0.8f\n',[T,VM]);
end
Più risposte (3)
Walter Roberson
il 16 Lug 2022
With the various exp(-x) terms as x approaches infinity the exp() terms go to 0. But the exp() terms are being multiplied by a polynomial in x, and as x goes to infinity the polynomial goes to either +inf or -inf (you would need further analysis to figure out which.) But when you are working in double precision, 0 times infinity gives nan.
If you switch over to the symbolic toolbox and use symbolic x, and use symbolic upper bound, then as the upper bound goes to infinity, VM goes to 0
Kumaresh Kumaresh
il 9 Ago 2022
2 Commenti
Walter Roberson
il 9 Ago 2022
With numeric integration you would encounter the 0 times infinity problem.
Consider for example,
f = @(x) exp(-x.^2) .* exp(x)
f(750)
f(sym(750))
exp(-x^2) obviously goes to 0 faster than exp(x) goes to infinity as x increases, so we can see that mathematically the limit has to be 0 -- but in floating point you are going to get 0 * infinity
Vedere anche
Categorie
Scopri di più su Calculus in Help Center e File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!