How to fix infinity limit in the integral equation
15 views (last 30 days)
Show older comments
Kumaresh Kumaresh
on 16 Jul 2022
Commented: Kumaresh Kumaresh
on 11 Aug 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 Comment
Bruno Luong
on 16 Jul 2022
Edited: Bruno Luong
on 16 Jul 2022
x = linspace(E0,Inf)
Nice try. You should think more about that.
Accepted Answer
Bruno Luong
on 16 Jul 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
More Answers (3)
Walter Roberson
on 16 Jul 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
on 9 Aug 2022
2 Comments
Walter Roberson
on 9 Aug 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
See Also
Categories
Find more on Calculus in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!