码字不易，觉得好，点个赞或收藏下吧！😊
数值分析编程作业3
预备知识
对于n次拉格朗日插值而言，对应点
(
x
i
,
y
i
)
(x_i,y_i)
(
x
i
,
y
i
)
的插值基函数
s.t.
l
i
(
x
k
)
=
{
y
i
k
=
i
0
k
≠
i
l_{i}(x_k)= \begin{cases} y_i& k=i\\ 0& k\not=i \end{cases}
l
i
(
x
k
)
=
{
y
i
0
k
=
i
k
=
i
也即
l
i
(
x
)
=
∏
j
=
0
n
(
x
−
x
j
)
∏
j
=
0
n
(
x
i
−
x
j
)
(
i
=
0
,
1
,
2
,
.
.
.
,
n
i
≠
j
)
l_i(x) = \frac{\prod_{j=0}^n(x-x_j)}{\prod_{j=0}^n(x_i-x_j)}\quad( i = 0,1,2,…,n \quad i\not=j)
l
i
(
x
)
=
∏
j
=
0
n
(
x
i
−
x
j
)
∏
j
=
0
n
(
x
−
x
j
)
(
i
=
0
,
1
,
2
,
.
.
.
,
n
i
=
j
)
可得到插值多项式为
L
n
(
x
)
=
∑
i
=
0
n
l
i
(
x
)
y
i
(
i
=
0
,
1
,
2
,
.
.
.
,
n
)
L_n(x) = \sum_{i=0}^{n}l_i(x)y_i\quad( i = 0,1,2,…,n )
L
n
(
x
)
=
i
=
0
∑
n
l
i
(
x
)
y
i
(
i
=
0
,
1
,
2
,
.
.
.
,
n
)
此外，对于不同的插值点的选取，也会影响到插值效果的好坏，也就是插值余项
R
n
(
x
)
=
f
(
n
+
1
)
(
ξ
)
(
n
+
1
)
!
ω
n
+
1
(
x
)
R_n(x)= \frac{f^{(n+1)}(\xi)}{(n+1)!}\omega_{n+1}(x)
R
n
(
x
)
=
(
n
+
1
)
!
f
(
n
+
1
)
(
ξ
)
ω
n
+
1
(
x
)
的大小。当插值次数越高时，会出现龙格现象，导致插值函数的插值结果偏离原函数很多。
程序思路
编写lagrange、chebyshev、omega函数脚本和test测试脚本。
-
lagrange脚本
用来进行拉格朗日插值函数的计算，有三个函数参数，x0,y0,x。其中x0和y0参数是插值节点对应的插值节点的x,y数值向量，x是插值区间上一定等距步长的坐标点向量，用于绘制多项式插值函数。
function y=lagrange(x0,y0,x)
ii=1:length(x0);
y=zeros(size(x));
for i=ii
ij=find(ii~=i);
y1=1;
for j=1:length(ij)
y1=y1.*(x-x0(ij(j)));
end
y=y+y1*y0(i)/prod(x0(i)-x0(ij));
end
-
chebyshev脚本
用于求取切比雪夫节点，参数n表示插值多项式的次数，对应地，所需插值节点数为（n+1）。对于
f(
x
)
ϵ
[
−
1
,
1
]
,
有
x
k
=
c
o
s
(
2
k
+
1
)
π
2
(
n
+
1
)
(
k
=
0
,
1
,
.
.
.
,
n
)
f(x)\epsilon[-1,1],有x_k=cos \frac{(2k+1)\pi}{2(n+1)}\quad (k=0,1,…,n)
f
(
x
)
ϵ
[
−
1
,
1
]
,
有
x
k
=
c
o
s
2
(
n
+
1
)
(
2
k
+
1
)
π
(
k
=
0
,
1
,
.
.
.
,
n
)
在实现过程中需要注意matlab下标从1开始，而k是从0开始，所以须将分母中的2k+1更改为2k-1。
function x=chebyshev(n)
x = zeros(1,n+1);
for k = 1:n+1
x(k) = cos(((2*k-1)*pi)/(2*(n+1)));
end
-
omega脚本
用于求取插值多项式的插值余项，有两个参数x和x0,其中x为插值区间上一定等距步长的坐标点向量，err向量的维度与其一致，x0为插值节点的x坐标值所对应的向量。由拉格朗日插值余项
Rn
(
x
)
=
f
(
x
)
−
L
n
(
x
)
=
f
(
n
+
1
)
(
ξ
)
(
n
+
1
)
!
ω
n
+
1
(
x
)
R_n(x)=f(x) – L_n(x) = \frac{f^{(n+1)}(\xi)}{(n+1)!}\omega_{n+1}(x)
R
n
(
x
)
=
f
(
x
)
−
L
n
(
x
)
=
(
n
+
1
)
!
f
(
n
+
1
)
(
ξ
)
ω
n
+
1
(
x
)
由于待插值函数的任意阶导数的最大值为e，则可化简为
Rn
(
x
)
=
f
(
x
)
−
L
n
(
x
)
≤
e
(
n
+
1
)
!
ω
n
+
1
(
x
)
R_n(x)=f(x) – L_n(x) \leq \frac{e}{(n+1)!}\omega_{n+1}(x)
R
n
(
x
)
=
f
(
x
)
−
L
n
(
x
)
≤
(
n
+
1
)
!
e
ω
n
+
1
(
x
)
按照以上公式实现该函数脚本。
function err = omega(x,x0)
err = zeros(1,length(x));
for i = 1:length(x)
err(i) = prod(x(i)*ones(1,length(x0))-x0);
end
err = err*exp(1)/factorial(length(x0));
- test测试脚本通过调用以上函数脚本，实现对插值函数以及插值余项的计算和绘制。
x = -1:0.01:1;
y = exp(1).^abs(x);
%均匀间距插值，使用等间距节点作为插值节点
subplot(2,2,1);
title('均匀间距插值');
x11 = -1:0.2:1; %n=10
y11 = exp(1).^abs(x11);
f11 = lagrange(x11,y11,x);
err11 = exp(1)*omega(x,x11);
plot(x,f11,'-r')
hold on
x12 = -1:0.1:1; %n=20
y12 = exp(1).^abs(x12);
f12 = lagrange(x12,y12,x);
err12 = exp(1)*omega(x,x12);
plot(x,f12,'-b')
hold on
plot(x,y,'-.k')
legend('n=10','n=20','原函数')
%切比雪夫插值，使用切比雪夫节点作为插值节点
subplot(2,2,2);
title('切比雪夫插值');
x21 = chebyshev(10);%n=10
y21 = exp(1).^abs(x21);
f21 = lagrange(x21,y21,x);
err21 = exp(1)*omega(x,x21);
plot(x,f21,'-r');
hold on
x22 = chebyshev(20);%n=20
y22 = exp(1).^abs(x22);
f22 = lagrange(x22,y22,x);
err22 = exp(1)*omega(x,x22);
plot(x,f22,'-b');
plot(x,y,'-.k')
legend('n=10','n=20','原函数')
%绘制插值余项的最大值
subplot(2,2,3);
plot(x,err11,'-r',x,err21,'-b');
title('当为10次插值多项式时的插值余项最值');
legend('等距节点','切比雪夫节点');
subplot(2,2,4);
plot(x,err12,'-r',x,err22,'-b');
title('当为20次插值多项式时的插值余项最值');
legend('等距节点','切比雪夫节点');
实验结果
更正的仿真图：
由结果图可得到，等距节点的拉格朗日插值和切比雪夫节点的拉格朗日插值在n=20时和n=10时区间边缘均有振荡现象，且前者振荡幅度更大，即出现了龙格现象。
同时，仔细观察可发现，切比雪夫节点插值振荡的幅度较等距节点插值振荡幅度较小，再进一步从插值余项的角度进行分析可以看出不论插值多项式的次数是10还是20，切比雪夫节点插值余项的最大值的取值范围均小于等距节点插值余项，这也进一步佐证了使用切比雪夫节点插值可以改善等距节点插值龙格现象的结论。