在进行方程待定系数求解时,MATLAB提供了多种解决方案。常用 的有:矩阵左除(超定方程求解)、转换为线性回归、曲线拟合、非线性回归等方法。
这里以求解如下方程系数为例:
a
p
(
v
s
,
v
w
,
F
n
)
=
C
v
s
α
v
w
β
(
F
n
R
c
)
γ
/
2
a_p(v_s,v_w,F_n) =C\frac{{v_s}^\alpha}{{v_w}^\beta}({\frac{F_n}{R_c}})^{\gamma/2}
ap(vs,vw,Fn)=Cvwβvsα(RcFn)γ/2
曲线拟合方法
即通过MATLAB拟合工具箱cftool进行拟合。
虽然MATLAB曲线拟合工具箱表面是只能是拟合二元曲线,实际上通过方程的变化可以拟合更多元的曲线,但始终还是存在无法拟合高元曲线的问题。因此针对模型1中的情况,无法使用曲线拟合工具箱。
超定方程组最小二乘/线性方程回归
对模型两边取对数可得将其转变为线性关系:
ln
a
p
=
ln
C
+
α
ln
v
s
−
β
ln
v
w
+
γ
2
ln
F
n
R
c
\ln a_p = \ln C+\alpha \ln v_s -\beta \ln v_w +\frac{\gamma}{2} \ln \frac{F_n}{R_c}
lnap=lnC+αlnvs−βlnvw+2γlnRcFn
由式可见,lnap与lnvs,lnvw,lnFn/Rc成线性关系。可以采用超定方程组最小二乘进行求解其系数,同样也可以采用线性回归的方法求解其系数。通过MATLAB求解超定方程组如下:
%%
%矩阵左除法(最小二乘)
A = [ones(size(vs_fit)) log(vs_fit) -log(vw_fit) log(Fn_fit/Rc)/2];
b = log(ap_fit);
coefs_matrix = A\b %左除,相当于求inv(A)
同样,可以基于线性回归函数regress()
对上式进行回归,具体方法如下:
%%
%曲线拟合方法/线性回归方法
x_curve_fit = [ones(size(vs_fit)) log(vs_fit) -log(vw_fit) log(Fn_fit/Rc)/2];
y_curve_fit = log(ap_fit);
[coefs_curve_fit, confidence_intervals, residual_error, RINT, STATS]...
= regress(y_curve_fit, x_curve_fit);
上述两种方法求得的系数一致。
非线性回归求解系数
MATLAB自带有非线性回归函数,基于迭代的方法求解待定系数。具体实现如下:
%%
%多元非线性回归方法
x_regress_fit = [vs_fit vw_fit Fn_fit/Rc];
y_regress_fit = ap_fit;
my_model = inline('beta(1)*x_regress_fit(:,1).^beta(2)./x_regress_fit(:,1).^beta(3).*x_regress_fit(:,3).^(beta(4)/2)',...
'beta','x_regress_fit'); % 创建内联函数,变量为beta和x_regress_fit
beta0 = [10000 -0.1 0.9 1]; % 设置迭代系数初始值
[beta, r] = nlinfit(x_regress_fit, y_regress_fit, my_model, beta0);
% 运行后警告雅可比矩阵为病态矩阵,且所得结果错误。
采用该方法时,虽然能得到相应结果,但是发现报出警告:雅可比矩阵病态,且得到的结果与实际相差很大。可见,这种方法有一定的局限性。
最终效果如下: