数值分析实验 - 函数逼近与曲线拟合

发表于 更新于 分类于 ， ， ， 本文字数： 633 阅读时长 ≈ 2 分钟

数值分析实验 3 - 函数逼近与曲线拟合

实验要求和目的

在某冶炼过程中，根据统计数据的含碳量与时间关系，试求含碳量 \(y\) 与时间 \(t\) 的拟合曲线

\(t~(\text{min})\)0510152025
\(y~(\times 10^{-4})\)01.272.162.863.443.87
\(t~(\text{min})\)303540455055
\(y~(\times 10^{-4})\)4.154.374.514.584.024.64

要求

  1. 用最小二乘法进行曲线拟合
  2. 近似解析表达式为 \(\varphi(t)=a_1t+a_2t^2+a_3t^3\)
  3. 打印出拟合函数 \(\varphi(t)\), 并打印出 \(\varphi(t_j)\) 与 \(y(t_j)\) 的误差，\(j=1,2,...,12\)
  4. 另外选取一个近似表达式，尝试拟合效果的比较
  5. * 绘制出曲线拟合图

计算公式

  • 最小二乘法原理 (参见 笔记 - 对称双线性度量空间与线性方程组)

    令 \(A=(\alpha_0,\alpha_1,\dots,\alpha_n),~\alpha_i\in\mathbb{R}^m,i=1,2,...,n\), 寻找 \(AX=B\) 的最小二乘解 \((x_1,x_2,\dots,x_n)^T\) 即为寻找一组实数 \(x_1,x_2,\dots,x_n\) 使得

    \[ \left|B-\sum_{i=1}^nx_i\alpha_i\right|^2=\sum_{i=1}^n\left(b_1-\sum_{j=1}^na_{ij}x_j\right)^2\tag{1} \]

    的值最小

    此时取的 \(x_1,x_2,\dots,x_n\) 只需使 \(\sum_{i=1}^nx_i\alpha_i\) 为 \(B\) 在 \(G[\alpha_0,\alpha_1,\dots,\alpha_n]\) 的正射影

    这样的 \(x_1,x_2,\dots,x_n\) 是且仅是

    \[ G(\alpha_0,\alpha_1,\dots,\alpha_n)X=((B,\alpha_0),(B,\alpha_1),\dots,(B,\alpha_n))^T\tag{2} \]

    的解，其中 \(G(\alpha_0,\alpha_1,\dots,\alpha_n)\) 为 \(\alpha_0,\alpha_1,\dots,\alpha_n\) 的 Gram 矩阵

    又

    \[ G(\alpha_0,\alpha_1,\dots,\alpha_n)=A^TA \]

    \[ ((B,\alpha_0),(B,\alpha_1),\dots,(B,\alpha_n))^T=A^TB \]

    因此式 \((2)\) 即为

    \[ A^TAX=A^TB\tag{3} \]

    由正射影的存在性可知该方程一定可解，其解为 \(X=(A^TA)^{-1}A^TB\)

程序设计

主程序

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
% Exp.3

% @Author: Tifa
% @LastEditTime: 2021-04-22 23:05:47

% Data
data_x = 0:5:55;
data_y = [0, 1.27, 2.16, 2.86, 3.44, 3.87, 4.15, 4.37, 4.51, 4.58, 4.02, 4.64] * 1e-4;
plot_x = linspace(0, 60, 1000);

% Calculate $\varphi(t)$ and errors
syms x
f1 = lsmfit([x, x^2, x^3], data_x, data_y);
plot_y = subs(f1, plot_x);
err1 = data_y - subs(f1, data_x);
% Print function and errors
disp('function 1:')
vpa(f1, 8)
disp('error:')
fprintf('%.8f ', err1)
fprintf('\n\n')

% Calculate $\varphi_2(t)=b_0+b_1t+b_2t^2+b_3t^3$ and errors
f2 = poly2sym(polyfit(data_x, data_y, 3));
plot_y2 = subs(f2, plot_x);
err2 = data_y - subs(f2, data_x);
% Print function and errors
disp('function 2:')
vpa(f2, 8)
disp('error:')
fprintf('%.8f ', err2)
fprintf('\n\n')

% Plot two functions
subplot(1, 2, 1)
hold on
plot(data_x, data_y, 'k*')
plot(plot_x, plot_y, 'k')
title('Fig. 1')
xlabel('t/minutes')
legend('data', '{f_1}')
hold off

subplot(1, 2, 2)
hold on
plot(data_x, data_y, 'k*')
plot(plot_x, plot_y2, 'k')
title('Fig. 2')
xlabel('t/minutes')
legend('data', '{f_2}')
hold off

最小二乘法函数文件

Show code

lsmfit.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
function output_function = lsmfit(input_functions, x, y)
    %lsmfit - Calculating linear least square method fit curve with given data
    %
    % Syntax: output_function = lsmfit(input_functions, x, y)

    % @Author: Tifa
    % @LastEditTime: 2021-04-22 23:05:47

    input_check(input_functions, x, y)

    len_x = length(x);
    len_funcs = length(input_functions);
    output_function = 0;

    X = zeros(len_x, len_funcs);

    for i = 1:len_funcs
        X(:, i) = subs(input_functions(i), x);
    end

    A = X' * X \ X' * y';

    for i = 1:len_funcs
        output_function = output_function + A(i) * input_functions(i);
    end

end

function input_check(input_functions, x, y)
    % Input check of Exp.2

    % @Author: Tifa
    % @LastEditTime: 2021-04-22 23:05:47

    if ~isvector(input_functions) ||~isvector(x) ||~isvector(y)
        error('All arguments should be vectors')
    end

    if ~isnumeric(x) ||~isnumeric(y)
        error('Input x and y should be numerals')
    end

    if ~isa(input_functions, 'sym')
        error('Input functions should be symbolic expressions')
    end

    if length(x) ~= length(y)
        error('The length of x and y should be equal')
    end

    var1 = symvar(input_functions(1));

    if length(var1) ~= 1
        error('All the input functions should be only one variant')
    end

    for now_func = input_functions(2:end)

        if symvar(now_func) ~= var1
            error('Variant in every input functions should be the same')
        end

    end

end

结果讨论和分析

结果

  • 拟合函数 1:

    \(3.5168482\times10^{-9}x^3 - 5.2948459\times10^{-6}x^2 + 2.6568799\times10^{-4}x\)

    • 误差:

      \(t\)0510152025
      \(\epsilon\)0.000000000.00000695-0.00000026-0.00000527-0.00000372-0.00000124
      \(t\)303540455055
      \(\epsilon\)-0.000000480.000004930.000010350.00001414-0.000042330.00001929

  • 拟合函数 2:

    \(3.4364154\times10^{-9}x^3 - 5.2155622\times10^{-6}x^2 + 2.6339853\times10^{-4}x + 1.7838828\times10^{-5}\)

    • 误差:

      \(t\)0510152025
      \(\epsilon\)-0.000001780.00000613-0.00000046-0.00000513-0.00000345-0.00000100
      \(t\)303540455055
      \(\epsilon\)-0.000000360.000004890.000010180.00001393-0.000042440.00001950

  • 图像

分析

  • 可以发现两个拟合函数图像几近重合
  • 在测试中发现当选用多项式拟合时，次数过小则误差较大，次数较大则对趋势的预测不够合理，且更易受到干扰