I want to numerically integrate an integral with infinite limit. Does anyone have any idea how should I do that?
int(x* exp (v*x + (1-exp(v*x))/v),x, o , inf) does not work.
Note that I will have values for v.
%n=10;
kappa=.5;
delta0=.5;
Vmax=500;
Vdep=2.2;
l=2.2;
kbT=4.1;
%xb=.4;
fb=10;
k=1;
V0=5;
e1=(fb*l/kbT)*(kappa/delta0);
e2=Vmax/V0;
e3=Vdep/V0;
w=zeros(1,25);
for v=1:25
w(:,v)=integral(@(x) x.*exp(v*x+((1-exp(v*x))/v)),0,inf);
end
e12=e2*exp(-e1*(1:25).*w.^2)-e3;
plot(e12);
ylim([0 25]);
hold on;
plot(0:25,0:25);
xlim([0 25]);
%hold off;
The plot is not matching the real data in the article!(for the e12 curve specially) I need to calculate the intersection of the 2 curves (which is ~13.8 based on the paper) and then in the second part I have to add a term in e12 which contains an independent variable:
v=13.8;
w= integral(@(x) x.*exp(v*x+((1-exp(v*x))/v)),0,inf)
e4 = zeros (1,180);
fl = 1:180;
e4(:,fl)= (fl*l/kbT)*(kappa/n);
e12=e2*exp(-e1*v*w^2-e4)-e3
But again the problem is that running this code I will end with a negative value for e12 which should be just approaching zero in large values of fl (fl>160)
to show how this code is different from the expected curve you can plot these data on the same figure:
fl = [0, 1, 4, 9, 15, 20, 25, 40, 60, 80, 100, 120, 140, 160, 180];
e12 = [66, 60, 50, 40, 30, 25.5, 20, 15.5, 10.5, 8.3, 6.6, 5, 2.25, 1.1, 0.5];
which obviously does not match the curve generated by the code.
Assuming the question is about this full code:
syms x;
v = 1; % For example
int(x*exp(v*x + (1-exp(v*x))/v),x, 0, Inf)
and the issue is that it returns itself (i.e., int doesn't find an analytic solution), one can set the 'IgnoreAnalyticConstraints' option to true (more details) to get a solution:
syms x;
v = 1; % For example
int(x*exp(v*x + (1-exp(v*x))/v),x, 0, Inf, 'IgnoreAnalyticConstraints', true)
returns -ei(-1)*exp(1), where ei is the exponential integral function (see also expint for numerical calculations). For negative values of v the solution will also be in terms of eulergamma, the Euler-Mascheroni constant. And of course the integral is undefined if v is 0.
Using Mathematica 10.0.2's Integrate yields a full solution for symbolic v.
Integrate[x Exp[v x - (Exp[v x] - 1)/v], {x, 0, Infinity}]
returns
ConditionalExpression[(E^(1/v) (EulerGamma + Gamma[0, 1/v] + Log[1/v]))/v, Re[v] < 0]
Applying Assumptions:
Integrate[x Exp[v x - (Exp[v x] - 1)/v], {x, 0, Infinity}, Assumptions -> v > 0]
Integrate[x Exp[v x - (Exp[v x] - 1)/v], {x, 0, Infinity}, Assumptions -> v < 0]
returns
(E^(1/v) Gamma[0, 1/v])/v
and
(E^(1/v) (2 EulerGamma - 2 ExpIntegralEi[-(1/v)] + Log[1/v^2]))/(2 v)
where Gamma is the upper incomplete gamma function. These appear match up with the results from Matlab.
To evaluate these numerically in Matlab:
% For v > 0
v_inv = 1./v;
exp(v_inv).*expint(v_inv).*v_inv
or
% For v < 0
v_inv = 1./v;
exp(v_inv).*(2*double(eulergamma)+2*(expint(v_inv)+pi*1i)+log(v_inv.^2)).*v_inv/2