一、简介
在科学计算中经常要建立实验数据的数学模型。给定函数的实验数据,需要用比较简单和合适的函数来逼近(拟合)实验数据。
曲线拟合问题的数学描述如下:
已知一组(二维)互不相同的数据 $( x_i , y_i ),i = 1,2,…,n $
寻求一个函数(曲线) y=f(x)y = f (x)y=f(x),使 f(x)f (x)f(x) 在某种准则下与所有数据点最为接近,即曲线拟合得最好。
实际中我们常采用最小二乘法作为上文提到的“某种准则”。
本文主要介绍线性与非线性最小二乘拟合的基本概念与求解方法,其他相关数学原理可结合线性代数知识参考书本介绍,在此不赘述。
二、最小二乘拟合
最小二乘拟合分为线性最小二乘拟合和非线性最小二乘拟合。
1、 线性最小二乘拟合
(1) 数学原理
线性最小二乘拟合是解决曲线拟合问题最常用的方法,基本思路是:
第一步: 先选定一组函数
第二步: 确定用来确定待定系数 a1,a2,…,ama1,a2, … ,ama1,a2,…,am 的准则。如这里我们选取最小二乘准则作为确定待定系数的准则。
最小二乘法:使 nnn 个点 (xi,yi)(x_i,y_i)(xi,yi) 与曲线 y=f(x)y=f(x)y=f(x) 的距离的平方和最小。 即:记J(a1,a2,…,am)=∑i=1nδi2=∑i=1n[f(xi)−yi]2=∑i=1n[∑k=1makrk(xi)−yi]2J(a_1,a_2,…,a_m)=\sum_{i=1}^n\delta_i^2=\sum_{i=1}^n[f(x_i)-y_i]^2\\=\sum_{i=1}^n[\sum_{k=1}^ma_kr_k(x_i)-y_i]^2J(a1,a2,…,am)=i=1∑nδi2=i=1∑n[f(xi)−yi]2=i=1∑n[k=1∑makrk(xi)−yi]2
寻找系数 a1,a2,…ama_1,a_2, …a_ma1,a2,…am 使 J(a1,a2,…am)J(a_1,a_2, …a_m)J(a1,a2,…am) 最小。
(2)求解方法
matlab中已有封装好的函数,用以求解满足上述要求的系数。
a=polyfit(x,y,m)
其中,aaa 为拟合多项式系数,x,yx,yx,y 为输入同长度的数组,mmm 为拟合多项式次数。
2、 非线性最小二乘拟合
(1)数学原理
非线性最小二乘拟合基本方法与线性最小二乘拟合相同,差别在于非线性最小二乘拟合的拟合函数 f(x)f(x)f(x) 为 xxx 的任意非线性函数。
(2)求解方法
用matlab求非线性最小二乘拟合的两个函数为:lsqcurvefit、lsqnonlin,可参考matlab中自带的help文档中的用法。(文末例题中有例子)
三、例题
1、对下面数据作二次多项式拟合
I | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 |
---|---|---|---|---|---|---|---|---|---|---|---|
X_i | 0 | 0.1 | 0.2 | 0.3 | 0.4 | 0.5 | 0.6 | 0.7 | 0.8 | 0.9 | 1 |
Y_i | 0.447 | 1.978 | 3.28 | 6.16 | 7.08 | 7.34 | 7.66 | 9.56 | 9.48 | 9.30 | 11.2 |
解法一:
% 解超定方程的方法
x=0:0.1:1;
y=[-0.447,1.978,3.28,6.16,7.08,7.34,7.66,9.56,9.48,9.30,11.2];
R=[(x.^2)', x', ones(11,1)];
A=R\y'
解法二:
%用多项式拟合的命令
x=0:0.1:1;
y=[-0.447,1.978,3.28,6.16,7.08,7.34,7.66,9.56,9.48,9.30,11.2];
A=polyfit(x,y,2)
z=polyval(A,x);
plot(x,y,'k+',x,z,'r')
拟合结果图像:
2、用下面一组数据拟合c(t)=a+be−0.02ktc(t)=a+be^{-0.02kt}c(t)=a+be−0.02kt 中的参数a,b,ka,b,ka,b,k
Tj | 100 | 200 | 300 | 400 | 500 | 600 | 700 | 800 | 900 | 1000 |
---|---|---|---|---|---|---|---|---|---|---|
Cj*0.001 | 4.54 | 4.99 | 5.35 | 5.65 | 5.90 | 6.10 | 6.26 | 6.39 | 6.50 | 6.59 |
**分析:**该问题即求解最优化问题:
min(12F(a,b,k))=12∑j=110[a+be−0.02ktj−cj]2min(\frac 1 2F(a,b,k))=\frac 1 2 \sum_{j=1}^{10}[a+be^{-0.02kt_j}-c_j]^2min(21F(a,b,k))=21j=1∑10[a+be−0.02ktj−cj]2
**解法一:**用函数lsqcurvefit
**解法二:**用函数lsqnonlin
**注:**两个命令都要先建立M文件定义函数f(x),但定义f(x)的方式不同。