Using vpaintegral to calculate integrations.

9 次查看（过去 30 天）
Qian Feng
Qian Feng 2017-9-22
评论： 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 个评论
Walter Oevel
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.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Walter Roberson
Walter Roberson 2017-9-25
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更早的评论
Walter Roberson
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.

请先登录，再进行评论。

更多回答（1 个）

Nicolas Schmit
Nicolas Schmit 2017-9-25
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
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
Walter Roberson
Walter Roberson 2017-9-25
Numeric processing of this with only 15 to 16 significant digits is not enough.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Calculus 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by