MATLAB 非线性隐函数拟合采坑记录(使用 fsolve solve nlinfit lsqcurvefit函数)
问题描述
MATLAB的 nlinfit 和 lsqcurvefit 函数可以实现显函数的非线性拟合,但是无法直接对隐函数进行拟合。所以,需要配合 fsolve 和 solve 函数建立自定义函数,然后实现隐函数的非线性拟合。
另外,理论上是可以实现多元非线性隐函数的拟合,但是目前我只尝试了形如 f ( x , y ) = 0 f(x,y)=0 f(x,y)=0 的隐函数,后续会对形如 f ( x 1 , x 2 , . . . , x n ) = 0 f(x_1,x_2,...,x_n)=0 f(x1,x2,...,xn)=0 多元函数进行尝试。
解决思路
由于我只是想解决拟合问题,但是并没有现成的数据给我拟合,所以我的解决思路是这样的:
- 人为定义一个形如 f ( x , y ) = 0 f(x,y)=0 f(x,y)=0 的隐函数(为了方便计算说明,这个隐函数也可以写成显函数的形式),随便定义一个函数如下所示。
f ( x 1 , x 2 ) = 1.2 x 1 3 + 3.6 x 2 5 = 0 f(x_1,x_2)=1.2x_1^3+3.6x_2^5=0 f(x1,x2)=1.2x13+3.6x25=0
- 然后手动设置 x 1 x_1 x1的范围,比如 x 1 = [ 0 , 0.01 , 0.02 , 0.03 , . . . , 4 ] x_1 = [0,0.01,0.02,0.03,...,4] x1=[0,0.01,0.02,0.03,...,4],并根据隐函数公式计算出 x 2 x_2 x2的值,然后给 x 1 x_1 x1和 x 2 x_2 x2加上一定的噪音,这样就得到了拟合需要的数据。
上面两步的MATLAB代码如下所示
% x1 与 x2 的值计算
b = [1.2 3.6];
x = zeros(length(0:0.01:4),2);
x(:,1) = (0:0.01:4)';
x(:,2) = - ( b(1)/b(2) .* x(:,1).^3 ) .^ ( 1/5 );
% 加入噪音
A = 0.01;
noise = (2 * rand(size(x)) - 1) * A;
x_noise = x + noise;
- 最后将加入噪音的 x 1 x_1 x1 x 2 x_2 x2 带入 nlinfit 或 lsqcurvefit 函数进行拟合,待拟合参数为 b b b ,预期拟合的结果为 b ( 1 ) = 1.2 , b ( 2 ) = 3.6 b(1)=1.2,b(2)=3.6 b(1)=1.2,b(2)=3.6,即预期的拟合结果应与我先前自己定义的隐函数的系数1.2和3.6相同。
错误示范1
由于在写代码过程中踩了不少的坑,所以现在这里记录一下两个错误的代码,以及发生错误的原因。
代码思路
由于进行拟合的两个函数只能直接拟合显函数,所以一开始我的思路是这样的:
- 令 z = f ( x 1 , x 2 ) z=f(x_1,x_2) z=f(x1,x2) ,把原来的隐函数当做是一个二元显函数;
- 然后令 z = 0 z=0 z=0 ,即函数值 z z z 恒等于0,这样即可输入 x 1 x_1 x1 x 2 x_2 x2和 z z z 的值来实现隐函数的拟合。
这一思路的代码如下:
clear
% 使用nlinfit 和 fsolve函数进行 非线性隐函数拟合的 错误示例1
%% 建立隐函数,隐函数公式为b1*x1^3+b2*x2^5 = 0
% x1 与 x2 的值计算
b = [1.2 3.6];
x = zeros(length(0:0.01:4),2);
x(:,1) = (0:0.01:4)';
x(:,2) = nthroot(- ( b(1)/b(2) .* x(:,1).^3 ), 5 );
% 假设二元函数b1*x1^3+b2*x2^5 = z ,并令z=0
z = zeros(size(x,1),1);
%% 加入噪音(-A到A之间的随机数)
A = 0.01;
noise = (2 * rand(size(x)) - 1) * A;
x_noise = x + noise;
%% 拟合
% 创建待拟合模型的函数句柄
modelfun = @(beta,x)( beta(1).*x(:,1).^3 + beta(2).*x(:,2).^5 );
% 待拟合参数的初始值(初始值的选取会在一定程度上影响拟合结果)
beta0 = [1;1];
% 进行拟合,beta即拟合结果
beta = nlinfit(x_noise,z,modelfun,beta0)
运行几次后的结果如图所示
我们发现拟合的结果并不是1.2和3.6,反而非常的小,几乎等于0,这是为什么呢?
原因解释
这里拟合出的结果,隐函数的两个系数 b ( 1 ) b ( 2 ) b(1)\ b(2) b(1) b(2)基本都等于0,这是因为我把隐函数当成了函数值恒等于0的二元显函数,即 f ( x 1 , x 2 ) = z = 0 f(x_1,x_2)=z=0 f(x1,x2)=z=0。
而在 f ( x 1 , x 2 ) = 1.2 x 1 3 + 3.6 x 2 5 = 0 f(x_1,x_2)=1.2x_1^3+3.6x_2^5=0 f(x1,x2)=1.2x13+3.6x25=0函数中,如果系数全等于0,那么等式仍然成立。
所以这导致了拟合的结果都是0。
模型更正
如果换一个待拟合的模型,那么会出现不同结果。
更正模型1
另外如果待拟合的模型函数是
f
(
x
1
,
x
2
)
=
1.2
x
1
3
+
3.6
x
2
5
+
1
=
0
f(x_1,x_2)=1.2x_1^3+3.6x_2^5+1=0
f(x1,x2)=1.2x13+3.6x25+1=0
那么拟合结果也会不准确,因为当 x 1 x 2 x_1\ x_2 x1 x2的值比较大时, x 1 3 x 2 5 x_1^3\ x_2^5 x13 x25的值远大于1,这就近似成 f ( x 1 , x 2 ) = 1.2 x 1 3 + 3.6 x 2 5 = 0 f(x_1,x_2)=1.2x_1^3+3.6x_2^5=0 f(x1,x2)=1.2x13+3.6x25=0,此时拟合的系数不会为0,但也不会是1.2和3.6。
而是满足 b ( 2 ) b ( 1 ) = 3.6 1.2 = 3 \frac{b(2)}{b(1)} =\frac{3.6}{1.2}=3 b(1)b(2)=1.23.6=3,且具体b(1)和b(2)的数值与常数项的取值有关,例如当常数项为1的时候,求解的结果基本稳定在 b ( 1 ) = 0.8 , b ( 2 ) = 2.4 b(1)=0.8,b(2)=2.4 b(1)=0.8,b(2)=2.4 左右。
更正模型2
如果待拟合的模型函数是
f
(
x
1
,
x
2
)
=
1.2
x
1
3
+
3.6
x
2
5
+
5
x
1
=
0
f(x_1,x_2)=1.2x_1^3+3.6x_2^5+5x_1=0
f(x1,x2)=1.2x13+3.6x25+5x1=0
那么最终拟合的结果会比较准确,即 b ( 1 ) = 1.2 , b ( 2 ) = 3.6 b(1)=1.2,b(2)=3.6 b(1)=1.2,b(2)=3.6
错误示范2
代码思路
之后我使用了fsolve函数来处理隐函数。代码如下
clear
% 使用nlinfit/lsqcurvefit配合fsolve进行拟合,拟合失败案例1
% 拟合失败原因b1*v^3+b2*u^5 = 0中,b1和b2的值的比例只要等于1.2/3.6,等式就成立
% 所以拟合出来的曲线是对的,但是参数并不一定等于1.2和3.6
% 这也是报错"警告: 解处的 Jacobian 矩阵为病态,而且某
% 些模型参数的估计值可能不准确(不可识别)。进行预测时要谨慎。"的原因
%% 原始数据
% 建立隐函数,隐函数公式为b1*u^3+b2*v^5 = 0
b = [1.2 3.6];
u = (0:0.01:4)';
% 使用这种方式开方,否则有可能返回复数值,而我们只需要实数值
v = nthroot(- ( b(1)/b(2) .* u.^3 ) , 5 );
% 加入噪音(-A到A之间的随机数,自变量u不加噪音,应变量v加噪音)
A = 0.01;
noise = (2 * rand(size(v)) - 1) * A;
v_noise = v + noise;
%% 拟合
% nlinfit中,beta是待拟合参数,u(x1)是自变量,函数输出结果v(x2_cal)是应变量,v_noise是fsolove用于求解方程组时候的初始值
modelfun = @(beta,x1)calfun(beta,x1,v_noise);
% 拟合初始值
beta0 = [1;1];
% 使用nlinfit拟合(用另一个拟合的时候就把这个注释掉)
beta = nlinfit(u,v_noise,modelfun,beta0)
% 使用lsqcruvefit拟合(用另一个拟合的时候就把这个注释掉)
% beta = lsqcurvefit(modelfun,beta0,u,v_noise)
%% 结果展示
figure
hold on
% 原始数据
scatter(u,v_noise);
% 拟合后的数据
v_fit = nthroot( -beta(1)/beta(2).*u.^3 , 5 );
plot(u,v_fit,'LineWidth',2);
%% test (这一节用于验证fsolve求解出的结果是否正确,不用时注释)
% b = [1.2 3.6];
% u = (0:0.01:4)';
% v = nthroot(- ( b(1)/b(2) .* u.^3 ) , 5 );
%
% % 加入噪音(-A到A之间的随机数)
% A = 0.01;
% noise = (2 * rand(size(v)) - 1) * A;
% v_noise = v + noise;
%
% % 构建模型函数
% modelfun = @(beta,x1)calfun(beta,x1,v_noise);
%
% % 计算结果
% v_test = modelfun([1.2 3.6],u);
%
% % 通过画图验证结果
% figure
% plot(v_noise);
% figure
% plot(v_test)
%% 由自变量计算隐函数的因变量
% x1 x2 分别为隐函数的两个变量
% 把x1当函数自变量,x2当函数因变量
% beta为需要拟合的参数
% x2_cal 为根据模型计算出的x2的值
% 输入x2的目的是,fsolve求解x2_cal的时候,需要使用x2当求解的初始值
function x2_cal = calfun(beta,x1,x2)
% 预分配
x2_cal = zeros(size(x1));
% 设置fsolve求解器为不显示输出
% 'TolFun'设置成1e-16貌似没啥用
options = optimoptions('fsolve','Display','none','TolFun',1e-16);
% 对每一个输入的u都计算出相应的v_cal
for i = 1:length(x1)
% 设置求解的初始值(fsolve函数从设置的初始值开始找方程的解)
% 把实验得到的v的值当做初始值。
v0 = x2(i);
% 设置求解方程组的函数句柄
fun = @(v)hidefunction(beta,x1(i),v);
% 求解v的值
x2_cal(i) = fsolve(fun,v0,options);
end
end
%% 隐函数模型
% b为需要拟合的参数
% x1和x2为变量
function F = hidefunction(b,x1,x2)
F = b(1).*x1.^3 + b(2).*x2.^5;
end
使用这种方法拟合的结果如下
我们可以从中看到,拟合会发出警告,提示我们雅各比矩阵为病态,并且拟合出来的参数也不对,但是拟合曲线却拟合的很好,这是为什么呢?
原因解释
观察我们拟合的模型,不难发现, b ( 1 ) b ( 2 ) b(1)\ b(2) b(1) b(2) 只要成任意比例都能满足等式成立,所以拟合的时候MATLAB会报警告,因为这个模型没有唯一解。所以返回任意一组解,即使系数值不是1.2 3.6,只要二者比例是对的,最终拟合的曲线就是正确的。
f ( x 1 , x 2 ) = 1.2 x 1 3 + 3.6 x 2 5 = 0 f(x_1,x_2)=1.2x_1^3+3.6x_2^5=0 f(x1,x2)=1.2x13+3.6x25=0
如果换一个待拟合的模型,那么会出现不同结果。
备注:在拟合 f ( x 1 , x 2 ) = 1.2 x 1 3 + 3.6 x 2 5 = 0 f(x_1,x_2)=1.2x_1^3+3.6x_2^5=0 f(x1,x2)=1.2x13+3.6x25=0时,使用 solve 或 lsqcurvefit 的结果是一样的因为系数 b ( 1 ) b ( 2 ) b(1)\ b(2) b(1) b(2) 只要成任意比例都能满足等式成立。
另外,使用solve求解速度非常慢,所以实际操作建议使用fsolve。
模型更正
更正模型1
另外如果待拟合的模型函数是
f
(
x
1
,
x
2
)
=
1.2
x
1
3
+
3.6
x
2
5
+
1
=
0
f(x_1,x_2)=1.2x_1^3+3.6x_2^5+1=0
f(x1,x2)=1.2x13+3.6x25+1=0
此时与上一种不使用fsolve和solve的方法不同,在这种情况下能得到正确的拟合参数,即 b ( 1 ) = 1.2 , b ( 2 ) = 3.6 b(1)=1.2,b(2)=3.6 b(1)=1.2,b(2)=3.6。
更正模型2
如果待拟合的模型函数是
f
(
x
1
,
x
2
)
=
1.2
x
1
3
+
3.6
x
2
5
+
5
x
1
=
0
f(x_1,x_2)=1.2x_1^3+3.6x_2^5+5x_1=0
f(x1,x2)=1.2x13+3.6x25+5x1=0
那么最终拟合的结果也会比较准确,即 b ( 1 ) = 1.2 , b ( 2 ) = 3.6 b(1)=1.2,b(2)=3.6 b(1)=1.2,b(2)=3.6
总结
拟合结果不正确,不一定是对MATLAB的函数使用错误,也可能是模型本身设置的有问题。
综上所述,不使用fsolve和solve函数,所求的拟合结果有局限性,不推荐使用。solve函数的求解时间过长,也不推荐使用。而nlinfit和lsqcurvefit性能目前来看差不多。
所以,推荐使用fsolve配合nlinfit或lsqcurvefit进行非线性隐函数的拟合。