多项式近似数据时出错，请求排查MATLAB数值分析代码问题

题目翻译 #

问题2：
考虑区间 [0,1] 上的函数 ( f(x) = \arctan(x(x+1)) )
a) 使用 [0,1] 上的8个等距节点构造7次插值多项式 ( P_7(x) )，计算 ( P_7(0.7) ) 的值以及相对误差 ( \frac{|f(0.7)-P_7(0.7)|}{|f(0.7)|} )。
b) 使用相同节点构造自然三次样条插值函数，计算该样条在 ( x=0.7 ) 处的值及相对误差。
c) 重复b部分，构造固支三次样条（端点导数等于 ( f'(0) ) 和 ( f'(1) )）。

代码问题与修正 #

1. 相对误差计算不符合题目要求 #

你的代码中相对误差的分母未取绝对值，虽然当前 ( f(0.7) ) 为正值，但严格不符合题目规范，需修改为：

err = abs(f(0.7)-p)/abs(f(0.7))

2. 冗余代码 #

y_val = f(xplot); 这一行未被后续代码调用，属于冗余代码，可直接删除。

3. 高次插值的数值稳定性问题 #

7次多项式插值属于高次插值，polyfit在处理这类插值时，可能因数值条件差导致系数出现舍入误差，进而影响插值结果的精度。可以改用更稳定的拉格朗日插值实现：

拉格朗日插值函数实现

function y = lagrange_interp(x_nodes, y_nodes, x)
    n = length(x_nodes);
    y = zeros(size(x));
    for i = 1:n
        l = ones(size(x));
        for j = 1:n
            if j ~= i
                l = l .* (x - x_nodes(j)) / (x_nodes(i) - x_nodes(j));
            end
        end
        y = y + y_nodes(i) * l;
    end
end

修正后的7次插值代码

close all; clear all; clc;
format long e

x = linspace(0,1,8);
f = @(x) atan(x.*(x+1));
y_nodes = f(x);
% 使用拉格朗日插值计算P7(0.7)
p = lagrange_interp(x, y_nodes, 0.7);
err = abs(f(0.7)-p)/abs(f(0.7))

若需完成样条插值部分（题目b/c） #

如果你的目标是完成题目中的样条插值要求，可使用MATLAB的csape函数（需样条工具箱）：

自然三次样条（题目b）

close all; clear all; clc;
format long e

x = linspace(0,1,8);
f = @(x) atan(x.*(x+1));
y_nodes = f(x);
% 构造自然三次样条
cs_natural = csape(x, y_nodes, 'variational');
p_spline_natural = fnval(cs_natural, 0.7);
err_natural = abs(f(0.7)-p_spline_natural)/abs(f(0.7))

固支三次样条（题目c）

close all; clear all; clc;
format long e

x = linspace(0,1,8);
f = @(x) atan(x.*(x+1));
% 计算f(x)的导数
f_prime = @(x) (2*x + 1)./(1 + (x.*(x+1)).^2);
y_nodes = f(x);
% 构造固支三次样条，指定端点导数
cs_clamped = csape(x, y_nodes, 'clamped', [f_prime(0), f_prime(1)]);
p_spline_clamped = fnval(cs_clamped, 0.7);
err_clamped = abs(f(0.7)-p_spline_clamped)/abs(f(0.7))

内容的提问来源于stack exchange，提问作者curiousityspark