How to Integrate this function?

I used the following matlab code to integrate f equation as follows:

syms v C sigma v0

f=((exp(-C*v)))*(1/(sigma*(pi/2)^0.5*erfc(-v0/(2^0.5*sigma))))*(exp(-(-((v-v0)^2)/(2*sigma^2))));

int(f,v,[0 inf])

to get this previously reported answer:

f=(erfc((-1/2^0.5)*((v0/sigma)-(C*sigma))))/(erfc((-1/2^0.5)*(v0/sigma)))*exp((-C*v0)+(0.5*sigma^2*C^2));

However, I get this, which is not correct:

ans =

-(4503599627370496*((2^(1/2)*pi^(1/2)*exp(- (C^2*sigma^2)/2 - v0*C)*limit(erf((2^(1/2)*(((C*sigma^2 + v0)*1i)/sigma^2 - (v*1i)/sigma^2))/(2*(1/sigma^2)^(1/2))), v, Inf)*1i)/(2*(1/sigma^2)^(1/2)) - (2^(1/2)*pi^(1/2)*exp(- (C^2*sigma^2)/2 - v0*C)*erf((2^(1/2)*(C*sigma^2 + v0)*1i)/(2*sigma^2*(1/sigma^2)^(1/2)))*1i)/(2*(1/sigma^2)^(1/2))))/(5644425081792261*sigma*(erfc((2^(1/2)*v0)/(2*sigma)) - 2))

Would you please help me on this?

Answer 1

[Walter Roberson]

The "previously reported answer" is incorrect, and MATLAB's answer is correct. The integral is infinite if C and sigma are nonnegative.

It appears to me that the "previously reported answer" might have been calculated by Maple. It appears to be to be a bug in Maple, which I will report to Maplesoft.

If you are using Maple, then if you substitute numeric values for v0, sigma, and C and then do the integration over 0 to infinity, then you will get infinity as the output; whereas in Maple when the bug is in effect if you integrate over v and then substitute in numeric values you will get a finite symbolic solution. It looks to me as if that finite solution that you obtain might be the integral from 0 to 1. If you do the same thing bug integrate from 1 to infinity in Maple you will get infinity -- and since the integral from 0 to 1 is finite and non-negative, the whole integral must be infinite.

[Walter Roberson]

Look very carefully at the Pneg they show first, Appendix A, belonging to "(6)". It is an integral from negative infinity to positive infinity, and the part with (v-v0)^2 has two negative signs that cancel out,

exp(- (-(v-v0)^2/(2*sigma^2)) )

when the quantities involved are real-valued (and possibly even when they are complex), this is

exp( ((v-v0)^2/(2*sigma^2)) )

The part inside the exp() gets indefinitely large as v approaches infinity, and exp() of that goes to infinity. The counter-balancing part has a -C*v but as v increases, v^2 is always going to overwhelm -C*v. So the integral from 0 to infinity would be infinite.

The key here is that they integrate from negative infinity to positive infinity, in which case there is the hypothetical potential for infinities "balancing out" to produce a finite result over the entire -inf to +inf range -- but when you integrate over 0 to inf there is no potential for canceling out.

Now look closely at the Pneg they give in Appendix B. It is an integral from 0 to positive infinity, and it has a single negative sign,

exp(- ((v-v0)^2/(2*sigma^2)) )

The part inside the exp() gets indefinitely negative as v approaches infinity, and exp() of that goes to 0, leading to a finite integral.

Now if you examine your code, you will see that you used the integral form with the two negatives (that cancel) but you used it over 0 to infinity.

[Fereshteh Emami Piri]

I have corrected it accordingly, with one negative over 0 and infinity. I still don't get the same answer with their result!

[Fereshteh Emami Piri]

Actually, my goal is to use "int" command to solve the following equation that is the previous equation they solved along with an extra ((C*v)^n)/factorial(n) part:

syms v C sigma v0 n 
f=((exp(-C*v))*(((C*v)^n)/factorial(n)))*(1/(sigma*(pi/2)^0.5*erfc(-v0/(2^0.5*sigma))))*(exp(-(((v-v0)^2)/(2*sigma^2)))); 
int(f,v,[0 inf])

And this is the answer I got:

int(-(4503599627370496*exp(-C*v)*exp(-(v - v0)^2/(2*sigma^2))*(C*v)^n)/(5644425081792261*sigma*factorial(n)*(erfc((2^(1/2)*v0)/(2*sigma)) - 2)), v, 0, Inf)

Any guide would be helpful.

[Walter Roberson]

After the change, the symbolic result I get is

(exp(-C*v0)*(exp((C^2*sigma^2)/2)*limit(erf((2^(1/2)*(C*sigma^2*1i + v*1i - v0*1i))/(2*sigma^2*(-1/sigma^2)^(1/2))), v, Inf)*1i + erf((2^(1/2)*(- C*sigma^2*1i + v0*1i))/(2*sigma^2*(-1/sigma^2)^(1/2)))*exp((C^2*sigma^2)/2)*1i))/(sigma*(erf((2^(1/2)*v0)/(2*sigma)) + 1)*(-1/sigma^2)^(1/2))

Although this is not in a nice form, I find it does produce the right answer.

For example,

syms v C sigma v0 nonnegative
Pi = sym('pi');
f=((exp(-C*v)))*(1/(sigma*(Pi/2)^sym(1/2)*erfc(-v0/(sym(2)^sym(1/2)*sigma))))*(exp((-((v-v0)^sym(2))/(sym(2)*sigma^sym(2)))));

simplify(int(f,v,[0 inf]))
simplify(subs(ans,[v0,sigma,C],[19,2,11]))
double(ans)

The formula reached after the subs() of those particular values is exactly the same as I get in Maple.

The symbolic form I get from Maple is

-exp((1/2)*C*(C*sigma^2-2*v0))*(erf(sqrt(2)*(C*sigma^2-v0)/(2*sigma))-1)/(1+erf(v0*sqrt(2)/(2*sigma)))

MATLAB is sometimes a bit weak in evaluating limits -- and this limit is a bit tricky since it involves erf and complex numbers.

[Walter Roberson]

For the version involving factorial, the following helps slightly:

syms v C sigma v0 n nonnegative
Pi = sym('pi');
f=((exp(-C*v))*(((C*v)^n)/factorial(n)))*(1/(sigma*(Pi/2)^sym(1/2)*erfc(-v0/(2^sym(1/2)*sigma))))*(exp(-(((v-v0)^2)/(2*sigma^2)))); 
simplify(int(f,v,[0 inf]))

but this still has an int() in it.

Maple is able to resolve it to

-(1/4)*Pi^(3/2)*exp(-(1/2)*v0^2/sigma^2)*2^((1/2)*n)*(sin((1/2)*Pi*n)*GAMMA(-(1/2)*n+3/2)*(C^4*sigma^8+(-4*C^3*v0+2*C^2*(n-1))*sigma^6+(6*C^2*v0^2-4*C*(n-1)*v0+n-1)*sigma^4+(-4*C*v0^3+(-2+2*n)*v0^2)*sigma^2+v0^4)*LaguerreL(-(1/2)*n+1/2, 1/2, (1/2)*(C*sigma^2-v0)^2/sigma^2)-(sin((1/2)*Pi*n)*GAMMA(-(1/2)*n+3/2)*(C^2*sigma^4+(-2*C*v0+n-1)*sigma^2+v0^2)*(C*sigma^2-v0)*LaguerreL(-(1/2)*n+1/2, 3/2, (1/2)*(C*sigma^2-v0)^2/sigma^2)-2^(1/2)*sigma*cos((1/2)*Pi*n)*GAMMA(-(1/2)*n+2)*((C^2*sigma^4+(-2*C*v0+n)*sigma^2+v0^2)*LaguerreL(-(1/2)*n, 1/2, (1/2)*(C*sigma^2-v0)^2/sigma^2)-LaguerreL(-(1/2)*n, 3/2, (1/2)*(C*sigma^2-v0)^2/sigma^2)*(C*sigma^2-v0)^2))*(C*sigma^2-v0))*sigma^(n-4)*C^n/(pi^(1/2)*GAMMA(1+n)*sin((1/2)*Pi*n)*(1+erf((1/2)*v0*2^(1/2)/sigma))*GAMMA(-(1/2)*n+3/2)*cos((1/2)*Pi*n)*GAMMA(-(1/2)*n+2))

MATLAB does not have LaguerreL .

The equivalent to the above in terms of hypergeom is

-(1/30)*((n-1)*GAMMA(-(1/2)*n+2)*(C^4*sigma^8+(-4*C^3*v0+2*C^2*(n-1))*sigma^6+(6*C^2*v0^2-4*C*(n-1)*v0+n-1)*sigma^4+(-4*C*v0^3+(-2+2*n)*v0^2)*sigma^2+v0^4)*sin((1/2)*Pi*n)*(C*sigma^2-v0)^2*hypergeom([1/2+(1/2)*n], [7/2], (1/2)*(C*sigma^2-v0)^2/sigma^2)-(1/2)*sigma*(-10*(n-1)*sigma*GAMMA(-(1/2)*n+2)*(C^4*sigma^8+(-4*C^3*v0+C^2*(n+2))*sigma^6+(3+6*C^2*v0^2-2*C*(n+2)*v0)*sigma^4+(-4*C*v0^3+(n+2)*v0^2)*sigma^2+v0^4)*sin((1/2)*Pi*n)*hypergeom([(1/2)*n-1/2], [5/2], (1/2)*(C*sigma^2-v0)^2/sigma^2)+GAMMA(-(1/2)*n+3/2)*2^(1/2)*(n-2)*n*((C^2*sigma^4+(-2*C*v0+n)*sigma^2+v0^2)*(C*sigma^2-v0)^2*hypergeom([1+(1/2)*n], [7/2], (1/2)*(C*sigma^2-v0)^2/sigma^2)+5*sigma^2*hypergeom([(1/2)*n], [5/2], (1/2)*(C*sigma^2-v0)^2/sigma^2)*(C^2*sigma^4+(-2*C*v0+3)*sigma^2+v0^2))*cos((1/2)*Pi*n)*(C*sigma^2-v0)))*Pi*sigma^(n-6)*exp(-(1/2)*v0^2/sigma^2)*2^((1/2)*n)*C^n/(pi^(1/2)*GAMMA(1+n)*(1+erf((1/2)*v0*2^(1/2)/sigma))*cos((1/2)*Pi*n)*GAMMA(-(1/2)*n+3/2)*sin((1/2)*Pi*n)*GAMMA(-(1/2)*n+2))

[Fereshteh Emami Piri]

Dear Walter,

The last solution for the version without factorial part works. However the solution for the version with factorial part you obtained, for any number of n (n=0, 1, 2,...,100), I get NaN!:

-(1/30)*((n-1)*GAMMA(-(1/2)*n+2)*(C^4*sigma^8+(-4*C^3*v0+2*C^2*(n-1))*sigma^6+(6*C^2*v0^2-4*C*(n-1)*v0+n-1)*sigma^4+(-4*C*v0^3+(-2+2*n)*v0^2)*sigma^2+v0^4)*sin((1/2)*Pi*n)*(C*sigma^2-v0)^2*hypergeom([1/2+(1/2)*n], [7/2], (1/2)*(C*sigma^2-v0)^2/sigma^2)-(1/2)*sigma*(-10*(n-1)*sigma*GAMMA(-(1/2)*n+2)*(C^4*sigma^8+(-4*C^3*v0+C^2*(n+2))*sigma^6+(3+6*C^2*v0^2-2*C*(n+2)*v0)*sigma^4+(-4*C*v0^3+(n+2)*v0^2)*sigma^2+v0^4)*sin((1/2)*Pi*n)*hypergeom([(1/2)*n-1/2], [5/2], (1/2)*(C*sigma^2-v0)^2/sigma^2)+GAMMA(-(1/2)*n+3/2)*2^(1/2)*(n-2)*n*((C^2*sigma^4+(-2*C*v0+n)*sigma^2+v0^2)*(C*sigma^2-v0)^2*hypergeom([1+(1/2)*n], [7/2], (1/2)*(C*sigma^2-v0)^2/sigma^2)+5*sigma^2*hypergeom([(1/2)*n], [5/2], (1/2)*(C*sigma^2-v0)^2/sigma^2)*(C^2*sigma^4+(-2*C*v0+3)*sigma^2+v0^2))*cos((1/2)*Pi*n)*(C*sigma^2-v0)))*Pi*sigma^(n-6)*exp(-(1/2)*v0^2/sigma^2)*2^((1/2)*n)*C^n/(pi^(1/2)*GAMMA(1+n)*(1+erf((1/2)*v0*2^(1/2)/sigma))*cos((1/2)*Pi*n)*GAMMA(-(1/2)*n+3/2)*sin((1/2)*Pi*n)*GAMMA(-(1/2)*n+2))

I expect the answer 0.0822 for n=0, v0=1; sigma=0.005, as it works for the one without factorial part. With solving the integrated solution, the expectation for any value of "n" is a value below 1, but NaN! Would you let me know your opinion?

[Walter Roberson]

What value of C are you using?

You are correct that the GAMMA / hypergeom form gives errors for all nonnegative integer n. Maple failed to take the possibility of integer n into account in building the formula. Unfortunately I am having difficulty finding a closed form version that takes integers into account.

[Fereshteh Emami Piri]

I use C=2.5. So, this means I cannot use MATLAB and even Maple to integrate the factorial version in a way that works for all nonnegative integer n?

[Walter Roberson]

I reported a bug (weakness) to Maplesoft about their option to notice special cases like these not having functioned properly. I will see if anyone happens to reply with a work around... though I don't think it likely I will get a reply, I admit.

[Fereshteh Emami Piri]

I understand, and will keep fingers crossed you will get a proper reply since it is very important to me. Please let me know if there is anything I should do. Anything that might help! In the meantime, thank you very much for all your helps and quick responses.

[Walter Roberson]

Unfortunately the reply that I did get pointed out that the Maple option I was reporting the problem on did not apply in this particular case.

The person who relied unfortunately failed to understand the importance of n being indefinite, so their suggestion was based around substituting in particular numeric n before doing the integration.

[Fereshteh Emami Piri]

OK, I can do that. However, the problem is I work with MATLAB, and when I use the code you suggested:

 n=1;
 syms v C sigma v0 nonnegative
pi = sym('pi');
f=((exp(-C*v))*(((C*v)^n)/factorial(n)))*(1/(sigma*(pi/2)^sym(1/2)*erfc(-v0/(sym(2)^sym(1/2)*sigma))))*(exp(-(((v-v0)^sym(2))/(sym(2)*sigma^sym(2)))));
simplify(int(f,v,[0 inf]))
simplify(subs(ans,[v0,sigma,C],[19,2,11]))
double(ans)

I don't get the same answer as yours, and it is very time consuming.

Answer 2

[Fereshteh Emami Piri]

Walter,

I appreciate your quick answer. In the attached file under Appendix B section, they explained how they ended up in the reported answer, which is equation (x). Actually, the authors derived the final equation manually. I use MATLAB.