巧妙地进行非线性拟合——非线性拟合转化为线性拟合

1. 前言（不在意来龙去脉的可忽略不看）

对于多项式函数，可以用最小二乘法求得精确的拟合结果，使得拟合函数具有全局最优的拟合误差；对于某些非线性函数，如指数函数y=e^(ax+b)，也可以对函数转化后，求得精确的拟合结果，如上述指数函数可转化为x=(ln y)/a -b/a，同样可以求得具有全局最优拟合误差的拟合函数。上述函数都可以用MATLAB的regress函数或者polyfit函数求得最优的拟合结果，或者可以用广义逆矩阵的最小二乘解计算而得。
但是，对于大多数非线性函数，难以求得拟合误差全局最优的拟合结果，形如y=c1*e^(px)+c2*e^(qx)。这一类函数，一般做法是首先预估拟合函数的参数，以此作为初始值，采用迭代优化的方法求出局部最优的拟合参数，MATLAB中用lsqcurvefit或者nonlinfit函数求出拟合函数，拟合结果的准确性由选取的初始点决定，初始点的选取极为困难。

2. 非线性拟合转化为线性拟合

对于大多数指数函数、三角函数、多项式函数通过四则运算或者复合得到的函数，通常可以用线性微分方程来表示，形如y''+ay'+by=0的常系数线性微分方程则可以表示以下形式的函数：①y=C1*e^(p*x)+C2*e^(q*x)；②y=C1*（1+C2*x）*e^(w*x)；③y=e^(p*x)*（C1*sin (w*x)+C2*cos (w*x)）。
对于上述非线性函数，则可以通过拟合线性微分方程的系数，再计算出对应的非线性函数的参数，并利用样本数据作为初始条件，拟合出函数中的其他常数项。
微分方程需要转化为差分方程：①y'(n)=(y(n+1)-y(n-1))/(x(n+1)-x(n-1))；②y''(n)=(y(n+1)-2*y(n)+y(n-1))/((x(n+1)-x(n-1))/2)^2。
通过拟合方程y''+ay'+by=0则可得到a和b，再通过求解二次方程r^2+a*r+b=0，得到两个解r1和r2。
当r1和r2不相等时，则可得到函数y=C1*e^(r1*x)+C2*e^(r2*x)，再利用样本数据线性拟合则可得到C1和C2。
当r1=r2时，可得到函数y=C1*e^(r1x)+C2*x*e^(r1x)，再利用样本数据线性拟合则可得到C1和C2。
当r1和r2为共轭复数时，可得函数y=e^(a*x)*（C1*sin (b*x)+C2*cos (b*x)），其中a和b分别为r1的实部和虚部，再利用样本数据线性拟合则可得到C1和C2。

事例MATLAB代码如下：

% 产生样本数据
x = linspace(0,98,50)';
y = 0.9*exp(-0.2*x)-0.3*exp(-0.06*x)+10.57;
plot(x,y)
% 计算差分
dyx = (y(3:end) - y(1:end-2))./(x(3:end) - x(1:end-2));
ddyx = (y(3:end) - 2*y(2:end-1) + y(1:end-2))./(x(3:end) - x(1:end-2)).^2*4;
% 拟合微分方程系数
A = [dyx,y(2:end-1),-1*ones(length(dyx),1)];
b = -ddyx;
a = (A.'*A)\A.'*b;
a(3) = a(3)/a(2);
% 求解二次方程
syms r
r0 = double(solve(r^2+a(1)*r+a(2)))
% 线性拟合得到常数项C1,C2
A = [exp(r0(1)*x),exp(r0(2)*x)];
C = (A.'*A)\A.'*(y-a(3))
% 作图
hold on
plot(x,C(1)*exp(r0(1)*x)+C(2)*exp(r0(2)*x)+a(3),'r.')

最后，还可以用lsqcurvefit函数进一步优化拟合结果，采用上述方法得到的参数作为初始点，可以修正差分方程拟合导致的误差。

f = @(b,x)b(1)*exp(b(2)*x)+b(3)*exp(b(4)*x)+b(5);
b0 = [C(1);r0(1);C(2);r0(2);a(3)];
b1 = lsqcurvefit(f,b0,x,y)
hold on
plot(x,f(b1,x),'k.')