使用切比雪夫点进行插值

Question

在n从10到170的切比雪夫点处插入示例10.6的龙格函数 以 10 为增量。计算统一评估的最大插值误差 网格 x = -1:.001:1 并绘制误差与多项式次数的关系，如图 10.8 所示，使用 符号学。观察光谱精度。

runge 函数由下式给出：f(x) = 1 / (1 + 25x^2)

到目前为止我的代码：

x = -1:0.001:1;
n = 170;
i = 10:10:170;
cx = cos(((2*i + 1)/(2*(n+1)))*pi); %chebyshev pts

y = 1 ./ (1 + 25*x.^2); %true fct
%chebyshev polynomial, don't know how to construct using matlab

yc = polyval(c, x); %graph of approx polynomial fct
plot(x, yc);
mErr =  (1 / ((2.^n).*(n+1)!))*%n+1 derivative of f evaluated at max x in [-1,1], not sure how to do this
%plotting stuff

我对matlab知之甚少，所以我在努力创建插值多项式。我做了一些谷歌工作，但我对当前的函数感到困惑，因为我没有找到一个只是简单地接受点和要插值的多项式的函数。在这种情况下，我也有点困惑我是否应该做i = 0:1:n 和n=10:10:170 或者n 是否固定在这里。感谢您的帮助，谢谢

Answer 1

由于你对 MATLAB 知之甚少，我将尝试一步一步地解释一切：

首先，要可视化 Runge 函数，您可以键入：

f = @(x) 1./(1+25*x.^2); % Runge function

% plot Runge function over [-1,1];
x = -1:1e-3:1;
y = f(x);
figure;
plot(x,y); title('Runge function)'); xlabel('x');ylabel('y');

代码的@(x) 部分是function handle，这是MATLAB 的一个非常有用的功能。请注意，该函数是正确的vecotrized，因此它可以接收变量或数组作为参数。绘图功能很简单。

要了解龙格现象，请考虑由 10 个元素组成的 [-1,1] 的 linearly spaced 向量，并使用这些点来获得插值（拉格朗日）多项式。您会得到以下信息：

% 10 linearly spaced points
xc = linspace(-1,1,10);
yc = f(xc);
p = polyfit(xc,yc,9); % gives the coefficients of the polynomial of degree 10
hold on; plot(xc,yc,'o',x,polyval(p,x)); 

polyfit 函数进行多项式曲线拟合 - 它获得插值多项式的系数，给定点 x,y 和多项式的次数 n。您可以使用polyval 函数轻松计算其他点的多项式。

请注意，靠近末端域，您会得到一个振荡多项式，并且插值不是该函数的良好近似值。事实上，您可以绘制绝对误差，比较函数f(x) 和插值多项式p(x) 的值：

plot(x,abs(y-polyval(p,x))); xlabel('x');ylabel('|f(x)-p(x)|');title('Error');

如果不使用线性空间向量，而是使用其他点进行插值，则可以减少此错误。一个不错的选择是使用Chebyshev nodes，这样可以减少错误。事实上，请注意：

% find 10 Chebyshev nodes and mark them on the plot
n = 10;
k = 1:10; % iterator
xc = cos((2*k-1)/2/n*pi); % Chebyshev nodes
yc = f(xc); % function evaluated at Chebyshev nodes
hold on;
plot(xc,yc,'o')

% find polynomial to interpolate data using the Chebyshev nodes
p = polyfit(xc,yc,n-1); % gives the coefficients of the polynomial of degree 10
plot(x,polyval(p,x),'--'); % plot polynomial
legend('Runge function','Chebyshev nodes','interpolating polynomial','location','best')

请注意错误是如何减少接近末端域的。你现在没有得到插值多项式的高振荡行为。如果您绘制错误，您将观察到：

plot(x,abs(y-polyval(p,x))); xlabel('x');ylabel('|f(x)-p(x)|');title('Error');

现在，如果您更改切比雪夫节点的数量，您将获得更好的近似值。对代码稍作修改，您就可以针对不同数量的节点再次运行它。您可以存储最大误差并将其绘制为节点数的函数：

n=1:20; % number of nodes

% pre-allocation for speed
e_ln = zeros(1,length(n)); % error for the linearly spaced interpolation
e_cn = zeros(1,length(n)); % error for the chebyshev nodes interpolation

for ii=1:length(n)
    % linearly spaced vector
    x_ln = linspace(-1,1,n(ii)); y_ln = f(x_ln);
    p_ln = polyfit(x_ln,y_ln,n(ii)-1);
    e_ln(ii) = max( abs( y-polyval(p_ln,x) ) );

    % Chebyshev nodes
    k = 1:n(ii); x_cn = cos((2*k-1)/2/n(ii)*pi); y_cn = f(x_cn);
    p_cn = polyfit(x_cn,y_cn,n(ii)-1);
    e_cn(ii) = max( abs( y-polyval(p_cn,x) ) );
end
figure
plot(n,e_ln,n,e_cn);
xlabel('no of points'); ylabel('maximum absolute error');
legend('linearly space','chebyshev nodes','location','best')