本文概述
MATLAB将多项式作为行向量执行, 包括按幂次降序排列的系数。例如, 等式P(x)= x4 + 7×3-5x + 9可以表示为-
p = [1 7 0 -5 9];
多项式函数
多元拟合
给定两个向量x和y, 命令a = polyfit(x, y, n)通过数据点(xi, yi)拟合n阶多项式, 并在行中返回(n + 1)个x的幂系数向量系数以x的幂的降序排列, 即, a = [a n a n-1…a 1 a 0]。
多才多艺
给定数据向量x和行向量中多项式的系数, 命令y = polyval(a, x)评估数据点xi处的多项式并生成值yi, 使得
yi = a(1)xn i + a(2)xi(n-1)+ … + a(n)x + a(n + 1)。
因此, 向量a的长度为n + 1, 因此, 求出的多项式的阶次为n。因此, 如果a是5个元素长, 则要求值的多项式将自动确定为四阶。
如果需要误差估计, 则polyfit和polyval均使用最佳参数。
示例:直线(线性)拟合
以下数据是从旨在测量给定弹簧的弹簧常数的实验中收集的。弹簧悬挂着不同的质量m, 并测量了弹簧从其未拉伸状态产生的相应偏移δ。
根据物理学, 我们有F =kδ, 这里F = mg。因此, 我们可以从关系式k = mg /δ中找到k。
在这里, 我们将通过绘制实验数据, 通过数据拟合最佳直线(我们知道δ和F之间的关系是线性的), 然后测量最佳拟合线的斜率来找到k。
m(g)
5.00 10.00 20.00 50.00 100.00
δ(mm)
15.5 33.07 53.39 140.24 301.03
通过数据拟合一条直线意味着我们要找到多项式系数a1和a0(一阶多项式), 以便a1 xi + a0给出yi的”最佳”估计。在步骤中, 我们需要执行以下操作:
步骤1:找出系数ak的值:
a = polyfit(x, y, 1)
步骤2:使用拟合多项式以更精细(更紧密分布)的xj评估y:
y_fitted = polyval(a, x_fine)
步骤3:绘制并查看。将给定输入绘制为点, 将拟合数据绘制为线:
情节(x, y, ‘o’, x_fine, y_fitted);
以下脚本文件显示了通过弹簧实验的给定数据进行直线拟合并找到弹簧常数所包含的所有步骤。
m= [5 10 20 50 100];% mass data (g)
d= [15 .5 33 .07 53 .39 140 .24 301 .03];% displacement data (mm)
g=9.81;% g= 9.81 m/s^2
F= m/1000*g;% compute spring force (N)
a= polyfit (d, F, 1);% fit a line (1st order polynomial)
d_fit=0:10:300; % make a finer grid on d
F_fit=polyval (a, d_fit);% calculate the polynomial at new points
plot (d, F, 'o', d_fit, F_fit)% plot data and the fitted curve
xlabel ('Displacement \delta (mm)'), ylabel ('Force (N)')
k=a(1);% Find the spring constant
text (100, .32, ['\leftarrow Spring Constant K=', num2str(k), 'N/mm']);
并绘制下图:
最小二乘曲线拟合
最小二乘曲线拟合的过程可以简单地在MATLAB中实现, 因为该技术导致需要求解的一组线性方程式。
大多数曲线拟合是多项式曲线拟合或指数曲线拟合(包括幂定律, 例如y = axb)。
两个最常用的功能是:
y = aebx
y = cxd
1.ln(y)= ln(a)+ bx或y = a0 + a1 x, 其中y = ln(y), a1 = b, 而a0 = ln?(a)。
2. ln(y)= ln(c)+ d ln(x)或y = a0 + a1 x, 其中y = ln(y), a1 = d, 而a0 = ln?(c)。
现在我们可以在两种方法中仅使用一阶多项式来使用polyfit来确定未知常数。
涉及的步骤如下:
步骤1:开发新输入:通过获取原始输入的对数, 开发新的输入向量y和x。例如, 要拟合y = aebx类型的曲线, 请创建
ybar = log(y)并保留x并拟合y = cxd类型的曲线, 创建ybar = log(y)和xbar = log(x)。
步骤2:执行线性拟合:使用polyfit找出线性曲线拟合的系数a0和a1。
步骤3:绘制曲线:根据曲线拟合系数, 计算原始常数的值(例如a, b)。根据获得的关系重新计算给定x处的y值, 并绘制曲线和原始数据。
下表显示了真空泵的时间与压力变化读数。我们将通过数据拟合一条曲线P(t)= P0 e-t /τ, 并确定未知常数P0和τ。
t
0 0.5 1.0 5.0 10.0 20.0
P
760 625 528 85 14 0.16
通过记录双方关系的日志, 我们可以
其中P = ln(P), a1 =
-并且a0 = ln?(P0)。因此, 一旦我们拥有a1和a0, 我们就可以轻松计算出P0和τ。
示例:创建脚本文件
% EXPFIT: Exponential curve fit example
% For the following input for (t, p) fit an exponential curve
%p=p0 * exp(-t/tau).
% The problem is resolved by taking a log and then using the linear.
% fit (1st order polynomial)
% original input
t=[0 0.5 1 5 10 20];
p=[760 625 528 85 14 0 .16];
% Prepare new input for linear fit
tbar= t;% no change in t is needed
pbar=log (p);% Fit a 1st order polynomial through (tbar, pbar)
a=polyfit (tbar, pbar, 1);% output is a=a1 a0]
% Evaluate constants p0 and tau
p0 = exp (a(2));% a(2)=a0=log (p0)
tau=-1/a(1);% a1=-1/tau
% (a) Plot the new curve and the input on linear scale
tnew= linspace (0, 20, 20);% generate more refined t
pnew=p0*exp(-tnew/tau); % calculate p at new t
plot (t, p, '0', tnew, pnew), grid
xlabel ('Time' (sec)'), ylabel ('Pressure ((torr)')
% (b) Plot the new curve and the input on semilog scale
lpnew = exp(polyval (a, tnew));
semilogy (t, p, 'o', tnew, lpnew), grid
xlabel ('Time (sec)'), ylabel ('Pressure (torr)')
% Note: We only need one plot; we can select (a) or (b).
并绘制下图: