How to fix infinity limit in the integral equation
28 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Kumaresh Kumaresh
el 16 de Jul. de 2022
Comentada: Kumaresh Kumaresh
el 11 de Ag. de 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 comentario
Bruno Luong
el 16 de Jul. de 2022
Editada: Bruno Luong
el 16 de Jul. de 2022
x = linspace(E0,Inf)
Nice try. You should think more about that.
Respuesta aceptada
Bruno Luong
el 16 de Jul. de 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
Más respuestas (3)
Walter Roberson
el 16 de Jul. de 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
el 9 de Ag. de 2022
2 comentarios
Walter Roberson
el 9 de Ag. de 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
Ver también
Categorías
Más información sobre Calculus en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!