How to fix infinity limit in the integral equation
15 ビュー (過去 30 日間)
古いコメントを表示
Kumaresh Kumaresh
2022 年 7 月 16 日
コメント済み: Kumaresh Kumaresh
2022 年 8 月 11 日
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!