利用matlab实现非线性拟合（上）

• 0 前言

一般而言，通过已有的数据点去推导其它数据点，常见的方法有插值和拟合。插值适用性较广，尤其是线性插值或样条插值已被广泛的应用。但是通过已知的函数去拟合数据，是连接理论与实验重要的桥梁，这一点是插值无法替代的。

日常学习工作中，经常会遇到下面这种问题：想要用某个具体函数去拟合自己的数据，明明知道这个函数的具体形式，却不知道其中的参数怎么选取。本文就简单介绍一下matlab环境下，如何进行非线性拟合。

由于篇幅有限，本章先以线性拟合为基础，非线性拟合放在下一篇文章中，敬请期待。

• 1 多项式拟合

多项式拟合就是利用下面形式的方程去拟合数据：

matlab中可以用polyfit()函数进行多项式拟合。下面举一个小例子：

对于已有的数据点，我们采用4阶多项式拟合。其中已知函数的的表达式为y=0.03 x^4 - 0.5 x^3 + 2 x^2 - 4，在此基础上添加了一些噪声点。拟合曲线依然采用4阶进行拟合，结果如下。

可以看到拟合曲线与理论曲线基本一致，说明这种方法能够较好的拟合出原始数据的趋势。源代码如下：

x=0:0.5:10;
y=0.03*x.^4-0.5*x.^3+2.0*x.^2-0*x-4+6*(rand(size(x))-0.5);
p=polyfit(x,y,4);
x2=0:0.05:10;
y2=polyval(p,x2);
figure();
subplot(1,2,1)
hold on
plot(x,y,'linewidth',1.5,'MarkerSize',15,'Marker','.','color','k')
plot(x,0.03*x.^4-0.5*x.^3+2.0*x.^2-0*x-4,'linewidth',1,'color','b')
hold off
legend('原始数据点','理论曲线','Location','southoutside','Orientation','horizontal')
legend('boxoff')
box on
subplot(1,2,2)
hold on
plot(x2,y2,'-','linewidth',1.5,'color','r')
plot(x,y,'LineStyle','none','MarkerSize',15,'Marker','.','color','k')
hold off
box on
legend('拟合曲线','数据点','Location','southoutside','Orientation','horizontal')
legend('boxoff')

对于多项式拟合，不是阶数越多越好。比如还是这个上面这个例子，阶数越多，曲线越能够穿过每一个点，但是曲线的形状与理论曲线偏离越大。所以选择最适合的才是最好的。

• 2 线性拟合

线性拟合就是能够把拟合函数写成下面这种形式的：

其中f(x)是关于x的函数，其表达式是已知的。p是常数系数，这个是未知的。

对于这种形式的拟合，matlab内部有一个及其强悍的函数，可以自动输出p的解，并且满足最小二乘。这个函数就是\。没错，就是斜杠。这个符号通常用于求解方程AX=B的情况，我们用X=A\B可以求出未知数X。我们利用当A行和列不等时，输出X的最小二乘这个特性，就可以求出相应的最佳拟合。

还是举个例子

上面这个函数够复杂吧，但是未知数满足线性拟合的要求，所以可以被非常简单的拟合出来。假设a=2.5，b=0.5，c=-1，加入随机扰动。拟合的最终效果为：

最终得到的拟合参数为：a=2.47，b=0.47，c=-0.66。考虑到随机噪声的影响，与原始数据相差不大，源代码如下：

%方法1
x=0:0.5:10;
a=2.5;
b=0.5;
c=-1;
%原函数
y=a*sin(0.2*x.^2+x)+b*sqrt(x+1)+c+0.5*rand(size(x));


%拟合部分
yn1=sin(0.2*x.^2+x);
yn2=sqrt(x+1);
yn3=ones(size(x));
X=[yn1;yn2;yn3];
%拟合，得到系数
Pn=X'\y';
yn=Pn(1)*yn1+Pn(2)*yn1+Pn(3)*yn3;


%绘图
figure()
subplot(1,2,1)
plot(x,y,'linewidth',1.5,'MarkerSize',15,'Marker','.','color','k')
legend('原始数据点','Location','southoutside','Orientation','horizontal')
legend('boxoff')
subplot(1,2,2)
hold on
x2=0:0.01:10;
plot(x2,Pn(1)*sin(0.2*x2.^2+x2)+Pn(2)*sqrt(x2+1)+Pn(3),'-','linewidth',1.5,'color','r')
plot(x,y,'LineStyle','none','MarkerSize',15,'Marker','.','color','k')
hold off
box on
legend('拟合曲线','数据点','Location','southoutside','Orientation','horizontal')
legend('boxoff')

如果细心一点，还可以发现，其实常用的拟合方式中，有很多都是线性拟合，比如多项式拟合，傅里叶拟合等。对于多项式拟合，只需要把拟合的部分替换成x，x.^2，x.^3这种形式。对于傅里叶级数，只需要把拟合的部分替换为sin(x)，sin(2*x)，sin(3*x)这种形式就行。

下面展示一个利用线性拟合，进行不同频率的三角波级数拟合正弦函数的例子：

clear;clc;close allt=0:0.001:2*pi;%原函数YS=sin(t);%基函数N=21;Yo=[];for k=1:N    Yn=sawtooth(k*(t+pi/2),0.5);    Yo=[Yo,Yn'];endYS=YS';%拟合a = Yo\YS;%绘图figure()for k=1:N    clf    plot(t,Yo(:,1:k)*a(1:k,:),t,YS,'LineWidth',1)    ylim([-1.3,1.3])    xlim([0,6.3])    pause(0.1)end

• 3 简单的非线性拟合

最后以一个简单的非线性拟合作为收尾。

如果一个非线性方程，可以化为上面线性方程中公式中给出的样子，那么我们不是也可以套用线性的方法去求解吗？

比如下面的方程：

经过取对数变换，那不就可以直接变为线性的形式了吗

这样求出来之后，再反变换回去，就可以得到原方程的系数了。

clear
clc
close all


%方法1
x=0:0.5:10;
a=2.5;
b=0.5;
%原函数
y=a*exp(-b*x);
y=y+0.5*y.*rand(size(y));%加噪声
%变形函数
%Lg_y=Lg_a+b*(-x)
Lg_y=log(y);
%拟合部分
yn1=ones(size(x));
yn2=-x;
X=[yn1;yn2];
%拟合，得到系数
Pn=X'\Lg_y';
%反变换
a_fit=exp(Pn(1));
b_fit=Pn(2);
%绘图
figure()
hold on
x2=0:0.01:10;
plot(x2,a_fit*exp(-b_fit*x2),'-','linewidth',1.5,'color','r')
plot(x,y,'LineStyle','none','MarkerSize',15,'Marker','.','color','k')
hold off

对于复杂的非线性方程如何求解，考虑到篇幅原因我们放在下集。下集高能，持续关注matlab爱好者公众号，学习matlab编程不迷路。

参考资料：

https://ww2.mathworks.cn/help/matlab/data_analysis/programmatic-fitting.html

图片来源：由 hyhhyh21  在Pixabay上发布