数值分析实验 - 函数插值方法

发表于 更新于 分类于 ， ， ， 本文字数： 1.4k 阅读时长 ≈ 5 分钟

数值分析实验 2 - 函数插值方法

实验要求

  • 问题 1

    对如下结点构造五次插值多项式和分段三次插值多项式

    \(x_i\) 0.4 0.55 0.65 0.80 0.95 1.05
    \(y_i\) 0.41075 0.57815 0.69675 0.90 1.00 1.25382

    并计算 \(f(0.596),f(0.99)\)

  • 问题 2

    对如下结点构造六次插值多项式和分段三次插值多项式

    \(x_i\) 1 2 3 4 5 6 7
    \(y_i\) 0.368 0.135 0.050 0.018 0.007 0.002 0.001

    并计算 \(f(1.8),f(6.15)\)

计算公式

本博客分别使用 MATLAB 实现了 Lagrange 插值，Neville 插值，Newton 插值和三次样条插值

  • Lagrange 插值 \[ L(x)=\sum_{k=0}^ny_k\prod_{i=0;~i\ne k}^n\frac{x-x_i}{x_k-x_i}\tag{1} \]

  • Neville 插值

    令 \(N\in\mathbb{R}^{(n+1)\times (n+1)}[x]\) 满足

    • \[ N_{i,1}(x)=y_{i-1} \]
    • \[ N_{i,j}(x)={(x-x_{i-1})N_{i-1,j-1}(x)+(x_{i-j}-x)N(i,j-1)\over x_{i-j}-x_{i-1}} \]

    则结果即为

    \[ N_{n+1,n+1}(x)\tag{2} \]

  • Newton 插值

    \[ N(x)=f(x_0)+\sum_{k=1}^nf[x_0,x_1,...,x_k]\prod_{i=0}^{k-1}(x-x_i)\tag{3} \]

    其中

    • \[ f[a,b]=\frac{f(a)-f(b)}{a-b}\tag{4} \]
    • \[ f[x_0,x_1,...,x_k]={f[x_0,x_1,...,x_{k-1}]-f[x_1,x_2,...,x_k]\over x_0-x_k}\tag{5} \]

  • 三次样条插值 (本文采用自然样条插值)

    设结果 \(S(x)\in C^2[x_0,x_n]\) 满足

    • \(S(x_i)=y_i\)
    • \(\partial S(x)\leqslant3,x\in[x_{i-1},x_i]\)

    令

    • \(h_i=x_i-x_{i-1}\)
    • \(\lambda_i={h_i\over h_i+h_{i+1}}\)
    • \(\mu_i=1-\lambda_i\)
    • \(g_i=6f[x_{i-1},x_i,x_{i+1}]\)
    • \(A_i=6y_{i-1}+M_{i-1}h_i^2\)
    • \(B_i=6y_i+M_ih_i^2\)
    • \(\displaystyle S_i(x)={M_{i-1}(x_i-x)^3+M_i(x-x_{i-1})^3\over 6h_i}+{A_i(x_i-x)+B_i(x-x_{i-1})\over 6h_i},~x\in[x_{i-1},x_i],i=1,2,...,n\)
    • \(\displaystyle S(x)=\sum_{i=1}^nS_i(x)\)

    由 \(\lim_{x\to k^+}S'(x)=\lim_{x\to k^-}S'(x),~\forall k\in[x_0,x_n]\) 知

    \[ \lambda_iM_{i-1}+2M_i+\mu_iM_{i+1}=g_i,~i=1,2,...,n-1\tag{6} \]

    又

    \[ S''(x_0)=S''(x_n)=0\tag{7} \]

    联立 \((6),(7)\), 有

    • \(M_0=M_n=0\)
    • \[ \begin{bmatrix} 2&\mu_1&&&&\\ \lambda_2&2&\mu_2&&&\\ &\lambda_3&2&\mu_3&&\\ &&\ddots&\ddots&\ddots&\\ &&&\lambda_{n-2}&2&\mu_{n-2}\\ &&&&\lambda_{n-1}&2\\ \end{bmatrix}\begin{bmatrix} M_1\\M_2\\\vdots\\M_{n-1} \end{bmatrix}=\begin{bmatrix} g_1\\g_2\\\vdots\\g_{n-1} \end{bmatrix}\tag{8} \]

解出 \(M_i,~i=0,1,...,n\), 即求得 \(S(x)\)

程序设计

主程序

Show code

main.mview raw
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
% Exp.2

% @Author: Tifa
% @LastEditTime: 2021-04-16 20:49:47

% Question flag
now_question = 1;

% Data
if now_question == 1
    x = [0.4, 0.55, 0.65, 0.80, 0.95, 1.05];
    y = [0.41075, 0.57815, 0.69675, 0.90, 1.00, 1.25382];
    xx = [0.596, 0.99];
    plot_x = linspace(0.2, 1.2, 1000);
else
    x = [1, 2, 3, 4, 5, 6, 7];
    y = [0.368, 0.135, 0.050, 0.018, 0.007, 0.002, 0.001];
    xx = [1.8, 6.15];
    plot_x = linspace(0, 9, 1000);
end

% Lagrange interpolating polynomial
p1 = calc_lagrange(x, y);
Y1 = polyval(p1, xx);
disp(vpa(poly2sym(p1), 8))

subplot(2, 2, 1)
plot_y = polyval(p1, plot_x);
hold on
plot(x, y, 'k*');
plot(plot_x, plot_y, 'k');
plot(xx, Y1, 'ko');
title('{\bf Fig. 1} Lagrange interpolating polynomial')
legend('Input points', 'Interpolating polynomial', 'Query points');
hold off

% Cubic spline interpolating polynomial
pp2 = calc_spline3(x, y);
Y2 = ppval(pp2, xx);

for i = 1:size(pp2.coefs, 1)
    disp(vpa(poly2sym(pp2.coefs(i, :)), 8))
end

subplot(2, 2, 2)
plot_y = ppval(pp2, plot_x);
hold on
plot(x, y, 'k*');
plot(plot_x, plot_y, 'k');
plot(xx, Y2, 'ko');
title('{\bf Fig. 2} Cubic spline interpolating polynomial')
legend('Input points', 'Interpolating polynomial', 'Query points');
axis
hold off

% Neville interpolating polynomial
p3 = calc_neville(x, y);
Y3 = polyval(p3, xx);
disp(vpa(poly2sym(p3), 8))

subplot(2, 2, 3)
plot_y = polyval(p3, plot_x);
hold on
plot(x, y, 'k*');
plot(plot_x, plot_y, 'k');
plot(xx, Y3, 'ko');
title('{\bf Fig. 3} Neville interpolating polynomial')
legend('Input points', 'Interpolating polynomial', 'Query points');
axis
hold off

% Newton interpolating polynomial
p4 = calc_newton(x, y);
Y4 = polyval(p4, xx);
disp(vpa(poly2sym(p4), 8))

subplot(2, 2, 4)
plot_y = polyval(p4, plot_x);
hold on
plot(x, y, 'k*');
plot(plot_x, plot_y, 'k');
plot(xx, Y4, 'ko');
title('{\bf Fig. 4} Newton interpolating polynomial')
legend('Input points', 'Interpolating polynomial', 'Query points');
axis
hold off

输入数据检查

Show code

input_check.mview raw
1
2
3
4
5
6
7
8
9
10
11
function input_check(input_x, input_y)
    % Input check of Exp.2

    % @Author: Tifa
    % @LastEditTime: 2021-04-16 20:49:47

    if ~isvector(input_x) ||~isvector(input_y) || length(input_x) ~= length(input_y)
        error('Invalid input!')
    end

end

Lagrange 插值

Show code

calc_lagrange.mview raw
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
function poly_out = calc_lagrange(input_x, input_y)
    %calc_lagrange - Calculating Lagrange interpolating polynomial with given points
    %
    % Syntax: poly_out = calc_lagrange(input_x, input_y)

    % @Author: Tifa
    % @LastEditTime: 2021-04-16 20:49:47

    input_check(input_x, input_y);

    len = length(input_x);
    base_funcs = ones(len);

    for i = 1:len
        now_func = 1;

        for j = 1:len

            if i ~= j
                now_func = conv(now_func, poly(input_x(j))) / (input_x(i) - input_x(j));
            end

        end

        base_funcs(i, :) = now_func;
    end

    poly_out = input_y * base_funcs;

end

Neville 插值

Show code

calc_neville.mview raw
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
function poly_out = calc_neville(input_x, input_y)
    %calc_neville - Calculating Neville interpolating polynomial with given points
    %
    % Syntax: poly_out = calc_neville(input_x, input_y)

    % @Author: Tifa
    % @LastEditTime: 2021-04-16 20:49:47

    input_check(input_x, input_y);

    if isrow(input_y)
        input_y = input_y';
    end

    len = length(input_x);
    A = sym(zeros(len)); A(:, 1) = sym(input_y)';

    syms x

    for i = 2:len

        for j = i:len
            A(j, i) =- ((x - input_x(j)) * A(j - 1, i - 1) + (input_x(j - i + 1) - x) * A(j, i - 1)) / (input_x(j) - input_x(j - i + 1));
        end

    end

    poly_out = sym2poly(A(len, len));

end

Newton 插值

Show code

calc_newton.mview raw
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
function poly_out = calc_newton(input_x, input_y)
    %calc_newton - Calculating Newton interpolating polynomial with given points
    %
    % Syntax: poly_out = calc_newton(input_x, input_y)
    input_check(input_x, input_y);

    % @Author: Tifa
    % @LastEditTime: 2021-04-16 20:49:47

    if isrow(input_y)
        input_y = input_y';
    end

    len = length(input_x);
    diff_mat = zeros(len); diff_mat(:, 1) = input_y';

    for i = 2:len

        for j = i:len
            diff_mat(j, i) = (diff_mat(j, i - 1) - diff_mat(j - 1, i - 1)) / (input_x(j) - input_x(j - i + 1));
        end

    end

    poly_out = diff_mat(len, len);

    for i = len - 1:-1:1
        poly_out = conv(poly_out, poly(input_x(i)));
        poly_out(end) = poly_out(end) + diff_mat(i, i);
    end

end

三次样条插值

Show code

calc_spline3.mview raw
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
function pp_out = calc_spline3(input_x, input_y)
    %calc_spline3 - Calculating cubic spline interpolating polynomial with given points
    %
    % Syntax: pp_out = calc_spline3(input_x, input_y)
    %
    % You can also use 'pp_out = spline(input_x, input_y)' instead

    % @Author: Tifa
    % @LastEditTime: 2021-04-16 20:49:47

    input_check(input_x, input_y);

    if ~issorted(input_x)
        error('X coordinates of input points should be in ascending order!')
    end

    if iscolumn(input_x)
        input_x = input_x';
    end

    if iscolumn(input_y)
        input_y = input_y';
    end

    len = length(input_x);

    diff_matrix = zeros(len);
    diff_matrix(:, 1) = input_y';

    for i = 2:len

        for j = i:len
            diff_matrix(j, i) = (diff_matrix(j, i - 1) - diff_matrix(j - 1, i - 1)) / (input_x(j) - input_x(j - i + 1));
        end

    end

    h = input_x(2:len) - input_x(1:len - 1);
    g = 6 * (diff_matrix(3:len, 2) - diff_matrix(2:len - 1, 2)) ./ (h(1:len - 2) + h(2:len - 1))';

    m = [0; (eye(len - 2) * 2 + ...
            diag(h(2:len - 2) ./ (h(1:len - 3) + h(2:len - 2)), 1) + ...
            diag(h(2:len - 2) ./ (h(2:len - 2) + h(3:len - 1)), -1)) \ g; 0];

    polys = zeros(len - 1, 4);

    for i = 1:len - 1
        syms x
        polys(i, :) = sym2poly(m(i, 1) * (h(i) - x)^3 / (6 * h(i)) + ...
            m(i + 1, 1) * x^3 / (6 * h(i)) + ...
            (diff_matrix(i, 1) - m(i, 1) / 6 * (h(i))^2) * (h(i) - x) / h(i) + ...
            (diff_matrix(i + 1, 1) - m(i + 1, 1) / 6 * (h(i))^2) * x / h(i));
    end

    pp_out = mkpp(input_x, polys);

end

结果讨论和分析

结果

  • 问题 1

    方法 多项式 \(f(0.596)\) \(f(0.99)\)
    Lagrange 插值 \(P_1\) \(0.62573238\) \(1.0542298\)
    Neville 插值 \(P_1\) \(0.62573238\) \(1.0542298\)
    Newton 插值 \(P_1\) \(0.62573238\) \(1.0542298\)
    三次样条插值 \(S_1\) \(0.62896167\) \(1.0842113\)

    其中

    • \[ P_1=121.62636\,x^5 - 422.75031\,x^4 + 572.56675\,x^3 - 377.25487\,x^2 + 121.97184\,x - 15.084523 \]
    • \[ S_1=\begin{cases} -0.56178215\,(x-0.4)^3 + 1.1286401\,(x-0.4) + 0.41075,&x\in[0.4, 0.55)\\ 12.056039\,(x-0.55)^3 - 0.25280197\,(x-0.55)^2 + 1.0907198\,(x-0.55) + 0.57815,&x\in[0.55, 0.65)\\ -24.508536\,(x-0.65)^3 + 3.3640098\,(x-0.65)^2 + 1.4018406\,(x-0.65) + 0.69675,&x\in[0.65, 0.80)\\ 47.096624\,(x-0.80)^3 - 7.6648315\,(x-0.80)^2 + 0.75671734\,(x-0.80) + 0.9,&x\in[0.80, 0.95)\\ -45.095498\,(x-0.95)^3 + 13.528649\,(x-0.95)^2 + 1.63629\,(x-0.95) + 1.0,&x\in[0.95, 1.05] \end{cases} \]

  • 问题 2

    方法 多项式 \(f(1.8)\) \(f(6.15)\)
    Lagrange 插值 \(P_2\) \(0.16476189\) \(0.0012658255\)
    Neville 插值 \(P_2\) \(0.16476189\) \(0.0012658255\)
    Newton 插值 \(P_2\) \(0.16476189\) \(0.0012658255\)
    三次样条插值 \(S_2\) \(0.17116591\) \(0.0016228947\)

    其中

    • \[ P_2=0.000058333333\,x^6 - 0.0016083333\,x^5 + 0.018583333\,x^4 - 0.11754167\,x^3 + 0.44185833\,x^2 - 0.96835\,x + 0.995 \]
    • \[ S_2=\begin{cases} 0.036229487\,(x-1)^3 - 0.26922949\,(x-1) + 0.368,&x\in[1,2)\\ -0.033147436\,(x-2)^3 + 0.10868846\,(x-2)^2 - 0.16054103\,(x-2) + 0.135,&x\in[2,3)\\ 0.0013602564\,(x-3)^3 + 0.0092461538\,(x-3)^2 - 0.04260641\,(x-3) + 0.05,&x\in[3,4)\\ -0.0042935897\,(x-4)^3 + 0.013326923\,(x-4)^2 - 0.020033333\,(x-4) + 0.018,&x\in[4,5)\\ 0.00081410256\,(x-5)^3 + 0.00044615385\,(x-5)^2 - 0.0062602564\,(x-5) + 0.007,&x\in[5,6)\\ -0.00096282051\,(x-6)^3 + 0.0028884615\,(x-6)^2 - 0.002925641\,(x-6) + 0.002,&x\in[6,7] \end{cases} \]

分析

  1. Lagrange 插值，Neville 插值和 Newton 插值得到的多项式相同

  2. 观察图像发现，Lagrange 插值，Neville 插值和 Newton 插值得到的多项式稳定性较差，尤其在问题 2 中，\(x>7\) 时的图像明显偏离预期

  3. 博客中所给出的 Lagrange 插值，Neville 插值和 Newton 插值程序的时间复杂度均为 \(O(n^2)\), 三次样条插值的时间复杂度为 \(O(n^3)\) , 后面将给出 \(O(n\log^2 n)\) 的 Lagrange 插值算法

  4. 四种方法的优缺点如下

    方法 优点 缺点
    Lagrange 插值 形式直观简洁，推导容易 增加新结点时，原有结果不能复用；数值稳定性问题
    Neville 插值 增加新结点时，原有结果可以复用 形式不直观；数值稳定性问题
    Newton 插值 增加新结点时，原有结果可以复用；形式直观简洁 数值稳定性问题
    三次样条插值 数值稳定性好 形式复杂，时间复杂度高

附：基于分治的 Lagrange 插值算法

令 \(\omega(x)=\prod_{i=0}^n(x-x_i)\) 我们改写式 \((1)\)

\[ \begin{aligned} L(x)&=\sum_{k=0}^ny_k\frac{\omega(x)}{x-x_k}\lim_{x\to x_k}\frac{x-x_k}{\omega(x)}\\ &=\sum_{k=0}^n\frac{y_k}{\omega'(x_k)}\frac{\omega(x)}{x-x_k} \end{aligned} \]

我们考虑分治求 \(L(x)\)

令

  • \[ L_{l,r}(x)=\sum_{k=l}^r\frac{y_k}{\omega'(x_k)}\prod_{i=l;i\ne k}^r(x-x_i) \]
  • \[ m=\left\lfloor\frac{l+r}{2}\right\rfloor \]

则有

\[ L_{l,r}(x)=L_{l,m}(x)\prod_{i=m+1}^r(x-x_i)+L_{m+1,r}(x)\prod_{i=l}^m(x-x_i) \]

基于 FFT 的多项式乘法的时间复杂度为 \(O(n\log n)\), 故该算法的时间复杂度为 \(O(n\log^2 n)\)

MATLAB 程序实现

Show code

calc_lagrange_fast.mview raw
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
function poly_out = calc_lagrange_fast(input_x, input_y)
    %calc_lagrange_fast - Calculating Lagrange interpolating polynomial with given points with O(n log^2 n)
    %
    % Syntax: poly_out = calc_lagrange_fast(input_x, input_y)

    % @Author: Tifa
    % @LastEditTime: 2021-04-16 20:49:47

    input_check(input_x, input_y);

    d_omega = domg(1, length(input_x));
    poly_out = llr(1, length(input_x));

    function output = domg(l, r)
    %domg - Calculating {\omega'(x_k)} based on divide and conquer

        if l == r
            output = input_x - input_x(l);
            output(l) = 1;
            return
        end

        mid = floor(l + (r - l) / 2);
        output = domg(l, mid) .* domg(mid + 1, r);
    end

    function [poly_ret, omega_l, omega_r] = llr(l, r)
    %llr - Calculating {L_{l,r}} based on divide and conquer

        if l == r
            poly_ret = input_y(l) / d_omega(l);
            omega_l = 1;
            omega_r = 1;
            return
        end

        mid = floor(l + (r - l) / 2);
        [poly_ret1, omega_l1, omega_r1] = llr(l, mid);
        [poly_ret2, omega_l2, omega_r2] = llr(mid + 1, r);
        omega_l = conv(omega_l1, omega_r1);
        omega_r = conv(omega_l2, omega_r2);

        if mid == l
            omega_l = [1, -input_x(l)];
        end

        if mid == r - 1
            omega_r = [1, -input_x(r)];
        end

        tmp_l = conv(poly_ret1, omega_r); tmp_r = conv(poly_ret2, omega_l);
        len_l = length(tmp_l); len_r = length(tmp_r);

        if len_l > len_r
            tmp = tmp_l; tmp_l = tmp_r; tmp_r = tmp;
            tmp = len_l; len_l = len_r; len_r = tmp;
        end

        tmp_r(len_r - len_l + 1:end) = tmp_r(len_r - len_l + 1:end) + tmp_l;
        poly_ret = tmp_r;

    end

end