Numerical values of integrals

I am dealing with a problem of finding accurate numerical values of integrals. Specifically, the integral is introduced by using the best approximation scheme (Legendre Polynomials) to approximate a vector valued function whose indefinite integral is not easy to be explicitly written down. The code is provided as follows:

 r1 = 0.7;  r2 = 1; r3 = r2 - r1; d = 8;

 syms y real 

 le = [];

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

 leg = [];

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

 syms x
 t = [];

 for i = 0 : d
    t = [t ; x^i];     
 end
 xp = t;

 lp = le*xp;  la = leg*xp; 

 fi1 = [exp(sin(x)); exp(cos(x))]; fi2 = [sin(x^2); cos(x^2)];

 ny1 = size(fi1,1); ny2 = size(fi2,1); ny = ny1 + ny2; 

 ga1 = fi1*lp';   ga2 = fi2*la';

 Ga1 = double(int(ga1,x,-r1,0)*diag(2.*(0:d)+1)*(1/r1));

 Ga2 = double(int(ga2,x,-r2,-r1)*diag(2.*(0:d)+1)*(1/r3));

 ep1 = fi1 - Ga1*lp;   ep2 = fi2 - Ga2*la;

 E1  = double(int(ep1*ep1',x,-r1,0)); E2  = double(int(ep2*ep2',x,-r2,-r1));

The code works fine until d = 8 when an error is returned to state that DOUBLE cannot convert the input expression into a double array. If the input expression contains a symbolic variable, use VPA.

I have tried vpa function but the same problem happens still.

One may suggest to use integral numerical integration instead of int. However, the numerical integration seems produce inaccurate result compared to the symbolic representation. Note that the error matrix E1 and E2 become extremely small when the approximation degree d becomes large.

To summarize, the problem here is how to extract, or accurately calculate if anyone has suggestions, the numerical values of E1 and E2.

Thanks a lot!

[John D'Errico]

You say that we can try it ourselves, but since there are undefined variables, we cannot do so.

Undefined function or variable 'n'.

So trying it ourselves stops at line 1.

[Qian Feng]

Sorry, now the code should be correct

Answer 1

[Walter Roberson]

If you change your first line to

 syms n
 r1 = 0.7;  r2 = 1; r3 = r2 - r1; d = 8; di = n*(d+1);

then you can complete down through E1. After that, unfortunately MATLAB does not know how to do the integral, even though there is a closed-form solution for it. You will need to switch to numeric:

FF = ep2*ep2';
FF1 = matlabFunction(simplify(FF(1));
FF2 = matlabFunciton(simplify(FF(2));
E2(1) = integral(FF1, -r2, -r1);
E2(2) = integral(FF2, -r2, -r1);

[Qian Feng]

but when d<8, why the value of E2 can be obtained in fact? This is one of the reason I posed this question since it seems weird that the code only fails to produce result when d become relatively higher. Try d = 1...7 then you can see what I am talking about.

Sorry, di and n are redundant in my code and they should not be there in the first place.

[Walter Roberson]

I believe the appropriate answer is "I don't know." It is not obvious to me why there might be problems with the numeric integral.

The other half of the answer is: considering you are converting to double() immediately after, is there a reason not to duck the issue by doing numeric integration over the function handle?

[Qian Feng]

Thank you very much for the answer. Now first I have a question. What is the mechanism which Matlab applies when I use the syntax double(Ga1) ? Does it automatically perform a numerical integration or use other scheme to get the numerical values of Ga1.

[Walter Roberson]

If you have an symbolic int() that it cannot find a closed form solutions for then when you double() it sends the problem to the mupad routine numeric::int for numeric integration. I do not recall at the moment what the default integration method is.

[Qian Feng]

The reason asked this question is because as I had said in the main question that It occurred that using symbolic produces more accurate results than applying numerical integrations.

[Walter Roberson]

https://www.mathworks.com/help/symbolic/mupad_ref/numeric-int.html

In cases where a closed-form solution can be found, double(int()) involves finding the closed form solution and evaluating it at the given locations. That might involve calculations with large or small magnitudes or have taken into account odd shapes that are difficult to calculate accurately numerically in a finite amount of time.

In the case where no closed form solution can be found, the double() triggers using numeric::int which uses quadrature . This still has advantages compared to numeric integration at the MATLAB level because of the extended range and extended digits that symbolics permits.

[Qian Feng]

OK. Based on the information you provided, then there must be some problems occurred to the scenario of numeric::int when d>=8. Do you have some idea about that ? Could the error happened simply because the numerical values become too small to be handled ?

[Karan Gill]

Use "vpaintegral" introduced in 16b: https://www.mathworks.com/help/symbolic/vpaintegral.html.

F2  = vpaintegral(ep2*ep2',x,-r2,-r1)
F2 =
[ 1.53919e-23,  2.0475e-23]
[  2.0475e-23, 2.73446e-23]

[Walter Roberson]

Workaround without using the new vpaintegral:

RRR = ep2*ep2';
E2 = zeros(size(RRR));
for K = 1 : numel(RRR)
  E2(K) = double( sum(int(children(expand(RRR(K))),x,-r2,-r1)) );
end

[Karan Gill]

Walter, that's very fancy ;)

[Qian Feng]

Thanks Walter and Karan, could you guys put your answer into a new thread so I can accept it ? There are many comments in this one.

[Walter Roberson]

You can Accept this one if I answered your question. If you ended up using vpaintegral() then Yes, Karan should create an Answer so that can be accepted.

[Walter Roberson]

(I just created a bug report about the integral failing.)

[Karan Gill]

FWIW I plonked down a separate answer.

[Qian Feng]

OK, I want to accept both answers of you guys but it seems there is only one answer I can accept.

[Walter Roberson]

Correct, only one Answer can be accepted at present. But you can Vote for the other.

[Qian Feng]

OK, So the integral failing can be considered as a bug ? Thanks a lot for highlighting it.

[Walter Roberson]

Mathworks tells me that double() fails for this because it uses fewer digits and apparently there might be some convergence problem in that case. double(vpa()) is another way of approaching it.

[Qian Feng]

double(vpa()) still generated error when d = 8

[Walter Roberson]

Use a higher digits for vpa() and allow double() to reduce it.

[Qian Feng]

Could you explain the meaning 'allow double() to reduce it' ? Is there an option for double() to do this ?

[Walter Roberson]

SomeLargeNumberOfDigits = 32;

double( vpa(YourExpression, SomeLargeNumberOfDigits) )

[Karan Gill]

About "Allow double() to reduce it", please see: https://www.mathworks.com/help/symbolic/choose-symbolic-or-numeric-arithmetic.html

[Qian Feng]

Ok, is this about reducing the Number of significant digits of vpa function ?

[Walter Roberson]

No, you need the number of digits of vpa to remain high so that vpa is able to converge. But then you can transform those higher number of symbolic digits into a numeric floating point value with double().

... or just use vpaintegral() directly if you have a new enough MATLAB.

[Qian Feng]

Right, so I need larger valued of digits for vpa function.

Answer 2

[Karan Gill]

Use "vpaintegral" introduced in 16b: https://www.mathworks.com/help/symbolic/vpaintegral.html.

F2  = vpaintegral(ep2*ep2',x,-r2,-r1)
F2 =
[ 1.53919e-23,  2.0475e-23]
[  2.0475e-23, 2.73446e-23]

[Qian Feng]

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;

Karan, see the attached code here. There are difficulties to get the value of Ga2. If we put the breakpoint at Ga2, you can see that Ga1 and E1 can be calculated within a reasonable period. The functions inside of the integrations contain no singularities I believe, so I am puzzled here since the type of functions inside of the integrations among Ga1, Ga2, are the same.

[Walter Roberson]

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]

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

[Walter Roberson]

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.