How to numerically integrate with infinite limit in MATLAB?

Question

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.

"Does not work" is not helpful. If there is an error, provide it in full. If the output does not match what you expect, show what it is, why you think it's wrong. Edit your question. Also provide fully runnable code, including all variable declarations. — horchler, Jun 05 '15 at 03:53
Also, is your lower integration bound really supposed to be a variable called `o` (letter) and not the number `0` (number)? — horchler, Jun 05 '15 at 04:01
Hi,Thanks for your help! 1. by does not work I meant it goes into error! I can't remember what was the error, but now I have changed the code and it is now successfully running but the results are somehow different from what I expect to be. (I will provide the code in the next comment) 2.. the lower integration ix 0 and it was a typo when I was typing it here! 3. the values for v are real integers ranging from 0 to 25. — Roza, Jun 07 '15 at 16:32

score -1 · Answer 1 · answered Jun 05 '15 at 04:34

Numerical integration is performed by summing the function at discrete points with distance dx. The smaller dx you choose, the better approximation you get. For example integrating from x=0 to x=10 is done by:

x = 0:dx:10;
I = sum(x.* exp (v*x + (1-exp(v*x))/v))*dx;

obviously, you can't do that for x=inf. But I believe you function decays rapidly. Therefore, you can assume that x* exp (v*x + (1-exp(v*x))/v) = 0 for large enough x. Otherwise the integral is divergent. So all you have to do is set the limit on x. If you are not sure what the limit should be, you can perform a loop with a stopping condition:

I = 0;
prevI = -1;
x = 0;
while abs(I-prevI)>err
  prevI = I;
  I = I + x.* exp (v*x + (1-exp(v*x))/v)*dx;
  x = x + dx;
end

Now, all you have to do is set the desired dx and err

While this does work, not using Matlab builtin function tend to be painfully slow. Not to mention using this kind of integration is very unprecise. Some methods have a much better precision without a lot more of computing cost. — meneldal, Jun 05 '15 at 05:11
Previously, I did the same method for approximation but the problem is that the results are different form what I expect it to be!I don't know what is the problem. How can I sent a photo here, so that I can show you what I mean by different results. — Roza, Jun 07 '15 at 16:48

score -1 · Answer 2 · edited Jun 05 '15 at 05:25

-1

You must read this:Mathwork Link

perhaps you are making a mistake in the function that you use. Also note that MATLAB syntax is case sensitive..

edited Jun 05 '15 at 05:25

Abdulla Nilam

36,589
17
64
85

answered Jun 05 '15 at 04:35

Devraj Dutt

1

horchler · Accepted Answer · 2015-06-07T20:06:55.270

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

How to numerically integrate with infinite limit in MATLAB?

3 Answers3