MATLAB 非线性隐函数拟合采坑记录（使用 fsolve solve nlinfit lsqcurvefit函数）

最新推荐文章于 2023-09-09 19:33:29 发布

EstrangedZ

最新推荐文章于 2023-09-09 19:33:29 发布

阅读量3.5k

点赞数 3

分类专栏： MATLAB 文章标签： matlab 开发语言

本文链接：https://blog.csdn.net/EstrangedZ/article/details/125442199

版权

MATLAB 专栏收录该内容

5 篇文章 5 订阅

订阅专栏

MATLAB 非线性隐函数拟合采坑记录（使用 fsolve solve nlinfit lsqcurvefit函数）

问题描述
解决思路
错误示范1
错误示范2
总结

问题描述

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);

end

end

%% 隐函数模型
% b为需要拟合的参数
% x1和x2为变量
function F = hidefunction(b,x1,x2)

F = b(1).*x1.^3 + b(2).*x2.^5;

end