How to solve for limit of double integral

I am new to MATLAB. Please help me. I want to solve for the limit of double integral.

syms u x
a = 5;
b = 0.2;
c = 10;
d = 0.1;
k=0.01;
m = 1./(2.^(a+c-1).*gamma(a).*gamma(c).*b.^a.*d.^c);
h(u,x) = m.*(u+x).^(a-1).*(u-x).^(c-1).*exp(-(u+x)./(2.*b)-(u-x)./(2.*d));
f(x) = int(h,u,abs(x), inf);
solPk = fzero(@(Pk) integral(f,-inf,Pk)-k,[-10,10]);

Thank you so much.

Answer 1

[Walter Roberson]

solPk = vpasolve(int(f,x,-inf,Pk) == k)

There is a closed form integral, but MATLAB is not able to find it. It is

piecewise(Pk < 0, -(770558495002939/5159780352000000000000)*exp(10*Pk)*(75937500*Pk^9-296156250*Pk^8+601425000*Pk^7-817897500*Pk^6+808258500*Pk^5-594641250*Pk^4+322528500*Pk^3-123369750*Pk^2+29996190*Pk-3512131), ...
          0 <= Pk, 1310719999999999239/1310720000000000000+(1/25798901760000000000000)*(-11650844444444437680000*Pk^4-40389594074074050624000*Pk^3-58409566814814780902400*Pk^2-41590925590123432642560*Pk-12267389871934149256192)*exp(-5*Pk))

You could fzero() or fsolve() on that after making a guess about whether Pk will be positive or negative.

[Mikie]

So, we are solving this problem with symbolic, right? I have a question. Can we solve this by numeric computation? How?

Thank you.

[Walter Roberson]

syms u x
a = 5;
b = 0.2;
c = 10;
d = 0.1;
k=0.01;
m = 1./(2.^(a+c-1).*gamma(a).*gamma(c).*b.^a.*d.^c);
h(u,x) = m.*(u+x).^(a-1).*(u-x).^(c-1).*exp(-(u+x)./(2.*b)-(u-x)./(2.*d));
f = matlabFunction(int(h,u,abs(x), inf));
solPk = fzero(@(Pk) integral(f,-inf,Pk)-k,[-10,10]);

It looks to me as if the function it produces is not correct for positive values, but it happens that for your purposes you do not need positive values.

[Mikie]

Thank you so much. Anyway, I do not need positive values.

[Mikie]

I have one more question. Why when

a = 30;
b = 0.2;
c = 20;
d = 0.1;

I cannot solve for solPk? It tells me that

Warning: Infinite or Not-a-Number value encountered. 
> In integralCalc/iterateScalarValued (line 349)
  In integralCalc/vadapt (line 132)
  In integralCalc (line 91)
  In integral (line 88)
  In @(Pk)integral(f,-inf,Pk)-k
  In fzero (line 230)
  In Percentile3FnNSH (line 42) 

    Warning: Infinite or Not-a-Number value encountered. 
    > In integralCalc/iterateScalarValued (line 349)
      In integralCalc/vadapt (line 132)
      In integralCalc (line 91)
      In integral (line 88)
      In @(Pk)integral(f,-inf,Pk)-k
      In fzero (line 240)
      In Percentile3FnNSH (line 42) 
    Error using fzero (line 242)
    Function values at interval endpoints must be finite and real.

    Error in Percentile3FnNSH (line 42)
    solPk = fzero(@(Pk) integral(f,-inf,Pk)-k,[-10,10])

[Walter Roberson]

The solution for that system is positive; it is possible that the integral is getting confused.

I would suggest that you split the integral at 0:

syms u x pk
a = 30;
b = 0.2;
c = 20;
d = 0.1;
k = 0.01;
m = 1./(2.^(a+c-1).*gamma(a).*gamma(c).*b.^a.*d.^c);
h(u,x) = m.*(u+x).^(a-1).*(u-x).^(c-1).*exp(-(u+x)./(2.*b)-(u-x)./(2.*d));
syms xm negative
syms xp positive
fm = int(h(u,xm),u,-xm, inf);  %xm is negative so -xm is abs(xm)
intm = int(fm, xm, -inf, xm);
int0 = subs(intm, xm, 0);
if int0 > k
  solpk = vpasolve(intm - k);
else
  fp = int(h(u,xp),u,xp,inf); %xp is positive so xp is abs(xp)
  intp = int(fp, xp, 0, xp);
  solpk = vpasolve(intp + int0 == k);
end

This should in theory be able to handle both cases, provided that MATLAB can get the integral right given the hint about whether x is positive or negative.

[Mikie]

Sorry, I mean the solution can be either positive or negative. I am so sorry. It's my bad. But why when I run the code didn't it give me any solution in case of a = 30; b = 0.2; c = 20; d = 0.1.

[Walter Roberson]

It gave me about 1.4 when I ran it.

[Mikie]

Could you give me any suggestion how to deal with this? This following happens when a = 30; b = 0.2; c = 20; d = 0.1 and when a = 25; b = 0.2; c = 40; d = 0.1. Thank you.

Warning: Infinite or Not-a-Number value encountered. 
> In integralCalc/iterateScalarValued (line 349)
  In integralCalc/vadapt (line 132)
  In integralCalc (line 91)
  In integral (line 88)
  In Percentile3FnNSH>@(Pk)integral(f,-inf,Pk)-k
  In fzero (line 228)
  In Percentile3FnNSH (line 30) 
Warning: Infinite or Not-a-Number value encountered. 
> In integralCalc/iterateScalarValued (line 349)
  In integralCalc/vadapt (line 132)
  In integralCalc (line 91)
  In integral (line 88)
  In Percentile3FnNSH>@(Pk)integral(f,-inf,Pk)-k
  In fzero (line 238)
  In Percentile3FnNSH (line 30) 
Error using fzero (line 240)
Function values at interval endpoints must be finite and real.

Error in Percentile3FnNSH (line 30)
solPk = fzero(@(Pk) integral(f,-inf,Pk)-k,[-10,10])

[Mikie]

Also, there is another problem when

a = 2.5
b = 0.2
c = 2
d = 0.1

I am wondering why it didn't work for some combination of a,b,c, and d. Thank you so much.

[Walter Roberson]

Don't do that. Split the integral at 0 like I advised.

syms u x pk
a = 30;
b = 0.2;
c = 20;
d = 0.1;
k = 0.01;
m = 1./(2.^(a+c-1).*gamma(a).*gamma(c).*b.^a.*d.^c);
h(u,x) = m.*(u+x).^(a-1).*(u-x).^(c-1).*exp(-(u+x)./(2.*b)-(u-x)./(2.*d));
syms xm negative
syms xp positive
fm = int(h(u,xm),u,-xm, inf);
intm = int(fm, xm, -inf, xm);
int0 = subs(intm, xm, 0);
if int0 > k
    disp('pk is negative, resolving it');
    solpk = vpasolve(intm - k);
    solpkn = fzero( matlabFunction(intm-k), [-10 0]);
else
    disp('pk is positive, resolving it');
    fp = int(h(u,xp),u,xp,inf);
    intp = int(fp, xp, 0, xp);
    solpk = vpasolve(intp + int0 == k);
    solpkn = fzero( matlabFunction(intp + int0 - k), [0 10]);
end
display(solpk)
display(solpkn)

[Mikie]

Oh Thank you so much. This is help me a lot. I still have another problem dealing with the integral. Could you please help me?

syms u x pk
k = 0.99
a = 5
b = 0.2
c = 40
d = 0.1
D1 = ((a.*b.^2+c.*d.^2).^3)/((a.*b.^3-c.*d.^3).^2);
D = (1/(2.*pi)).*((4./pi)-1)^2.*D1;
syms t
fun = (D+(2./pi).^3).*t.^3 - (3.*(2./pi).^2).*t.^2 + 3.*(2./pi).*t - 1;
theta = double(vpasolve(fun,t,[-Inf Inf]));
lambda = sign(a.*b.^3 - c.*d.^3).*sqrt(theta./(1-theta)); %sign
sigma = sqrt((a.*b.^2+c.*d.^2)./(1-(2./pi).*theta));
mu = a.*b-c.*d-sigma.*sign(a.*b.^3 - c.*d.^3).*sqrt(2./pi.*theta);
p = matlabFunction(2./sigma.*normpdf((x-mu)./sigma).*normcdf(lambda.*(x-mu)./sigma));
solPkS = fzero(@(PkS) integral(p,-Inf,PkS)-k,[-10,10]);

Here, I cannot solve for solPkS. It said that

Error using fzero (line 246)
FZERO cannot continue because user-supplied function_handle ==> @(PkS)integral(p,-Inf,PkS)-k failed with the error below.

   Too many input arguments.

Error in Percentile3FnNSH (line 56)
solPkS = fzero(@(PkS) integral(p,-Inf,PkS)-k,[-10,10]);

Your help is greatly appreciated.

[Walter Roberson]

The denominator of D1 = ((a.*b.^2+c.*d.^2).^3)/((a.*b.^3-c.*d.^3).^2); is 0, so D1 is inf, leading to empty theta, leading to matlabFunction being applied to a constant empty array. When you matlabFunction something that is a constant and you do not specify the 'vars' then it generates a function handle with no arguments.

[Mikie]

Thank you so much. I still have another problems. Is it possible for these cases:

a=40; b=0.1;c=60;d=0.0667;

and

a=30;b=0.1;c=75;d=0.0667;

How can I fix this? Thank you.

[Mikie]

Also, when

a = 50
b = 0.1
c = 60
d = 0.0667

Why does not it work?

[Walter Roberson]

With the old equations or with the D1 and D equations?

[Mikie]

For finding solpk and solPkS

[Mikie]

a=40;b=0.1;c=60;d=0.0667;k=0.99;
m = 1./(2.^(a+c-1).*gamma(a).*gamma(c).*b.^a.*d.^c);
h(u,x) = m.*(u+x).^(a-1).*(u-x).^(c-1).*exp(-(u+x)./(2.*b)-(u-x)./(2.*d));
syms xm negative
syms xp positive
fm = int(h(u,xm),u,-xm, inf);
intm = int(fm, xm, -inf, xm);
int0 = subs(intm, xm, 0);
if int0 > k
    disp('pk is negative, resolving it');
    solpk = vpasolve(intm - k);
    solpkn = fzero(matlabFunction(intm-k), [-10 0]);
else
    disp('pk is positive, resolving it');
    fp = int(h(u,xp),u,xp,inf);
    intp = int(fp, xp, 0, xp);
    solpk = vpasolve(intp + int0 == k);
    solpkn = fzero(matlabFunction(intp + int0 - k), [0 10]);
end

D1 = ((a.*b.^2+c.*d.^2).^3)/((a.*b.^3-c.*d.^3).^2);
D = (1/(2.*pi)).*((4./pi)-1)^2.*D1;
syms t
fun = (D+(2./pi).^3).*t.^3 - (3.*(2./pi).^2).*t.^2 + 3.*(2./pi).*t - 1;
theta = double(vpasolve(fun,t,[-Inf Inf]));
lambda = sign(a.*b.^3 - c.*d.^3).*sqrt(theta./(1-theta)); %sign
sigma = sqrt((a.*b.^2+c.*d.^2)./(1-(2./pi).*theta));
mu = a.*b-c.*d-sigma.*sign(a.*b.^3 - c.*d.^3).*sqrt(2./pi.*theta);
p = matlabFunction(2./sigma.*normpdf((x-mu)./sigma).*normcdf(lambda.*(x-mu)./sigma));
solPkS = fzero(@(PkS) integral(p,-Inf,PkS)-k,[-10,10]);

display(solpk)
display(solpkn)
display(solPkS)

[Walter Roberson]

The integral fm involves pretty much every combination of u and x with total power not exceeding (a+c-1), so in this case 40+60-1 = 99. That is 99*98/2 terms for that portion. It takes MATLAB a while to compute... and then it has to do the integral over that for fm. Takes a while :(

[Mikie]

But the process is stopped:

Error using symengine
Out of memory.

Error in sym/int (line 162)
   rSym = mupadmex('symobj::intdef',f.s,x.s,a.s,b.s,options);

Error in FnNSH (line 29)
intm = int(fm, xm, -inf, xm);

How can I fix it?

[Walter Roberson]

Ah, that does not surprise me, Now that you have mentioned it, I see that MuPAD silently died underneath me; I have 32 gigabytes of RAM but it must have filled it up and croaked.

You are probably going to need to switch to a more numeric path.

[Walter Roberson]

Bleh. The int() that it creates for fm involves a limit() as u approaches infinity, which is an operation that effectively cannot be converted using matlabFunction :(

[Mikie]

Oh what should I do then? :'(

[Walter Roberson]

I am still experimenting. It's pretty slow :(

[Mikie]

Thank you so much. :))

[Walter Roberson]

You might be interested to know that c = 56 does not produce the limit() inside fm, but that c = 57 does.