matlab-函數(shù)逼近與擬合法課件

科學(xué)計(jì)算與MATLAB主講：唐建國(guó)中南大學(xué)材料科學(xué)與工程學(xué)院2013.09第五講函數(shù)逼近與擬合法內(nèi)容提要引言函數(shù)逼近傅里葉逼近最小二乘法擬合最小二乘法多元線性擬合非線性擬合MATLAB的擬合函數(shù)小結(jié)2023/7/213例：考察某種纖維的強(qiáng)度與其拉伸倍數(shù)的關(guān)系,下表是實(shí)際測(cè)定的24個(gè)纖維樣品的強(qiáng)度與相應(yīng)的拉伸倍數(shù)的記錄:纖維強(qiáng)度隨拉伸倍數(shù)增加而增加。1、引言2023/7/21424個(gè)點(diǎn)大致分布在一條直線附近。故可認(rèn)為強(qiáng)度y與拉伸倍數(shù)x的主要關(guān)系應(yīng)為線性關(guān)系：必須找到一種度量標(biāo)準(zhǔn)來(lái)衡量什么曲線最接近所有數(shù)據(jù)點(diǎn)。2023/7/2152023/7/216在一個(gè)包含有很多數(shù)據(jù)點(diǎn)的區(qū)間內(nèi)構(gòu)造插值函數(shù)，必然使用高次多項(xiàng)式。而高次插值多項(xiàng)式是不穩(wěn)定的。由于數(shù)據(jù)本身存在誤差，利用插值方法得到的插值多項(xiàng)式必然保留了所有的測(cè)量誤差，導(dǎo)致插值函數(shù)與物理規(guī)律差異較大。

實(shí)驗(yàn)數(shù)據(jù)的擬合可以克服插值方法在處理這類問(wèn)題中存在的缺點(diǎn)。

對(duì)這樣的數(shù)據(jù)采用上一講介紹的插值方法近似求描述物理規(guī)律的解析函數(shù)，必然存在下列缺點(diǎn)：2023/7/217實(shí)驗(yàn)數(shù)據(jù)擬合的基本思想：

使近似函數(shù)盡量靠近數(shù)據(jù)點(diǎn)，而不要求近似函數(shù)一定通過(guò)所有數(shù)據(jù)點(diǎn)。

實(shí)驗(yàn)數(shù)據(jù)擬合可以在一定精度內(nèi)找出反映物理量間客觀函數(shù)關(guān)系的解析式。如果實(shí)驗(yàn)數(shù)據(jù)存在誤差，這種做法可以部分抵消原來(lái)數(shù)據(jù)中的測(cè)量誤差，從而使所得到的擬合函數(shù)更好地反映物理規(guī)律。2023/7/218利用擬合可以解決兩類物理問(wèn)題：物理規(guī)律已知，但描述物理規(guī)律的解析式中某些系數(shù)未知，可以利用實(shí)驗(yàn)方法獲得了物理量之間的關(guān)系，通過(guò)擬合的方法，求出這些系數(shù)的近似值。物理規(guī)律未知，利用實(shí)驗(yàn)方法獲得了物理量之間的關(guān)系，通過(guò)擬合的方法，得到一個(gè)近似的解析式，用于描述物理規(guī)律。擬合函數(shù)盡量靠近數(shù)據(jù)點(diǎn)如何實(shí)現(xiàn)？2023/7/2192、函數(shù)逼近

在區(qū)間[a,b]上已知一連續(xù)函數(shù)f(x)，如果該函數(shù)表達(dá)式太過(guò)于復(fù)雜不利于進(jìn)行計(jì)算機(jī)運(yùn)算，就會(huì)利用一個(gè)簡(jiǎn)單函數(shù)去近似f(x)，這就是函數(shù)逼近問(wèn)題。如果f(x)的表達(dá)式未知，只知道描述f(x)的一條曲線，這就是曲線擬合問(wèn)題。和插值問(wèn)題不同，逼近和擬合并不要求逼近函數(shù)在已知點(diǎn)上的值一定等于原函數(shù)的函數(shù)值，而是按照某種標(biāo)準(zhǔn)使得二者的差值達(dá)到最小。2023/7/2110逼近方法：Chebyshev(切比雪夫)逼近：連續(xù)函數(shù)，多項(xiàng)式。F=Chebyshev(y,k,x0)Legendre(勒讓德)逼近：多項(xiàng)式。F=Legendre(y,k,x0)Pade(帕德)逼近：有理分式。F=Pade(y,k,x0)傅里葉逼近：周期函數(shù)，三角多項(xiàng)式。連續(xù)周期函數(shù)，[A0,A,B]=FZZ(func,T,n)離散周期函數(shù)，c=DFF(f,N)2023/7/2111Chebyshev(切比雪夫)逼近當(dāng)一個(gè)連續(xù)函數(shù)定義在區(qū)間[-1，1]上時(shí)，可以展開(kāi)成為切比雪夫級(jí)數(shù)。2023/7/2112functionf=Chebyshev(y,k,x0)%用切比雪夫多項(xiàng)式逼近已知函數(shù)%已知函數(shù)：y%逼近已知函數(shù)所需項(xiàng)數(shù)：k%逼近點(diǎn)的x坐標(biāo)：x0%求得的切比雪夫逼近多項(xiàng)式或在x0處的逼近值：fsymst;T(1:k+1)=t;T(1)=1;T(2)=t;c(1:k+1)=0.0;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;……fori=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)/sqrt(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);elsef=vpa(f,6);endendend2023/7/2113離散周期函數(shù)的傅里葉逼近functionc=DFF(f,N)%用傅里葉級(jí)數(shù)逼近已知的離散周期函數(shù)%離散數(shù)據(jù)點(diǎn)：f%展開(kāi)項(xiàng)數(shù)：N%離散傅里葉逼近系數(shù)：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;end2023/7/2114例N123456Y0.84150.90930.1411-0.7568-0.9589-0.2794>>y=[0.84150.90930.1411-0.7568-0.9589-0.2794];>>c=DFF(y,6)c=Columns1through4-0.0926-0.5003i-0.0260-0.0194i-0.0251+0.0000i-0.0260+0.0194iColumns5through6-0.0926+0.5003i-0.0172-0.0000i2023/7/2115例：基于神經(jīng)網(wǎng)絡(luò)的高爐鐵水硅含量預(yù)報(bào)模型

根據(jù)RBF神經(jīng)網(wǎng)絡(luò)具有收斂速度快和全局優(yōu)化的特點(diǎn),建立了RBF網(wǎng)絡(luò)模型,并將其應(yīng)用對(duì)高爐鐵水硅含量預(yù)報(bào)。監(jiān)于鐵水硅含量與爐缸溫度之間的密切相關(guān)性,通過(guò)鐵水硅含量來(lái)間接地反映爐內(nèi)溫度變化。采用MATLAB中的Newrbe函數(shù)進(jìn)行函數(shù)逼近,對(duì)高爐一段連續(xù)時(shí)期內(nèi)正常生產(chǎn)的數(shù)據(jù)的歸一化處理后進(jìn)行訓(xùn)練和仿真,提高了鐵水硅含量預(yù)報(bào)的命中率。徑向基函數(shù)RBF神經(jīng)網(wǎng)絡(luò)是一個(gè)只有一個(gè)隱藏層的三層前饋神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)。它不僅具有很強(qiáng)的非線性映射能力,并且具有收斂速度快,全局優(yōu)化的特點(diǎn)2023/7/21162023/7/21172023/7/21182023/7/21193.1最小二乘法

首先，從一個(gè)簡(jiǎn)單的例子來(lái)討論一元線性擬合與最小二乘法問(wèn)題。為了具有一般性，把上式改寫(xiě)為：

通過(guò)實(shí)驗(yàn)測(cè)量，求金屬銅電阻溫度系數(shù)α，金屬電阻與溫度關(guān)系如下：3、最小二乘法擬合2023/7/2120

通過(guò)實(shí)驗(yàn)測(cè)得金屬銅溫度x與電阻y數(shù)據(jù)如下：xi(℃)Yi(Ω)xi(℃)Yi(Ω)xi(℃)Yi(Ω)04.38705.581406.74104.56805.741506.94204.70905.961607.12304.861006.061707.28405.081106.261807.42505.241206.441907.60605.401306.582007.782023/7/2121

設(shè)一元線性擬合函數(shù)為：將實(shí)驗(yàn)數(shù)據(jù)代入擬合函數(shù)，得到方程組21個(gè)線性方程矛盾方程組2023/7/2122

由于以上矛盾方程組不能確定一組唯一的A0和A1，也就是說(shuō)，由方程組可求得A0和A1的多組解，那么究竟哪一組解最接近客觀真實(shí)值呢?

按照擬合的思想，應(yīng)當(dāng)使在每一個(gè)測(cè)量點(diǎn)擬合函數(shù)的函數(shù)值盡量接近測(cè)量值，這樣的擬合函數(shù)才是滿足要求的，即：定義偏差：2023/7/2123

按照擬合的思想，必須在每一個(gè)測(cè)量點(diǎn)的偏差都很小，如何達(dá)到這一要求？

但是由于偏差有正有負(fù)，求和時(shí)可能互相抵消，這并不能保證在每一個(gè)測(cè)量點(diǎn)的偏差都很小。方法一：偏差之和最小

盡管這種方法可以保證在每一個(gè)測(cè)量點(diǎn)的偏差都很小，但這種方法數(shù)學(xué)處理比較困難。方法二：偏差絕對(duì)值之和最小2023/7/2124

這種方法既可以保證在每一個(gè)測(cè)量點(diǎn)的偏差都很小，又方便數(shù)學(xué)處理，所以這種方法是可行的。方法三：偏差的平方和最小-----最小二乘法2023/7/2125殘差向量的各分量平方和記為：最小二乘法：以殘差平方和最小問(wèn)題的解來(lái)確定擬合函數(shù)的方法。令－－在回歸分析中稱為殘差(i=1,2,…m)殘差向量：2023/7/2126由多元函數(shù)求極值的必要條件，有可得即2023/7/2127上式為由n+1個(gè)方程組成的方程組，稱正規(guī)方程組。由得即2023/7/2128引入記號(hào)則由內(nèi)積的概念可知顯然內(nèi)積滿足交換律正規(guī)方程組便可化為2023/7/2129將其表示成矩陣形式：其系數(shù)矩陣為對(duì)稱陣。所以正規(guī)方程組的系數(shù)矩陣非奇異,即根據(jù)Crame法則,正規(guī)方程組有唯一解，稱其為最小二乘解。2023/7/21302023/7/2131實(shí)現(xiàn)流程圖function[a,b]=LZXEC(x,y)%離散試驗(yàn)數(shù)據(jù)點(diǎn)的線性最小二乘擬合%試驗(yàn)數(shù)據(jù)點(diǎn)的x坐標(biāo)向量：X%試驗(yàn)數(shù)據(jù)點(diǎn)的y坐標(biāo)向量：Y%擬合的一次項(xiàng)系數(shù)：a%擬合的常數(shù)項(xiàng)：bif(length(x)==length(y))n=length(x);elsedisp('x和y的維數(shù)不相等！');return;end%維數(shù)檢查A=zeros(2,2);A(2,2)=n;B=zeros(2,1);fori=1:nA(1,1)=A(1,1)+x(i)*x(i);A(1,2)=A(1,2)+x(i);B(1,1)=B(1,1)+x(i)*y(i);B(2,1)=B(2,1)+y(i);endA(2,1)=A(1,2);s=A\B;a=s(1);b=s(2);2023/7/2132例X12345y1.51.843.45.7>>x=1:5;>>y=[1.51.843.45.7];>>[a,b]=LZXEC(x,y)a=1.0000b=0.28002023/7/2133例.回到引言的實(shí)例，從散點(diǎn)圖可以看出，纖維強(qiáng)度和拉伸倍數(shù)之間近似線性關(guān)系，故可選取線性函數(shù)為擬合函數(shù)建立正規(guī)方程組，其基函數(shù)為根據(jù)內(nèi)積公式，可得正規(guī)方程組為2023/7/2134解得殘差平方和：擬合曲線與散點(diǎn)的關(guān)系如右圖:即為所求的最小二乘解。故2023/7/2135例：金屬銅溫度x與電阻y，線性擬合matlab程序2023/7/21362023/7/2137

線性擬合在物理實(shí)驗(yàn)中應(yīng)用十分廣泛，例如彈性介質(zhì)楊氏模量測(cè)量中應(yīng)變與應(yīng)力的關(guān)系，電阻電路中電流與電壓的關(guān)系等。

有些物理量之間在一定范圍內(nèi)是線性關(guān)系，也可使用線性擬合的方法，只是要注意其適用范圍。

還有一種情況是量物理量之間并不存在線性關(guān)系，但經(jīng)過(guò)適當(dāng)變換后可轉(zhuǎn)化為線性關(guān)系。2023/7/2138

常用的線性變換

函數(shù)變換后的函數(shù)2023/7/21393.2多元線性擬合

影響一個(gè)物理量的因素，很多情況下不止一個(gè)，為了得到描述物理規(guī)律的近似函數(shù)，就必須采取多元線性擬合。

設(shè)物理量y隨x1，x2，…，xk等k個(gè)物理量而變化，即：

為了尋求描述物理規(guī)律的近似函數(shù)，通過(guò)實(shí)驗(yàn)測(cè)得n組數(shù)據(jù)(一般n>k)，利用擬合的方法近似函數(shù)。2023/7/2140測(cè)得n組實(shí)驗(yàn)數(shù)據(jù)如下：i123……nx1ix11x12x13……x1nx2ix21x22x23……x2nx3ix31x32x33……x3n.....xki.....xk1.....xk2.....xk3.....Xknyiy1y2y3……yn2023/7/2141

設(shè)近似函數(shù)為

與一元線性擬合的思路相同，由n組實(shí)驗(yàn)數(shù)據(jù)代入上式，得到n個(gè)方程式構(gòu)成的k元線性方程組，用最小二乘原理確定函數(shù)系數(shù)A0,A1,…,Ak,使yi與Yi的偏差的平方和最小。2023/7/2142偏差為即(i=1,2,…,n)令2023/7/2143求極值2023/7/2144得到2023/7/2145化簡(jiǎn)整理后可得

將第一個(gè)方程式中的A0提出，代入其它各式后，關(guān)于A0,A1,A2,…,Ak的正規(guī)方程組則為2023/7/21462023/7/2147上式中(r,s=1,2,…,k)2023/7/2148

物理學(xué)及各科學(xué)技術(shù)領(lǐng)域都普遍存在非線性函數(shù)關(guān)系，對(duì)此不能直接使用線性擬合的方法，對(duì)這類問(wèn)題可以采取不同的方法解決。方法一：變換為線性擬合有些非線性函數(shù)可以經(jīng)過(guò)變量替換，轉(zhuǎn)化成線性函數(shù)關(guān)系，然后對(duì)替換變量進(jìn)行線性擬合，最后再還原為原始的物理量。但不是所有的函數(shù)關(guān)系都可經(jīng)過(guò)變量替換而轉(zhuǎn)化成線性函數(shù)關(guān)系的。方法二：多項(xiàng)式擬合

任何一個(gè)連續(xù)函數(shù)，在一個(gè)比較小的鄰域內(nèi)均可用多項(xiàng)式任意逼近。所以在物理學(xué)的許多問(wèn)題中，不論物理量直接存在何種函數(shù)關(guān)系，都可用多項(xiàng)式進(jìn)行數(shù)據(jù)擬合。2023/7/2149直線：多項(xiàng)式：有理函數(shù)：指數(shù)：3.3非線性擬合2023/7/2150Log-linear:Gaussians:2023/7/2151

設(shè)有n組實(shí)驗(yàn)數(shù)據(jù)xi，yi，（i=1,2,…,n），用k次多項(xiàng)式擬合，設(shè)擬合方程為：即：偏差為2023/7/2152偏差的平方和為：取極小值即2023/7/2153整理得：即即得到了正則方程組2023/7/2154MATLAB程序?qū)崿F(xiàn)functionA=multifit(X,Y,m)%離散試驗(yàn)數(shù)據(jù)點(diǎn)的多項(xiàng)式曲線擬合%試驗(yàn)數(shù)據(jù)點(diǎn)的x坐標(biāo)向量：X%試驗(yàn)數(shù)據(jù)點(diǎn)的y坐標(biāo)向量：Y%擬合多項(xiàng)式的次數(shù)：m%擬合多項(xiàng)式的系數(shù)向量：AN=length(X);M=length(Y);if(N~=M)disp('數(shù)據(jù)點(diǎn)坐標(biāo)不匹配！');return;endc(1:(2*m+1))=0;b(1:(m+1))=0;。。。。。。forj=1:(2*m+1)%求出c和bfork=1:Nc(j)=c(j)+X(k)^(j-1);if(j<(m+2))b(j)=b(j)+Y(k)*X(k)^(j-1);endendendC(1,:)=c(1:(m+1));fors=2:(m+1)C(s,:)=c(s:(m+s));endA=b'\C;%直接求解法求出擬合系數(shù)2023/7/2155多項(xiàng)式擬合也可轉(zhuǎn)化為多元擬合，只要令一元非線性擬合轉(zhuǎn)化為多元線性擬合解正則方程組2023/7/2156x1

2

3

4y4101826例用多項(xiàng)式擬合函數(shù)：解：設(shè)得即2023/7/2157記系數(shù)矩陣為，則故正規(guī)方程組為解得2023/7/2158擬合曲線：2023/7/2159MATLAB程序解法：>>x=1:4;>>y=[4101826];>>A=multifit(x,y,2)A=0.04890.16120.5672>>m=1:0.01:4;>>n=-1.5+0.49.*m+0.5.*m.^2;>>l=0.0489+0.1612.*m+0.5672.*m.^2;>>plot(x,y,'-*',m,n,'--y',m,l,'-r')>>P=polyfit(x