Generalized exponential integral with negative argument

Is there a way to compute the generalized exponential integral E[a,x]=int_1^infty e^(-x t) t^a dt with negative x? In particular, I need to compute the integral int_1^z Exp[-y^b/b] dy where b<0 and z>1. (I need to do this for 1000’s of z-values so numerical integration is not an option) Mathematica tells me this equals an expression involving ExpIntegralE[(-1 + b)/b, 1/b] (to be precise ExpIntegralE[(-1 + b)/b, 1/b] - z ExpIntegralE[(-1 + b)/b, z^b/b])/b. So since b<0, x=1/b is negative, and this causes a problem.

There seems to be a related thread but with no answer https://www.mathworks.com/matlabcentral/newsreader/view_thread/289192

Finally, there is an implementation of E[a,x] on file exchange https://www.mathworks.com/matlabcentral/fileexchange/52694-generalised-exponential-integral but it cannot handle negative x-values.

Answer 1

[Star Strider]

See if expint will do what you want. There’s also a Symbolic Math Toolbox function by the same name, so search for it if you want it instead.

The gammainc function computes the incomplete gamma function.

[Benjamin Moll]

Thanks for looking into this! The problem is that the matlab expint function only does the plain-vanilla exponential integral function expint(x) and not the generalized expint(a,x) one, see http://www.mathworks.com/help/matlab/ref/expint.html

Instead, the symbolic version of expint http://www.mathworks.com/help/symbolic/expint.html does indeed have an option to complete the generalized exponential integral expint(a,x). But I need to embed the computation in a larger code with many loops etc (solving some PDEs via finite difference methods) and so I was hoping for a standard matlab implementation.

[Star Strider]

My pleasure.

The only other option I can see is to code it as a series expansion. This may require some research, since my ancient copy of Abramowitz and Stegun doesn’t even mention the generalized exponential integral. There’s an exhaustive discussion of them in Introduction to the exponential integrals. That’s the best I can do.

[Benjamin Moll]

Thanks for the reference, that looks useful. Quick thought: is there a way to call the symbolic expint(a,x) from within a standard matlab script? Sorry if this is a stupid question -- I have zero experience with the symbolic toolbox.

[Star Strider]

You would have to do at least some of the calculations in the Symbolic Math Toolbox. Probably the easiest way to do this is to create a separate function file to invoke the Symbolic Math Toolbox function, do the calculation, and return the numerical result to your calling script.

Example:

function y = SymbolicEi(n,x)
ns = sym(n);
xs = sym(x);
y  = vpa(expint(ns,xs));
end

I tested this and it returns a complex numeric result for negative ‘x’. Save it as SymbolicEi.m (change the function and file names to the same name if you want to save it as a different function name), then call it from your main script as you would any other function, for example:

y = SymbolicEi(5,-1)

[Walter Roberson]

I suggest using double() instead of vpa() if you want to return a numeric result.

[Benjamin Moll]

Great thank you Star Strider and Walter. Both formulas work and give the same answer. One question left: I need to do this for a vector of 1000 x's, e.g. x=linspace(1,2,1000). The way you guys do it, it cannot handle vectors. Is there a way to get this to work without writing a loop?

[Star Strider]

My (our) pleasure!

The inability to handle vectors is probably inherent in the symbolic expint code. The Symbolic Math Toolbox functions code are probably not vectorised because symbolic computations are usually one-offs that yield one specific result. Sometimes using the matlabFunction function works, but this time it isn’t terribly helpful:

syms ns xs
y  = matlabFunction(expint(ns,xs))

y = @(ns,xs)expint(sym(ns),sym(xs))

(I put the function on the output on the same line for convenience.)

So the loop is probably the only option.

[Benjamin Moll]

Sorry for the late reply, I was traveling. It all works now, thanks so much. I ended up doing a sort of convex combination of Walter’s and Star Striders suggestions. Here's what I did:

1) I wrote a function

function y = SymbolicEi(n,x)
syms ns nx
G = expint(ns,nx);
GG = subs(G, {ns, nx}, {n, x});
y = double(GG);
end

2) I called it as

b=-1.01; z = linspace(1,2,100);
(SymbolicEi((-1 + b)/b, 1/b)*ones(1,length(z)) - z.*SymbolicEi((-1 + b)/b, z.^b/b))/b

The advantage over Walter’s is that it’s less messy and works also for integer 1/b. I'm giving randomly giving the answer credit to StarStrider. But really it should go to both of you. Thanks so much!

[Star Strider]

My (our) pleasure!

Add your vote to mine for Walter’s Answer, and Walter gets the same 4 RP’s I got.

[Walter Roberson]

The defining formula for Ei with two arguments is in terms of the Gamma function, which has singularities at every negative integer, so if your 1/b is a negative integer you should be expecting a singularity.

Answer 2

[Walter Roberson]

This formula appears to work. It requires the symbolic toolbox

(-z^(b+1)*hypergeom([(1/2)*(b+1)/b], [3/2, (1/2)*(3*b+1)/b], (1/4)*z^(2*b)/b^2)+b*z*(b+1)*hypergeom([(1/2)/b], [1/2, (1/2)*(2*b+1)/b], (1/4)*z^(2*b)/b^2)+hypergeom([(1/2)*(b+1)/b], [3/2, (1/2)*(3*b+1)/b], (1/4)/b^2)+(-b^2-b)*hypergeom([(1/2)/b], [1/2, (1/2)*(2*b+1)/b], (1/4)/b^2))/(b*(b+1))

Note: this might fail if 1/b is an integer

Reference: Maple

[Benjamin Moll]

Great thank you Star Strider and Walter. Both formulas work and give the same answer. One question left: I need to do this for a vector of 1000 z's, e.g. z=linspace(1,2,1000). The way you guys do it, it cannot handle vectors. Is there a way to get this to work without writing a loop?

[Walter Roberson]

syms b z
G = (-z^(b+1)*hypergeom([(1/2)*(b+1)/b], [3/2, (1/2)*(3*b+1)/b], (1/4)*z^(2*b)/b^2)+b*z*(b+1)*hypergeom([(1/2)/b], [1/2, (1/2)*(2*b+1)/b], (1/4)*z^(2*b)/b^2)+hypergeom([(1/2)*(b+1)/b], [3/2, (1/2)*(3*b+1)/b], (1/4)/b^2)+(-b^2-b)*hypergeom([(1/2)/b], [1/2, (1/2)*(2*b+1)/b], (1/4)/b^2))/(b*(b+1));
GG = subs(G, {z, b}, {linspace(1,2,1000), -sqrt(pi)});
GGn = double(GG);