Matlab函数逼近与拟合法

2017/5/9 10
逼近方法： Chebyshev(切比雪夫)逼近：连续函数，多项式。 F= Chebyshev(y,k,x0) Legendre(勒让德)逼近：多项式。 F= Legendre(y,k,x0) Pade(帕德)逼近：有理分式。 F= Pade(y,k,x0) 傅里叶逼近：周期函数，三角多项式。 连续周期函数，[A0,A,B]=FZZ(func,T,n) 离散周期函数，c=DFF(f,N)
i 1
k ( xi ) yi
i 1
m
(k 0,1,, n)
28
上式为由 n+1个方程组成的方程组，称正规方程组。 2017/5/9
引入记号
r (r ( x1),r ( x2 ),,r ( xm ))
f ( y1, y2 ,, ym )
4
9
24个点大致分布在 一条直线附近。
8 7 6
故可认为强度y与 拉伸倍数x的主要 关系应为线性关系：
5 4 3 2 1
y( x) 0 1 x
1
2
3
4
5
6
7
8
9
10
其中0 , 1为待定参数
希望y( x) 0 1 x与所有的数据点 (样本点)(xi , yi ) 越接近越好 .
某些系数未知，可以利用实验方法获得了物 理量之间的关系，通过拟合的方法，求出这 些系数的近似值。 2. 物理规律未知，利用实验方法获得了物理量 之间的关系，通过拟合的方法，得到一个近 似的解析式，用于描述物理规律。
拟合函数尽量靠近数据点如何实现？
2017/5/9 9
2、函数逼近
在区间[a,b]上已知一连续函数f(x)，如果该函 数表达式太过于复杂不利于进行计算机运算，就 会利用一个简单函数去近似f(x)，这就是函数逼近 问题。 如果f(x)的表达式未知，只知道描述f(x)的一 条曲线，这就是曲线拟合问题。 和插值问题不同，逼近和拟合并不要求逼近 函数在已知点上的值一定等于原函数的函数值， 而是按照某种标准使得二者的差值达到最小。
c(1)=int(subs(y,findsym(sym(y)),sym('t'))*T(1)/sqrt(1-t^2),t,-1,1)/pi; c(2)=2*int(subs(y,findsym(sym(y)),sym('t'))*T(2)/sqrt(1-t^2),t,-1,1)/pi; f= c(1)+c(2)*t; 2017/5/9 ……
2017/5/9
25
最小二乘法：
以残差平方和最小问题的解来确定拟合函数的方法。
令
i ( xi ) yi
(i=1,2,…m) －－在回归分析中称为残差
残差向量：
残差向量的各分量平方和记为：
S (a0 , a1 ,, an ) [ a j j ( xi ) yi ]
m
n
m
i
得
[
j 0 i 1
n
m
k
( xi ) j ( xi )]a j k ( xi ) yi
i 1
m
即
m m
(k 0,1,, n)
m i 1 i 1
a0 k ( xi )0 ( xi ) a1 k ( xi )1 ( xi ) an k ( xi ) n ( xi )
科学计算与MATLAB
主讲：唐建国 中南大学材料科学与工程学院 2013.09
第五讲 函数逼近与拟合法
内容提要





2017/5/9
引言 函数逼近 傅里叶逼近 最小二乘法拟合 最小二乘法 多元线性拟合 非线性拟合 MATLAB的拟合函数 小结
3
1、引言
例：考察某种纤维的强度与其拉伸倍数的关系,下表是实 际测定的24个纤维样品的强度与相应的拉伸倍数的记录:
必须找到一种度量标准来衡量什么曲线最接近所有数据点。
2017/5/9 5
15 14 1ቤተ መጻሕፍቲ ባይዱ 12 11 10 9 8
数 拟 插
据 合 值
点 曲 曲
线 线
Y
7 6 5 4 3 2 1 0 0 1 2 3 4 5 6 7 8 9 10 11
X
2017/5/9 6
对这样的数据采用上一讲介绍的插值方法近 似求描述物理规律的解析函数，必然存在下列缺 点：
2017/5/9 15
例：基于神经网络的高炉铁水硅含量预报模型
根据RBF神经网络具有收敛速度快和全局优化的 特点,建立了RBF网络模型,并将其应用对高炉铁水硅 含量预报。监于铁水硅含量与炉缸温度之间的密切 相关性,通过铁水硅含量来间接地反映炉内温度变化。 采用MATLAB中的Newrbe函数进行函数逼近,对高 炉一段连续时期内正常生产的数据的归一化处理后 进行训练和仿真,提高了铁水硅含量预报的命中率。 径向基函数RBF神经网络是一个只有一个隐藏 层的三层前馈神经网络结构。它不仅具有很强的非 线性映射能力,并且具有收敛速度快,全局优化的特点
2017/5/9
12
function f = Chebyshev(y,k,x0) %用切比雪夫多项式逼近已知函数 %已知函数：y
%逼近已知函数所需项数：k
%逼近点的x坐标：x0 %求得的切比雪夫逼近多项式或在x0处
的逼近值：f
syms t; T(1:k+1) = t; T(1) = 1; T(2) = t; c(1:k+1) = 0.0;
① 在一个包含有很多数据点的区间内构造插值
函数，必然使用高次多项式。而高次插值多 项式是不稳定的。 ② 由于数据本身存在误差，利用插值方法得到 的插值多项式必然保留了所有的测量误差， 导致插值函数与物理规律差异较大。
实验数据的拟合可以克服插值方法在处理这 类问题中存在的缺点。
2017/5/9 7
实验数据拟合的基本思想：
i 1 j 0 j k i j i k i i
m
n
a ( x ) ( x ) ( x )y
2017/5/9
m
n
m
i 1 j 0
j
k
i
j
i
i 1
k
i
i
27
由
a ( x ) ( x ) ( x )y
i 1 j 0 j k i j i i 1 k i
y R0 (1 x)
2017/5/9 20
通过实验测得金属铜温度x与电阻y数据如下：
xi(℃) 0 10 20 Yi(Ω) xi(℃) 70 80 90 Yi(Ω) xi(℃) 140 150 160 Yi(Ω)
30 40 50 60
2017/5/9
4.38 4.56 4.70 4.86 5.08 5.24 5.40
使近似函数尽量靠近数据点，而不要求近似 函数一定通过所有数据点。 实验数据拟合可以在一定精度内找出反映物 理量间客观函数关系的解析式。如果实验数据 存在误差，这种做法可以部分抵消原来数据中 的测量误差，从而使所得到的拟合函数更好地 反映物理规律。
2017/5/9
8
利用拟合可以解决两类物理问题：
1. 物理规律已知，但描述物理规律的解析式中
2017/5/9
14
例
N
Y
1
0.8415
2
0.9093
3
0.1411
4
5
6
-0.7568 -0.9589 -0.2794
>> y=[0.8415 0.9093 0.1411 -0.7568 -0.9589 -0.2794]; >> c=DFF(y,6)
c=
Columns 1 through 4 -0.0926 - 0.5003i -0.0260 - 0.0194i -0.0251 + 0.0000i -0.0260 + 0.0194i Columns 5 through 6 -0.0926 + 0.5003i -0.0172 - 0.0000i
100 110 120 130
5.58 5.74 5.96 6.06 6.26 6.44 6.58
170 180 190 200
6.74 6.94 7.12 7.28 7.42 7.60 7.78
21
设一元线性拟合函数为：
Y A0 A1 x
4.38 A0 4.56 A 10 A 0 1 4.70 A0 20 A1 4.86 A0 30 A1 7.78 A0 200 A1
Yi A0 A1 xi
定义偏差：
yi
23
i yi Yi
2017/5/9
yi A0 A1 xi
按照拟合的思想，必须在每一个测量点的偏差都 很小，如何达到这一要求？ 方法一：偏差之和最小
y
i 1
n
n
i
Y ( xi )
但是由于偏差有正有负，求和时可能互相抵消，这并 不能保证在每一个测量点的偏差都很小。 方法二：偏差绝对值之和最小
± à Å º À ­ ì É ± ¶ Ê ý 1 1.9 2 2 3 2.1 4 2.5 5 2.7 6 2.7 7 3.5 8 3.5 9 4 10 4 11 4.5 12 4.6
2017/5/9
xi
¿ Ç
È ¶ 1.4 1.3 1.8 2.5 2.8 2.5 3 2.7 4 3.5 4.2 3.5
2017/5/9 16
2017/5/9
17
2017/5/9
18
2017/5/9
19
3、最小二乘法拟合
3.1 最小二乘法
首先，从一个简单的例子来讨论一元线性拟合与 最小二乘法问题。 通过实验测量， 求金属铜电阻温度系数α，金属 电阻与温度关系如下：
R R0 (1 t )
为了具有一般性，把上式改写为：
y R0 (1 x)
将实验数据代入拟合函数，得到方程组
矛盾方程组
21个线性方程
2017/5/9
22
由于以上矛盾方程组不能确定一组唯一的A0和A1， 也就是说，由方程组可求得A0和A1的多组解，那么究 竟哪一组解最接近客观真实值呢 ? 按照拟合的思想，应当使在每一个测量点拟合函 数的函数值尽量接近测量值，这样的拟合函数才是满 足要求的，即：
y
i 1
i
Y ( xi )
尽管这种方法可以保证在每一个测量点的偏差都很小， 但这种方法数学处理比较困难。
2017/5/9 24
方法三：偏差的平方和最小-----最小二乘法
y
i 1
n
i
Y ( xi )
2
这种方法既可以保证在每一个测量点的偏差都 很小，又方便数学处理，所以这种方法是可行的。
for i=3:k+1 T(i) = 2*t*T(i-1)-T(i-2); c(i) = 2*int(subs(y,findsym(sym(y)),sym('t'))*T(i)/s qrt(1-t^2),t,-1,1)/2; f = f + c(i)*T(i); f = vpa(f,6); if(i==k+1) if(nargin == 3) f = subs(f,'t',x0); else f = vpa(f,6); end end end
2017/5/9
11
Chebyshev(切比雪夫)逼近
当一个连续函数定义在区间[-1，1]上时，可以展开成为切比雪夫 级数 。

f ( x) f nTn ( x)
n 0
T0 ( x) 1, T1 ( x) x Tn 1 ( x) 2 xTn ( x) Tn 1 ( x)
13
离散周期函数的傅里叶逼近
y ci eikx
k 0
n 1
function c=DFF(f,N) %用傅里叶级数逼近已知的离散周期函数 %离散数据点：f %展开项数：N %离散傅里叶逼近系数：c c(1:N)=0; for(m=1:N) for(n=1:N) c(m)=c(m)+f(n)*exp(-i*m*n*2*pi/N); end c(m)=c(m)/N; end
2
2017/5/9
m
n
i 1
j 0


2
26 2
由多元函数求极值的必要条件，有
S 0 ak
m n
(k 0,1,, n)
可得
2 k ( xi )[ a j j ( xi ) yi ] 0
i 1 j 0
即
[ a ( x ) ( x ) ( x ) y ] 0
yi ± à Å º
13 14 15 16 17 18 19 20 21 22 23 24
­ É À ì ± ¶ Ê ý 5 5.2 6 6.3 6.5 7.1 8 8 8.9 9 9.5 10
xi
¿ Ç È ¶ 5.5 5 5.5 6.4 6 5.3 6.5 7 8.5 8 8.1 8.1
yi
纤维强度随拉伸倍数增加而增加。