**You are now following this question**

- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.

# Integral inside integral is taking too much time

45 views (last 30 days)

Show older comments

Hi,

I am trying to evaluate integral inside integral using the following script but it is taking too much time, is there anyway I can reduce the time of evaluation:

f =@(x,y) x+y; % "x+y is just a simple example, the actual function is way more complicated"

integral(@(y) exp(-integral(@(x) f(x,y),a,b,'RelTol',1e-8,'AbsTol',1e-13),aa,bb,'RelTol',1e-8,'AbsTol',1e-13,'ArrayValued',true)

##### 5 Comments

Waseem Akhtar
on 20 Jul 2021

f is a function of both x and y, I'm not sure if analytical integration is possible for such a case. Can you give some example scriptt, please.

Thank you

Torsten
on 20 Jul 2021

Edited: Walter Roberson
on 20 Jul 2021

syms x y a b aa bb

f = x + y;

Fx = int(f,x,a,b);

Fxy = int(exp(-Fx),y,aa,bb)

Fxy =

integrates your function from above symbolically.

Waseem Akhtar
on 23 Jul 2021

Edited: Waseem Akhtar
on 23 Jul 2021

@Torsten Thank you for your reply. I tried symbolic integration but it is even taking longer than the numerical integration

Thank you for your help!

### Accepted Answer

Walter Roberson
on 20 Jul 2021

Use worse tolerances to reduce the computation time, at the expense of decreased accuracy.

In some cases, switching to integral2 might help. You are using arrayvalued so that would require multiple calls to integral2.

The fastest approach can depend on what your function is sensitive to. Numeric integration is adaptive, so if part of the range has rapid changes then it can slow down. In some cases, the rapid change only applies for one part of a range, and if you are using arrayvalued you can end up having to do the detailed integration for everything according to the worst part of the range, so sometimes arrayfun of non-arrayvalued can be faster. Sometimes .

##### 14 Comments

Waseem Akhtar
on 23 Jul 2021

@Walter Roberson Thank you for your reply. Actually, the integral I'm trying to evaluate is not a simple double integral rather, one of the integrals (the inner integral) is inside the expression for the second integral. My code looks like the following:

tic

aa=0; bb=Inf; a=0.000001; b=3.8e-04;

value =integral(@(l) exp(-(3.8e-04-l.*0.97).*2e+05).*(integral(@(omega) (exp(-omega-((2e+05.*sqrt((3.8e-04-l.*0.97).^2+(1.61e-05+l.*0.26).^2)).^2./(4*omega))))./(2*omega),aa,bb)),a,b,'ArrayValued',true);

toc

I'm not sure if integral2 can be used for such a case.

Thank you!

Walter Roberson
on 24 Jul 2021

Q = @(v) sym(v); %convert to rational

syms L omega AA BB A B real

assume(A<=B)

aa = Q(0); bb = Q(Inf); a = Q(0.000001); b = Q(38)*10^-5;

value = int(exp(-(b-L.*Q(0.97)).*Q(2e+05)).*(int((exp(-omega-((Q(2e+05).*sqrt((b-L.*Q(0.97)).^2+(Q(161)*10^-7+L.*Q(0.26)).^2)).^2./(4*omega))))./(2*omega),omega,AA,BB)),L,A,B)

value =

tic

double(vpa(subs(value,[AA,BB,A,B],[aa,bb,a,b])))

ans = 1.4551e-17

toc

Elapsed time is 1.090468 seconds.

valF = matlabFunction(value, 'vars', [AA,BB,A,B])

valF = function_handle with value:

@(AA,BB,A,B)integral(@(L)exp(L.*1.94e+5-7.6e+1).*integral(@(omega)exp(-omega-((L.*(9.7e+1./1.0e+2)-3.8e-4).^2.*4.0e+10+(L.*(1.3e+1./5.0e+1)+1.61e-5).^2.*4.0e+10)./(omega.*4.0))./(omega.*2.0),AA,BB),A,B)

tic

valF(double(aa),double(bb),double(a),double(b))

ans = 7.2568e-18

toc

Elapsed time is 0.070170 seconds.

If you can let aa and bb be constants for your purpose, then the inner integral has a definite form, so the computation can be faster.

Waseem Akhtar
on 24 Jul 2021

@Walter RobersonThank you for your detailed reply. I was wondering why the result/ans of value fn (ans = 1.4551e-17) is different than that of valF (ans = 7.2568e-18), aren't these the result of evaluation of same expression by two different ways?

I was also curious how can we represent the integrals in the following form in Matlab:

Thank you once again!

Walter Roberson
on 24 Jul 2021

I was also curious how can we represent the integrals in the following form in Matlab:

What you see in my comment above is actual MATLAB output; you can get it by using the LiveScript feature.

I was wondering why the result/ans of value fn (ans = 1.4551e-17) is different than that of valF (ans = 7.2568e-18)

When specific AA BB values are substituted into the int() expression, MATLAB would try to evaluate the inner integral to a closed form expression, which would succeed in this particular case, giving a theoretically exact solution for the inner integral. Then with specific A and B values substituted in, MATLAB would try to find a closed form solution for the outer integral, which would fail, and the int() would be returned symbolically. At that point, the vpa() would take effect, and would convert the integral from int() to vpaintegral() and then vpaintegral() would have a go at executing. vpaintegral() has its own set of default absolute and relative tolerances and works according to the current digits() setting, producing a symbolic solution. The double() around the vpa() then converts the symbolic floating point number into a double precision floating point number.

For the other expression, matlabFunction() applied to the symbolic expression with symbolic limits, converts the symbolic nested int() into nested integral() calls, converting all of the symbolic constants into double precision and combining them where possible, and returning back code for an anonymous function. In some cases, converting the symbolic constants into floating point and combining them can cause significant loss of precision

When valF() is invoked with specific number limits, the anonymous function is executed, so integral() is invoked; integral() has its own default absolute and relative tolerances; I believe that integral() also uses a different integration method than vpaintegral() does.

The anonymous function that is generated can have moved expressions around: symbolic expressions are re-arranged for symbolic processing according to what is mathematically equivalent, but order of operations is important for double precision functions, because of round-off error.

You are doing integral of terms that look as if they might have exp() of a notable negative number; if so then the order of additions could be fairly important: adding smallest to largest could result in the small values adding up to something large, but adding large to small could end up with the values of the small items being neglected. Consider for example the harmonic sequence: if you add 1/1 + 1/2 + 1/3 + 1/4 ... then around 1/50000 the value being added in double precision becomes too small to change the output... and yet mathematically 1/50000 + 1/50001 + ... 1/inf adds up to infinity.

It is therefore not unexpected that the two versions will produce different answers.

Waseem Akhtar
on 25 Jul 2021

Edited: Walter Roberson
on 25 Jul 2021

@Walter RobersonThank you for your detailed reply. Your comprehensive explanation helped alot to get to the core of these ideas. I do have one more questions though:

When I use my original expression as below, the result/ans and time is somewhat similar to that of valF case in your comments but I further need to reduce the time:

tic

aa=0; bb=Inf; a=0.000001; b=3.8e-04;

value =integral(@(l) exp(-(3.8e-04-l.*cos(0.27)).*2e+05).*(integral(@(omega) (exp(-omega-((2e+05.*sqrt((3.8e-04-l.*cos(0.27)).^2+(1.61e-05+l.*sin(0.27)).^2)).^2./(4*omega))))./(2*omega),aa,bb)),a,b,'ArrayValued',true);

toc

%value = 5.1023e-18

%Elapsed time is 0.106522 seconds.

Actually, I have to use a double loop to evaluate this integral at multiple locations as shown in my code below. This evaluation takes minutes to run which is too much. Basically, the inner integral is a modified form of 'besselk(0,z)' function in matlab. When I use the original 'besselk(0,z)' function, the whole calculation takes only secs to evaluate. Since, I must modify some parameters of this function, it seems that I don't have anyother choice but to use the integral form. I wonder if there is a way to significantly reduce the time of computation. The code is as under:

tic

j=0; %counter

for U=0.0001e-3:0.01e-3:0.2e-3

j=j+1;

i=0; %counter

for T=-1.14e-3:0.0057e-3:1.14e-3

i=i+1;

aa=0; bb=Inf; a=0.000001; b=3.8e-04;

value(i,j) =integral(@(l) exp(-(T-l.*cos(0.27)).*2e+05).*(integral(@(omega) (exp(-omega-((2e+05.*sqrt((T-l.*cos(0.27)).^2+(U+l.*sin(0.27)).^2)).^2./(4*omega))))./(2*omega),aa,bb)),a,b,'ArrayValued',true);

end

end

toc

THANK YOU!

Waseem Akhtar
on 25 Jul 2021

@Walter RobersonDid you edited a comment after my last comment? I see the status that you edited a comment but I don't see the comment.

Thank you

Walter Roberson
on 25 Jul 2021

If I recall correctly, I edited your comment to format your code as code blocks.

Waseem Akhtar
on 25 Jul 2021

@Walter Roberson I am trying to evaluate the following integral in a loop but it's giving error. Is it possible to avoid this error. The code is as under:

syms L omega AA BB A B U T

value = int(exp(-(T-L.*(0.97)).*(2e+05)).*(int((exp(-omega-(((2e+05).*sqrt((T-L.*(0.97)).^2+(U+L.*(0.26)).^2)).^2./(4*omega))))./(2*omega),omega,AA,BB)),L,A,B);

valF = matlabFunction(value, 'vars', [AA,BB,A,B,U,T]);

aa = 0; bb = Inf; a = 0.000001; b = 38*10^-5;

tic

j=0; %counter

for U=0.0001e-3:0.01e-3:0.2e-3

j=j+1;

i=0; %counter

for T=-1.14e-3:0.0057e-3:1.14e-3

i=i+1;

Qty(i,j) =valF(double(aa),double(bb),double(a),double(b),double(U),double(T));

end

end

toc

Thank you!

Walter Roberson
on 25 Jul 2021

tic

j=0; %counter

for U=0.0001e-3:0.01e-3:0.2e-3

j=j+1;

i=0; %counter

for T=-1.14e-3:0.0057e-3:1.14e-3

i=i+1;

aa=0; bb=Inf; a=0.000001; b=3.8e-04;

value(i,j) =integral(@(l) exp(-(T-l.*cos(0.27)).*2e+05).*(integral(@(omega) (exp(-omega-((2e+05.*sqrt((T-l.*cos(0.27)).^2+(U+l.*sin(0.27)).^2)).^2./(4*omega))))./(2*omega),aa,bb,'ArrayValued', true)),a,b,'ArrayValued',false);

end

end

toc

Elapsed time is 9.459138 seconds.

Walter Roberson
on 25 Jul 2021

Comparison:

Your code, with ArrayValued set on the outside integral but not the inside integral, took 151 seconds on my system.

My code, with ArrayValued set on the inside integral but not the outside integral, took 7 seconds on my system.

I am in the process of comparing the results to see if they are substantially the same.

Waseem Akhtar
on 26 Jul 2021

@Walter Roberson Bravo! That works wonders, Thank you!

Could you please explain a bit the logic behind setting the Arrayvalued "true" for inner integral and viceversa.

Also, could you please give a look to the last code I shared.

Thanks

Walter Roberson
on 14 Aug 2021

The outer integral() passes a vector of values into the inner interval. When the inner integral is set to ArrayValued, the inner integral only passes in one trial omega value at a time, and for that omega does the value calculation that is vectorized with respect to the vector of values from the outer integral, l . The inner integral continues to pass in one value at a time until it is satisfied that it has met integration tolerances for all elements corresponding to the positions in l .

When the outer integral is marked as array valued instead, but the inner integral is not, then the outer integral would pass a single value at a time to the inner integral. The inner integral would then construct a vector of values to pass to be calculated with the scalar l from the outer integral.

What is the difference, since either way the inner calculation is going to be evaluate with a vector? Well, there is not necessarily an inherent advantage to one or the other. However... depending on the calculations for each of the variables, it might happen to be the case that numerically one of the two versions happens to have its integral converge faster than the other way.

Walter Roberson
on 14 Aug 2021

value = int(exp(-(T-L.*(0.97)).*(2e+05)).*(int((exp(-omega-(((2e+05).*sqrt((T-L.*(0.97)).^2+(U+L.*(0.26)).^2)).^2./(4*omega))))./(2*omega),omega,AA,BB)),L,A,B);

valF = matlabFunction(value, 'vars', [AA,BB,A,B,U,T]);

That is going to fail because matlabFunction does not generate integral2() automatically when it sees nested integral(), and matlabFunction() does not know it should generate 'arrayvalued' at any point.

### More Answers (0)

### See Also

### Categories

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!**An Error Occurred**

Unable to complete the action because of changes made to the page. Reload the page to see its updated state.

Select a Web Site

Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .

You can also select a web site from the following list:

## How to Get Best Site Performance

Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.

### Americas

- América Latina (Español)
- Canada (English)
- United States (English)

### Europe

- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)

- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)

### Asia Pacific

- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)