Using vpaintegral to calculate integrations.

clear all
tic 

my = 2; 

d = 15; de = 15; 

r1 = 2;  r2 = 4.05; r3 = r2 - r1;

dai1 = d+1;  dai2 = de+1;  

ro = (dai1 + dai2)*my; 

syms y real 

le = [];
for i = 0:d
    l = [coeffs(legendreP(i,(2/sym(r1))*y+1))   zeros(1,d-i)];
le = [le; l];
end

leg = [];
for i = 0:d
    l = [coeffs(legendreP(i,(2*y + sym(r1 + r2))/ sym(r3)  ))   zeros(1,d-i)];
leg = [leg; l];
end

syms x real

xp = (x.^(0:d))';

l1 = le*xp; l2 = leg*xp; 

a = 18; 

fi1 = [sin(a*x); cos(a*x); exp(sin(a*x)); exp(cos(a*x))]; fi2 = fi1;

ga1 = fi1*l1'; ga2 = fi2*l2';

Dd  = diag(2.*(0:d)+1); Dde  = diag(2.*(0:de)+1);

Ga1 = vpaintegral(ga1,x,-r1,0, 'AbsTol',1e-20, 'RelTol',1e-20, 'MaxFunctionCalls', Inf)*Dd*(r1^(-1)); 

ep1 =   simplify((fi1 - Ga1*l1)*(fi1 - Ga1*l1)', 'Steps', 100);    

E1  = vpaintegral(ep1,x, -r1, 0, 'AbsTol',1e-20, 'RelTol',1e-20);   

Ga2 = vpaintegral(ga2,x,-r2,-r1, 'MaxFunctionCalls', Inf)*Dde*(r3^(-1));

ep2 = simplify((fi2 - Ga2*l2)*(fi2 - Ga2*l2)', 'Steps', 100); 

E2  = vpaintegral(ep2,x, -r2,-r1, 'MaxFunctionCalls', Inf);
toc
eta1 = 1; eta2 = 1;

I am trying to calculate the numerical values of Ga1, Ga2 and E1,E2 via the function vpaintegral. The functions inside of the integrations contain no singularities I believe, but it seems Ga2 demanding an very long time (I did not get the value of Ga2) to be calculated. If we put the breakpoint at Ga2, you can see that Ga1 and E1 can be calculated within a reasonable period. The reason I am puzzled here is that the type of functions among Ga1 Ga2, are the same, so I do not know why the calculation of Ga2 is so difficult.

Could somebody give me some suggestions? Thanks a lot!

[Walter Oevel]

The numerical evaluation of the higher degree Legendre polynomials is not very stable. When no AbsTol/RelTol values are specified in the call of vpaintegral, a working precision of only 8 decimal digits is used internally by the vpa evaluation in the symbolic engine. Unfortunately, with this low number of digits, the convergence of the integral is marred by the round-off in the evaluation of the Legendre polynomials .

I recommend to use the AbsTol/RelTol values that you specified for Ga1 for the Ga2 computation as well (note that you did not specify any tolerances for the computation of Ga2). Specifying the tolerances make the computation much faster.

Answer 1

[Walter Roberson]

This was also posted https://www.mathworks.com/matlabcentral/answers/311003-numerical-values-of-integrals#comment_486666. I replied there,

In your Ga2 computation, you left out the AbsTol and RelTol options. The ga2(end-1) is quite sensitive to the number of digits of computation near -4.05

Qian Feng then asked,

"Thanks Walter, so what number do you suggest that I should choose for the AbsTol and RelTol in this case?"

I replied,

My tests with a different package suggested that 20 or so digits of precision was needed to get a decent result, so 1e-20 like you used for Ga1 might be enough.

[Walter Roberson]

The Legendre polynomials are ill conditioned; https://www.hindawi.com/journals/jam/2014/742538/

[Qian Feng]

Thanks Walter, could do tell me what 'different package' you have applied to find the information in your previous post ?

[Walter Roberson]

I used Maple and pushed up the Digits until I got a reasonable looking graph of the function being integrated.

[Qian Feng]

Thanks Walter, so basically in this situation, we cannot expect Matlab can quickly produce a numerical result via vpaintegral. But based on you experiment Maple may get an acceptable result within a reasonable period, is that right?

[Walter Roberson]

Pass in the tolerance values like you do for ga1 and the result is fast.

Ga2 = vpaintegral(ga2,x,-r2,-r1, 'AbsTol',1e-20, 'RelTol',1e-20, 'MaxFunctionCalls', Inf)*Dde*(r3^(-1));

Or

 NumDig = 20;
 RR = int(ga2,x,-r1,r1)*Dde*(r3^(-1));    %it is able to do some of them exactly
 GA2 = vpa(RR, NumDig);

The difference between these two approaches appears to be on the order of 1E-16

[Qian Feng]

Thanks Walter. Are you sure about the command Ga2 = vpaintegral(ga2,x,-r2,-r1, 'AbsTol',1e-20, 'RelTol',1e-20, 'MaxFunctionCalls', Inf)*Dde*(r3^(-1)); can produce result in a reasonable period ? I think I had tried that previously but it seems very slow to get the result.

[Walter Roberson]

Checking again I see that the fast integration was when I used -r1,r1 for the bounds instead of -r2,-r1 . Even -r2,r1 runs quickly. One would expect that integrating over a larger range would take longer, but I guess with the larger range it takes larger steps and does not notice the integration difficulties.

I will need to do more checking.

[Qian Feng]

Thanks Walter.

Answer 2

[Nicolas Schmit]

Since you are looking for a numerical solution, a workaround to speed up the calculation is to convert the symbolic function into a MATLAB function then use a numerical integration routine to calculate the integral.

Replace

Ga2 = vpaintegral(ga2,x,-r2,-r1, 'MaxFunctionCalls', Inf)*Dde*(r3^(-1));

With

Ga2 = zeros(size(ga2));
for k=1:numel(ga2)
    Ga2(k) = integral(matlabFunction(ga2(k)),-r2,-r1);
end
Ga2 =  Ga2*Dde*(r3^(-1));

[Qian Feng]

Hi Nicolas, thanks for the answer. However, I got the following warning concerning some integrals I believe after I run the program.

Reached the limit on the maximum number of intervals in use. Approximate bound on error is 7.0e-02. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.

The Approximate bound error 7.0e-02 may not be accurate enough for my problem. That is the reason why I appeal to the usage of the function vpaintegral

[Walter Roberson]

Numeric processing of this with only 15 to 16 significant digits is not enough.