How to fix infinity limit in the integral equation
28 次查看(过去 30 天)
显示 更早的评论
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 个评论
Bruno Luong
2022-7-16
编辑:Bruno Luong
2022-7-16
x = linspace(E0,Inf)
Nice try. You should think more about that.
采纳的回答
Bruno Luong
2022-7-16
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
更多回答(3 个)
Walter Roberson
2022-7-16
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
2022-8-9
2 个评论
Walter Roberson
2022-8-9
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
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Calculus 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!