Numerical values of integrals

9 次查看（过去 30 天）
Qian Feng
Qian Feng 2016-11-7
编辑： Qian Feng 2021-9-1
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!
  2 个评论
John D'Errico
John D'Errico 2016-11-7
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
Qian Feng 2016-11-17
Sorry, now the code should be correct

请先登录，再进行评论。

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

采纳的回答

Walter Roberson
Walter Roberson 2016-11-7
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);
  26 个评论
显示 24更早的评论
Walter Roberson
Walter Roberson 2016-12-30
编辑：Walter Roberson 2017-9-24
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
Qian Feng 2017-1-4
Right, so I need larger valued of digits for vpa function.

请先登录，再进行评论。

更多回答（1 个）

Karan Gill
Karan Gill 2016-11-30
编辑：Karan Gill 2017-10-17
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]
  4 个评论
显示 2更早的评论
Qian Feng
Qian Feng 2017-9-24
Thanks Walter, so what number do you suggest that I should choose for the AbsTol and RelTol in this case?
Walter Roberson
Walter Roberson 2017-9-24
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.

请先登录，再进行评论。

类别

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

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by