The computer hangs while evaluating the following integral and code....pls help

w=15;a=100;b=0.01;OUT=1/30;deltaOUT=5.2;
syms x;
pOUT=@(x)max(0,(1-exp(-OUT*x+deltaOUT)));
pNL=@(x)1-pOUT(x);

err=45.*2.*pi./180;mu=err;sg=1;
pdnormal = makedist('Normal',mu,sg);pd= truncate(pdnormal,0,inf);u0=cdf(pd,err);u1=1-u0;

syms y z t;
f1=u1.*dirac(y-b)+u0.*dirac(y-a);f2=u1.*dirac(z-b)+u0.*dirac(z-a);
fGTi=(w/(2*pi)*(a))+(1-(w/(2*pi))*b);fGRi=(w/(2*pi)*(a))+(1-(w/(2*pi))*b);

alpha=2.1;beta=4;d=10;neta=10.^(-12);Pt=1;lambda = 10.^(-3);product=-2*pi*lambda;

intfn=(1-(1/(1+((beta*d^(alpha)*fGTi*fGRi*t.^(-alpha))/(Pt*y*z))))).*t.*product;
innerintegral=int(intfn,t,0,1);
mainfun=exp(innerintegral)*pNL(d)*exp((-beta*d^(alpha)*neta)/(Pt*y*z))*f1*f2;
resultintermediate=int(mainfun,y,0.01,inf);
result=int(resultintermediate,z,0.01,inf);
double(result) 

[Amit dot207]

It is pls. help....not held

[Greg Heath]

Thanks for the correction:

My arms are so, so tired!!

Greg

[Amit dot207]

That was a witty comment, Greg and quite entertaining. However, MATLAB requires superior wits, as you would have discovered by now trying to solve the problem.

[Jan]

@Amit: You can edit your question directly instead of appending comments to fix typos.

What does "computer hangs" mean? Do you get an error message? Does the computation take longer than you expect? How much longer and what do you expect? Does Matlab crash or the computer melt down?

When I run your code, the line:

innerintegral=int(intfn,t,0,1);

takes a long time. It is useful, if such details are mentioned in the question directly. Do not let the readers find out, where the problem is, if you can provide this information already.

By the way, your code would be much easier to read, if you write one command per line and insert spaces. A standard scheme is: Spaces around all operators and one space after a comma.

Strange: fGTi and fGRi are exactly the same:

fGTi=(w/(2*pi)*(a))+(1-(w/(2*pi))*b);
fGRi=(w/(2*pi)*(a))+(1-(w/(2*pi))*b);

This looks simpler:

fGTi = 1 + (a - b) * w / (2 * pi);

This will not influence the processing speed remarkably, but teh debugging gets much easier.

While neta=10.^(-12) is a very expensive power operation, neta=1e-12 is a cheap constant.

[Amit dot207]

Thanks, Jan for answering. There is no error report and the status bar keeps on saying "Busy". I tried out the exponential thing also and the formula which you suggested but result is the same.

[Jan]

If you omit the y and z in the function to be integrated, you get a rather ugly and huge result already. A simplification of the symbolic formula does not decrease the runtime also.

[Amit dot207]

Yes. Is there a way around it?

[Jan]

A way around what? Do y and z depend on t? Is there a need to solve the equation symbolically or would a numerical solution match your needs - at least you request a numerical output finally.

It is hard to test for me. I do not have the Symbolic Toolbox, so I run you code online. Unfortunately I cannot break the execution there, such that the window is busy until a timeout occurs. Each trial costs 5 or 10 minutes and my experiments are tedious. What did you try so far? maybe some new insights? Can you post the output of:

simplify(intfn)

? Maybe some other users understand, why the integration is such expensive.

[Amit dot207]

simplify(intfn) gives

-(7767242930207539*pi*t)/(500*(268435456*t^(21/10)*y*z + 7767242930207539))

I have tried out for half an hour but the computer says only busy

[Amit dot207]

changed the code to the following:

syms x;
pOUT=@(x)max(0,(1-exp(-OUT*x+deltaOUT)));
pNL=@(x)1-pOUT(x);

err=45.*2.*pi./180;mu=err;sg=1;
pdnormal = makedist('Normal',mu,sg);pd= truncate(pdnormal,0,inf);u0=cdf(pd,err);u1=1-u0;

syms y z t;
f1=@(y)u1.*dirac(y-b)+u0.*dirac(y-a);f2=@(z)u1.*dirac(z-b)+u0.*dirac(z-a);
fGTi=(w./(2.*pi).*(a))+(1-(w./(2.*pi)).*b);fGRi=(w./(2.*pi).*(a))+(1-(w./(2.*pi)).*b);

alpha=2.1;beta=4;d=10;neta=1e-12;Pt=1;lambda = 1e-3;product=-2*pi*lambda;

intfn=@(y,z,t)(1-(1./(1+((beta.*d.^(alpha).*fGTi.*fGRi.*t.^(-alpha))./(Pt.*y.*z))))).*t.*product;
innerintegral=@(y,z)integral(@(t)intfn(y,z,t),0,1);
mainfun=@(y,z)simplify(exp(innerintegral(y,z)).*pNL(d).*exp((-beta*d.^(alpha).*neta)./(Pt.*y.*z)).*f1(y).*f2(z));
result=integral2(@(y,z)mainfun(y,z),0.01,inf,0.01,inf);
double(result) 

[Amit dot207]

The above code shows the following error:

Error using ./

Matrix dimensions must agree.

Error in @(y,z,t)(1-(1./(1+((beta.*d.^(alpha).*fGTi.*fGRi.*t.^(-alpha))./(Pt.*y.*z))))).*t.*product

Error in @(t)intfn(y,z,t)

Error in integralCalc/iterateScalarValued (line 314)

fx = FUN(t);

Error in integralCalc/vadapt (line 132)

[q,errbnd] = iterateScalarValued(u,tinterval,pathlen);

Error in integralCalc (line 75)

[q,errbnd] = vadapt(@AtoBInvTransform,interval);

Error in integral (line 88)

Q = integralCalc(fun,a,b,opstruct);

Error in @(y,z)integral(@(t)intfn(y,z,t),0,1)

Error in @(y,z)simplify(exp(innerintegral(y,z)).*pNL(d).*exp((-beta*d.^(alpha).*neta)./(Pt.*y.*z)).*f1(y).*f2(z))

Error in @(y,z)mainfun(y,z)

Error in integral2Calc>@(y)fun(xi*ones(size(y)),y) (line 18)

@(y)fun(xi*ones(size(y)),y),y1i,y2i,opstruct.integralOptions), ...

Error in integralCalc/iterateScalarValued (line 314)

fx = FUN(t);

Error in integralCalc/vadapt (line 132)

[q,errbnd] = iterateScalarValued(u,tinterval,pathlen);

Error in integralCalc (line 83)

[q,errbnd] = vadapt(@AToInfInvTransform,interval);

Error in integral2Calc>@(xi,y1i,y2i)integralCalc(@(y)fun(xi*ones(size(y)),y),y1i,y2i,opstruct.integralOptions) (line 17)

innerintegral = @(x)arrayfun(@(xi,y1i,y2i)integralCalc( ...

Error in integral2Calc>@(x)arrayfun(@(xi,y1i,y2i)integralCalc(@(y)fun(xi*ones(size(y)),y),y1i,y2i,opstruct.integralOptions),x,ymin(x),ymax(x))

(line 17)

innerintegral = @(x)arrayfun(@(xi,y1i,y2i)integralCalc( ...

Error in integralCalc/iterateScalarValued (line 314)

fx = FUN(t);

Error in integralCalc/vadapt (line 132)

[q,errbnd] = iterateScalarValued(u,tinterval,pathlen);

Error in integralCalc (line 83)

[q,errbnd] = vadapt(@AToInfInvTransform,interval);

Error in integral2Calc>integral2i (line 20)

[q,errbnd] = integralCalc(innerintegral,xmin,xmax,opstruct.integralOptions);

Error in integral2Calc (line 7)

[q,errbnd] = integral2i(fun,xmin,xmax,ymin,ymax,optionstruct);

Error in integral2 (line 106)

Q = integral2Calc(fun,xmin,xmax,yminfun,ymaxfun,opstruct);

Error in mymatlab_script_10_1_2019 (line 18)

result=integral2(@(y,z)mainfun(y,z),0.01,inf,0.01,inf);

[Walter Roberson]

Your line

fGTi=(w./(2.*pi).*(a))+(1-(w./(2.*pi)).*b);fGRi=(w./(2.*pi).*(a))+(1-(w./(2.*pi)).*b);

uses w without it having been defined.

[Amit dot207]

w=15; defined in the first line

[Walter Roberson]

Please recheck your line

f1=u1.*dirac(y-b)+u0.*dirac(y-a);f2=u1.*dirac(z-b)+u0.*dirac(z-a);

That defines f1 as u1 if y == b, and as u0 if y == a, and as 0 otherwise, and f2 does the same except with z as the variable. Therefore f1 and f2 are non-zero only when y == a exactly or y == b exactly, or likewise in z.

Then in mainfun, you multiply the overall result by f1 and f2. The effect of that is to have mainfun be 0 exactly when y == a and z == a, or y == a and z == b, or y == b and z == a, or y == b and z == b -- so 0 exactly at four exact points. The integral of that over values that include a and b is going to be the the sum of the formula evaluated only at those four exact combinations. And if that is the intention then you would save yourself a lot of integration effort to just do the sums directly.

[Amit dot207]

Agreed. In that case also can you pls suggest the modification in the code so that the result is obtained?

[Walter Roberson]

Is it correct that you only want the sum of those four combinations? Or is your dirac in error ? Perhaps you want a heaviside combination, in order to create a window ?

[Amit dot207]

I want only a sum of those four combinations. That is how the formula goes. A window is not the idea. It is only a chance that the formula has come like this.

[Walter Roberson]

resultintermediate=int(mainfun,y,0.01,inf);
result=int(resultintermediate,z,0.01,inf);

would be

result = subs(mainfun, [y, z], [a, a]) + subs(mainfun, [y z], [a, b]) + subs(mainfun, [y z], [b a]) + subs(mainfun, [y z], [b b])

[Amit dot207]

Accepted. Thanks a lot for your time spent on the matter and deeply grateful to you. I will request you to give a suggestion on how to put a waitbar on the original integration problem.

[Walter Roberson]

The original integration problem will give you poor results. The symbolic toolbox is not well adapted to noticing the conditions under which the dirac is non-zero, so it calculates 0 as the answer. If you were to switch to numeric integration you would need to add waypoints at a and b, but waypoints are only available for 1D integration.

Symbolic integration by int() and by vpaintegral() have no built-in hooks for progress bars, and because they are adaptive they cannot easily estimate how much longer it will take to execute. Most of the numeric integrators are adaptive as well and so cannot easily estimate remaining time.

About the closest you could get for symbolic integration would be to use the Parallel Computing Toolbox in order to run the computation in a worker while the client ran a timer that updated a "busy" bar rather than a progress bar; even then it would have to pretty much guess that useful progress was being made at the symbolic engine level.

With numeric integration you could move one of the computations from anonymous function into a full function that updated a busy bar. The difference would be that under this circumstance that you would know that some real work was being done, whereas with the symbolic case you just have to trust that the symbolic engine is doing something useful.

[Walter Roberson]

Note: the symbolic engine is better at processing piecewise() than dirac(), but potentially much slower.