使用切比雪夫点进行插值

Question

Interpolate the Runge function of Example 10.6 at Chebyshev points for n from 10 to 170 in increments of 10. Calculate the maximum interpolation error on the uniform evaluation mesh x = -1:.001:1 and plot the error vs. polynomial degree as in Figure 10.8 using semilogy. Observe spectral accuracy.

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

Answer 1

由于您对 MATLAB 知之甚少，我将尝试逐步解释所有内容：

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

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, a very useful feature of MATLAB. Notice the function is properly 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 function does a polynomial curve fitting - it obtains the coefficients of the interpolating polynomial, given the poins x,y and the degree of the polynomial n. You can easily evaluate the polynomial at other points with the polyval 函数。

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

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

如果不使用线性 space 向量，而是使用其他点进行插值，则可以减少此错误。一个好的选择是使用 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')

使用切比雪夫点进行插值

Interpolation using chebyshev points

matlab

numerical-methods