ATLAB 脚本进行最小二乘法拟合

2744557306

这是一个 MATLAB 脚本，用于进行最小二乘法拟合。脚本首先要求用户输入已知点的 x 和 y 坐标，然后输入拟合的多项式次数 n。脚本使用最小二乘法拟合数据，并绘制了原始数据点和拟合曲线的图表。以下是对代码的主要部分的解释：
function fp = fitpt()
% P6 [; L- r% t. Z0 r! H    % 最小二乘
! P1 z$ }2 W' s! {* N    % 基取 {1, x, ...}
/ Y5 Q6 j$ E: u4 c- `    % fitpt.m

' Y0 I9 A- }/ ?; Q6 o& t6 M1 _    % 默认算例为课本：P65，例3.2
    % x = [0,1,2,3,4,5,6,7]
    % y = [3.95,6.82,9.78,12.91,15.74,19.26,21.73,24.07]
    % 结果：P(x) = 4.005 + 2.936x  平方误差=0.6162
& v- X% j# m+ i: j
! L7 t8 q9 i& t2 y    % MatLab函数：polyfit(x, y, n)
9 o; m: x$ {, [* v
% ^2 b3 |8 M. o1 k% H6 l+ S    s = input('<最小二乘>\n输入已知点的x坐标：（回车表示[0,1,2,3,4,5,6,7]）\n', 's');
    if isempty(s)
/ N7 v+ d" ?8 g: X# W        s = '[0,1,2,3,4,5,6,7]';
    else
7 c3 S& P8 l' _6 V# ?# a        if (s(1) ~= '[')
4 V+ ~5 t/ D3 L0 ~1 f7 A) n! R+ m% B            s = strcat('[', s);
            s = strcat(s, ']');
0 U. _& Y+ K3 z( l        end
    end
! V- _: b' P3 N/ U# ~    x = sym(s);
+ c5 R# u2 [6 {* k0 m; E+ H& f
- h- p. R- ^7 e  J    s = input('输入已知点的y坐标：（回车表示[3.95,6.82,9.78,12.91,15.74,19.26,21.73,24.07]）\n', 's');
0 U! p* Q8 y% ~- k    if isempty(s)
. c$ Q2 `4 }, \1 _3 I9 v        s = '[3.95,6.82,9.78,12.91,15.74,19.26,21.73,24.07]';
* B& k; R7 e9 W. |- l3 j% D    else
+ _. I. ]5 E% T        if (s(1) ~= '[')
            s = strcat('[', s);
            s = strcat(s, ']');
        end
) M" V; x( ^1 j) }    end
$ t2 _5 p6 l- B' j    y = sym(s);
. v0 H$ w  B& ^% F# G    sz = size(x);
    sz = sz(2);
! j/ j' K1 U8 K( N; X    n = input('输入多项式次数n：');
    if (n + 1 > sz)
        n = input('多项式次数需要小于已知点个数，请重新输入n：');
    end
    if (n + 1 > sz)
5 W* ^8 j9 D$ |        error('多项式次数不能小于已知点个数！');
* \, C9 ]: K$ o& C: Q* Q% `    end
- M! q9 C+ u" U+ W: r    fp = s_fitpt_p(x, y, n);

    % 绘制原始数据点和拟合曲线
, @( H4 U7 x9 _- l- r    plot(double(x), double(y), 'r*')
    hold on
    a = double(x(1));
+ P- q+ L' J+ p+ B  P# W    b = double(x(sz));
    x = a:abs(b - 1)/100:b;
+ X3 b. t/ Z& N0 T5 a% @    y = subs(fp, x);
- k; H* S: U* B8 i    plot(x, y)
end
* t3 H, {' M3 K- N- j
/ Y7 W$ d& G* q# d8 O& K3 W  O%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function f = s_fitpt_p(x, y, n)
    % 用 n 次多项式实现的最小二乘法

$ v5 w3 c% U' F) z' q. ]: x& h    sz = size(x);
    sz = sz(2);
1 f% P! T4 T( R5 X    A = zeros(sz, n + 1);
    v = vh(n);
    for i = 1:sz
1 k( e% H- G4 O9 @; z+ ^7 V# i( C; j        A(i, = subs(v, double(x(i)));
! m' }' p4 H) v    end
    f = linsolve(A' * A, A' * y');
+ M% t) m  M) H    f = vpa(f, 4);
/ V3 n: {3 z; Q1 x    f = v * f;
* @# J5 X/ A% }- {end
1 ~4 s7 s- f# \/ Z0 ]- U
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function v = vh(n)
    % Create vector in horizontal style, such as
- B- O- J. f, ^) _    % v = [1, x, x^2, ..., x^n]
6 B* V. ]. J8 K8 i
    if (n < 0 || n > 9)
        error('Make sure ''n'' is in range of [0, 9]')
    end
$ _4 k, F$ ~7 d1 r# X    s = '';
    for i = 0:n
        s = strcat(s, ',x^');
        s = strcat(s, num2str(i));
    end
    s(1) = '[';
4 p4 x' `+ Z  @2 x    sz = size(s);
0 V* \* w; l  _1 V9 F4 X    s(sz(2) + 1) = ']';
0 k2 t9 E' [+ I/ _4 \: r
! \/ J  a5 B, b. L  g, `3 ~/ P    v = simplify(sym(s));
end
8 ^3 X; y8 K; @2 k3 _% x1 ^
% A, E# Y; m) ?这个脚本首先获取用户输入的已知点的 x 和 y 坐标，然后使用最小二乘法进行拟合。最后，脚本绘制了原始数据点和拟合曲线的图表。
% J3 S+ k9 Y# m# K5 c
% Z  S' \! @! L
6 \: p+ v  J' `; {& \" n