问题描述

MATLAB的 nlinfit 和 lsqcurvefit 函数可以实现显函数的非线性拟合，但是无法直接对隐函数进行拟合。所以，需要配合 fsolve 和 solve 函数建立自定义函数，然后实现隐函数的非线性拟合。

另外，理论上是可以实现多元非线性隐函数的拟合，但是目前我只尝试了形如 $f (x, y) = 0$ 的隐函数，后续会对形如 $f(x_1,x_2,...,x_n)=0$ 多元函数进行尝试。

解决思路

由于我只是想解决拟合问题，但是并没有现成的数据给我拟合，所以我的解决思路是这样的：

人为定义一个形如 $f (x, y) = 0$ 的隐函数（为了方便计算说明，这个隐函数也可以写成显函数的形式），随便定义一个函数如下所示。

$f(x_1,x_2)=1.2x_1^3+3.6x_2^5=0$

然后手动设置 $x_1$ 的范围，比如 $x_1 = [0,0.01,0.02,0.03,...,4]$ ，并根据隐函数公式计算出 $x_2$ 的值，然后给 $x_1$ 和 $x_2$ 加上一定的噪音，这样就得到了拟合需要的数据。

上面两步的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_2$ 带入 nlinfit 或 lsqcurvefit 函数进行拟合，待拟合参数为 $b$ ，预期拟合的结果为 $b (1) = 1.2, b (2) = 3.6$ ，即预期的拟合结果应与我先前自己定义的隐函数的系数1.2和3.6相同。

错误示范1

由于在写代码过程中踩了不少的坑，所以现在这里记录一下两个错误的代码，以及发生错误的原因。

代码思路

由于进行拟合的两个函数只能直接拟合显函数，所以一开始我的思路是这样的：

令 $z=f(x_1,x_2)$ ，把原来的隐函数当做是一个二元显函数；
然后令 $z = 0$ ，即函数值 $z$ 恒等于0，这样即可输入 $x_1$ $x_2$ 和 $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)$ 基本都等于0，这是因为我把隐函数当成了函数值恒等于0的二元显函数，即 $f(x_1,x_2)=z=0$ 。

而在 $f(x_1,x_2)=1.2x_1^3+3.6x_2^5=0$ 函数中，如果系数全等于0，那么等式仍然成立。

所以这导致了拟合的结果都是0。

模型更正

如果换一个待拟合的模型，那么会出现不同结果。

更正模型1

另外如果待拟合的模型函数是
$f(x_1,x_2)=1.2x_1^3+3.6x_2^5+1=0$

那么拟合结果也会不准确，因为当 $x_1\ x_2$ 的值比较大时， $x_1^3\ x_2^5$ 的值远大于1，这就近似成 $f(x_1,x_2)=1.2x_1^3+3.6x_2^5=0$ ，此时拟合的系数不会为0，但也不会是1.2和3.6。

而是满足 $\frac{b(2)}{b(1)} =\frac{3.6}{1.2}=3$ ，且具体b(1)和b(2)的数值与常数项的取值有关，例如当常数项为1的时候，求解的结果基本稳定在 $b (1) = 0.8, b (2) = 2.4$ 左右。

更正模型2

如果待拟合的模型函数是
$f(x_1,x_2)=1.2x_1^3+3.6x_2^5+5x_1=0$

那么最终拟合的结果会比较准确，即 $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);endend%% 隐函数模型
% b为需要拟合的参数
% x1和x2为变量
function F = hidefunction(b,x1,x2)F = b(1).*x1.^3 + b(2).*x2.^5;end