9 views (last 30 days)

Hi,

The function "funct" takes up almost 0.06 seconds and it is the first responsible for a slow execution of an integral (more than 3 hours) when it is used in it with regards to y.

Any help please to reduce the computation time of the following function

function [out] = funct(mu,lmax,wmax,hmax,lambdab,hu,hb,y,x)

nx=length(x);

out=zeros(1,nx);

A=(mu./pi).*(wmax+lmax);

B=(mu.*wmax.*lmax)./4;

C1=(hb-(2.*hmax))./(2*hmax);

C2=(hu+hb-(2.*hmax))./(2*hmax);

for i=1:nx

S1=@(u,theta) u.*(1-exp(C1.*((A.*sqrt(u.^2+x(i).^2-(2.*x(i).*u.*cos(theta))))+B)));

S2=@(theta) (1-exp(C1.*((A.*sqrt(y.^2+x(i).^2-(2.*x(i).*y.*cos(theta))))+B)));

out(i)=lambdab.*y.*integral(S2,0,2*pi,'ArrayValued',true).*exp(-lambdab.*quad2d(S1,0,y,0,2*pi))./(1-exp(-lambdab.*quad2d(S1,0,1000,0,2*pi)));

end

end

Walter Roberson
on 25 May 2020

By going back and replacing that 1000 limit with infinity, and adding the assumption that all values are non-negative, and that hb < 2*hmax, Maple is able to combine terms to remove one integral. The symbolic form is

-exp(lambdab * int( int( u * (exp(1/2 * (1/4 * lmax * wmax * pi + (lmax+wmax) * ...

(u^2 - 2 * cos(theta) * u * x + x^2)^(1/2)) * mu * (hb-2*hmax) / pi / hmax)-1), u, 0, y), ...

theta, 0, 2*pi)) * lambdab * y * int(-1 + exp(1/8 * ((4 * lmax + 4 * wmax) * ...

(x^2 - 2 * cos(theta) * x * y + y^2)^(1/2) + lmax * wmax * pi) * (hb-2*hmax) * mu / pi / hmax), ...

theta, 0, 2*pi)

This should be more efficient

Walter Roberson
on 25 May 2020

Walter Roberson
on 25 May 2020

per isakson
on 24 May 2020

Edited: per isakson
on 25 May 2020

On R2018b and a new vanilla PC I found a significant decrease in computation time by replacing quad2d() by integral2() (with the same signature, i.e. defaults are used). With your input data this replacement didn't affect the result, out. (Caveat, I do this because I'm curious, not because I know much about the topic.) There are two posts on Double Integration in MATLAB ( and here ) at Loren on the Art of MATLAB. I don't find backing for the difference in speed that I've seen, which makes me a bit nervous. What's the catch? There ought to exist a discussion on how to chose between the two functions, but I fail to find it.

Have you tried to simplify the expressions with the help of the Symbolic Toolbox? See https://se.mathworks.com/matlabcentral/answers/33399-problem-with-double-integral-dblquad-and-quad2d#answer_42126 it's a bit old but shows the idea. Maybe that toolbox sees something that I don't see.

In response to comments

The difference in speed between integral2(), introduced in R2012a, and quad2d(), introduced in R2009a, that I see is irking.

I've compared your function, funct(), to funct_int2(), in which I replaced quad2d() by integral2(). (And split the long statement into three to make the code and the result of profile() easier to read.)

>> cssm % first call to cssm after starting Matlab

Elapsed time is 0.364935 seconds.

Elapsed time is 0.140186 seconds.

0

>> cssm

Elapsed time is 0.104307 seconds.

Elapsed time is 0.021930 seconds.

0

>> cssm

Elapsed time is 0.088423 seconds.

Elapsed time is 0.016356 seconds.

0

...

...

>> cssm

Elapsed time is 0.091209 seconds.

Elapsed time is 0.017180 seconds.

0

where cssm.m is a script containing

%%

mu=1; lmax=4; wmax=6; hmax=30; lambdab=1; hu=40; hb=10; y=4; x=2;

%%

tic

for jj = 1 : 10

out1 = funct(mu,lmax,wmax,hmax,lambdab,hu,hb,y,x);

end

toc

%%

tic

for jj = 1 : 10

out2 = funct_int2(mu,lmax,wmax,hmax,lambdab,hu,hb,y,x);

end

toc

disp( out2-out1 )

and

function [out] = funct_int2(mu,lmax,wmax,hmax,lambdab,hu,hb,y,x)

nx=length(x);

out=zeros(1,nx);

A=(mu./pi).*(wmax+lmax);

B=(mu.*wmax.*lmax)./4;

C1=(hb-(2.*hmax))./(2*hmax);

C2=(hu+hb-(2.*hmax))./(2*hmax); %#ok<NASGU>

for ii=1:nx

S1=@(u,theta) u.*(1-exp(C1.*((A.*sqrt(u.^2+x(ii).^2-(2.*x(ii).*u.*cos(theta))))+B)));

S2=@(theta) (1-exp(C1.*((A.*sqrt(y.^2+x(ii).^2-(2.*x(ii).*y.*cos(theta))))+B)));

f1 = lambdab.*y.*integral(S2,0,2*pi,'ArrayValued',true);

f2 = exp(-lambdab.*integral2(S1,0,y,0,2*pi));

f3 = (1-exp(-lambdab.*integral2(S1,0,1000,0,2*pi)));

out(ii) = f1.*f2./f3;

end

end

Comments

- Fine print; don't blame me for mistakes.
- In this comparison both integral2() and quad2d() use the "tiled method"
- The fact that the two return exactly the same result makes me believe that deep inside the same calculation is made.

per isakson
on 25 May 2020

Opportunities for recent engineering grads.

Apply Today
## 4 Comments

## Direct link to this comment

https://fr.mathworks.com/matlabcentral/answers/531213-how-to-reduce-computation-time#comment_859878

⋮## Direct link to this comment

https://fr.mathworks.com/matlabcentral/answers/531213-how-to-reduce-computation-time#comment_859878

## Direct link to this comment

https://fr.mathworks.com/matlabcentral/answers/531213-how-to-reduce-computation-time#comment_859893

⋮## Direct link to this comment

https://fr.mathworks.com/matlabcentral/answers/531213-how-to-reduce-computation-time#comment_859893

## Direct link to this comment

https://fr.mathworks.com/matlabcentral/answers/531213-how-to-reduce-computation-time#comment_861523

⋮## Direct link to this comment

https://fr.mathworks.com/matlabcentral/answers/531213-how-to-reduce-computation-time#comment_861523

## Direct link to this comment

https://fr.mathworks.com/matlabcentral/answers/531213-how-to-reduce-computation-time#comment_862028

⋮## Direct link to this comment

https://fr.mathworks.com/matlabcentral/answers/531213-how-to-reduce-computation-time#comment_862028

Sign in to comment.