Using vpaintegral to calculate integrations.

39 views (last 30 days)
Qian Feng
Qian Feng on 22 Sep 2017
Commented: Qian Feng on 27 Sep 2017
clear all
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];
leg = [];
for i = 0:d
l = [coeffs(legendreP(i,(2*y + sym(r1 + r2))/ sym(r3) )) zeros(1,d-i)];
leg = [leg; l];
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);
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 Comment
Walter Oevel
Walter Oevel on 26 Sep 2017
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.

Sign in to comment.

Accepted Answer

Walter Roberson
Walter Roberson on 25 Sep 2017
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.

Sign in to comment.

More Answers (1)

Nicolas Schmit
Nicolas Schmit on 25 Sep 2017
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.
Ga2 = vpaintegral(ga2,x,-r2,-r1, 'MaxFunctionCalls', Inf)*Dde*(r3^(-1));
Ga2 = zeros(size(ga2));
for k=1:numel(ga2)
Ga2(k) = integral(matlabFunction(ga2(k)),-r2,-r1);
Ga2 = Ga2*Dde*(r3^(-1));
Walter Roberson
Walter Roberson on 25 Sep 2017
Numeric processing of this with only 15 to 16 significant digits is not enough.

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!

Translated by