Using vpaintegral to calculate integrations.

Question

Qian Feng 2017-9-22

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/357855-using-vpaintegral-to-calculate-integrations

评论： Qian Feng 2017-9-27

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!

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Walter Oevel 2017-9-26

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 2017-9-25

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/357855-using-vpaintegral-to-calculate-integrations#answer_282960

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.

8 个评论
显示 6更早的评论隐藏 6更早的评论

Walter Roberson 2017-9-26

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 2017-9-27

Thanks Walter.

请先登录，再进行评论。

Answer 2

Nicolas Schmit 2017-9-25

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/357855-using-vpaintegral-to-calculate-integrations#answer_282953

在 MATLAB Online 中打开

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));

2 个评论
显示无隐藏无

Qian Feng 2017-9-25

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