正交多项式的数值评估
Numerical evaluation of orthogonal polynomials
我已经编写了一些评估正交多项式的 Matlab 程序,作为完整性检查,我试图确保它们的点积为零。
但是,虽然我相当确定不会出错,但我发现自己的行为有点奇怪。我的测试很简单:
x = -1:.01:1;
for i0=0:9
v1 = poly(x, i0);
for i1=0:i0
v2 = poly(x,i1);
fprintf('%d, %d: %g\n', i0, i1, v1*v2');
end
end
(注意点积 v1*v2'
需要这样,因为 x
是水平向量。)
现在,要切入故事的结尾,对于加起来为奇数的成对度数(即 i0+i1=2k+1
). i0==i1
的时候我希望点积不为0,但是i0+i1=2k
的时候也是这样,没想到
为了给你更多的细节,我最初用第一类切比雪夫多项式做了这个测试。现在,它们与权重正交
1 ./ sqrt(1-x.^2)
当 x
变为 1 时,它变为无穷大。所以我认为将这一项排除在外可能是非零点积的原因。
但是后来,我用勒让德多项式做了同样的测试,得到了完全相同的结果:当度数和为偶数时,点积肯定远离0(数量级1e2)。
最后一个细节,我使用了三角函数公式cos(n*acos(x))
来计算切比雪夫多项式,我尝试了递归公式以及其中一个涉及二项式系数的公式来计算勒让德多项式。
谁能解释这种奇怪的(双关语)行为?
你被对称误导了。两个 Chebyshev and Legendre 多项式都是奇偶运算符的特征函数,这意味着它们都可以归类为 either 奇函数或偶函数。我想您的自定义正交多项式也是如此。
由于这种对称性,如果将多项式 P_n(x)
乘以 P_m(x)
,则如果 n+m
为奇数,则结果将为奇函数,否则为偶数.您正在计算 sum_k P_n(x_k)*P_m(x_k)
围绕原点的一组对称 x_k
值。这意味着对于奇数 n+m
你将 总是 得到零。尝试用 P
勒让德和 Q
切比雪夫多项式计算 sum_k P_n(x_k)*Q_m(x_k)
。我的观点是,对于 n+m=odd
,结果不会告诉您任何有关正交性或积分准确性的信息。
问题是您可能没有足够准确地集成。这些在 [-1,1]
上定义的正交多项式在其定义域上变化非常快,尤其是在边界附近 (x==+-1
)。尝试增加积分点,使用非等距网格,或使用 integral
.
进行适当的积分
最后说明:我建议您不要调用您的函数 poly
,因为 that's a MATLAB built-in. (And so is legendre
。)
我已经编写了一些评估正交多项式的 Matlab 程序,作为完整性检查,我试图确保它们的点积为零。
但是,虽然我相当确定不会出错,但我发现自己的行为有点奇怪。我的测试很简单:
x = -1:.01:1;
for i0=0:9
v1 = poly(x, i0);
for i1=0:i0
v2 = poly(x,i1);
fprintf('%d, %d: %g\n', i0, i1, v1*v2');
end
end
(注意点积 v1*v2'
需要这样,因为 x
是水平向量。)
现在,要切入故事的结尾,对于加起来为奇数的成对度数(即 i0+i1=2k+1
). i0==i1
的时候我希望点积不为0,但是i0+i1=2k
的时候也是这样,没想到
为了给你更多的细节,我最初用第一类切比雪夫多项式做了这个测试。现在,它们与权重正交
1 ./ sqrt(1-x.^2)
当 x
变为 1 时,它变为无穷大。所以我认为将这一项排除在外可能是非零点积的原因。
但是后来,我用勒让德多项式做了同样的测试,得到了完全相同的结果:当度数和为偶数时,点积肯定远离0(数量级1e2)。
最后一个细节,我使用了三角函数公式cos(n*acos(x))
来计算切比雪夫多项式,我尝试了递归公式以及其中一个涉及二项式系数的公式来计算勒让德多项式。
谁能解释这种奇怪的(双关语)行为?
你被对称误导了。两个 Chebyshev and Legendre 多项式都是奇偶运算符的特征函数,这意味着它们都可以归类为 either 奇函数或偶函数。我想您的自定义正交多项式也是如此。
由于这种对称性,如果将多项式 P_n(x)
乘以 P_m(x)
,则如果 n+m
为奇数,则结果将为奇函数,否则为偶数.您正在计算 sum_k P_n(x_k)*P_m(x_k)
围绕原点的一组对称 x_k
值。这意味着对于奇数 n+m
你将 总是 得到零。尝试用 P
勒让德和 Q
切比雪夫多项式计算 sum_k P_n(x_k)*Q_m(x_k)
。我的观点是,对于 n+m=odd
,结果不会告诉您任何有关正交性或积分准确性的信息。
问题是您可能没有足够准确地集成。这些在 [-1,1]
上定义的正交多项式在其定义域上变化非常快,尤其是在边界附近 (x==+-1
)。尝试增加积分点,使用非等距网格,或使用 integral
.
最后说明:我建议您不要调用您的函数 poly
,因为 that's a MATLAB built-in. (And so is legendre
。)