九色国产,午夜在线视频,新黄色网址,九九色综合,天天做夜夜做久久做狠狠,天天躁夜夜躁狠狠躁2021a,久久不卡一区二区三区

插值與擬合的Matlab實(shí)現(xiàn)

王 正 盛   編寫

       在科技工程中，除了要進(jìn)行一定的理論分析外，通過(guò)實(shí)驗(yàn)、觀測(cè)數(shù)據(jù)，做分析、處理也是必不可少的一種途徑。由于實(shí)驗(yàn)測(cè)定實(shí)際系統(tǒng)的數(shù)據(jù)具有一定的代表性，因此在處理時(shí)必須充分利用這些信息；又由于測(cè)定過(guò)程中不可避免會(huì)產(chǎn)生誤差，故在分析經(jīng)驗(yàn)公式時(shí)又必須考慮這些誤差的影響。兩者相互制約。據(jù)此合理建立實(shí)際系統(tǒng)數(shù)學(xué)模型的方法成為數(shù)值逼近法。

一、        插值法

1、    數(shù)學(xué)原理

           工程實(shí)踐和科學(xué)實(shí)驗(yàn)中，常常需要從一組實(shí)驗(yàn)觀測(cè)數(shù)據(jù) 中，求自變量 與因變量 的一個(gè)近似的函數(shù)關(guān)系式 。

例如：觀測(cè)行星的運(yùn)動(dòng)，只能得到某時(shí)刻 所對(duì)應(yīng)的行星位置 （用經(jīng)緯度表示），想知道任何時(shí)刻 的行星位置。

例如：大氣壓測(cè)定問(wèn)題；導(dǎo)彈發(fā)射問(wèn)題；程序控制銑床加工精密工件問(wèn)題；飛機(jī)船舶制造問(wèn)題等等。都屬于此類問(wèn)題。因?yàn)榭紤]到代數(shù)多項(xiàng)式既簡(jiǎn)單又便于計(jì)算，所以人們就用代數(shù)多項(xiàng)式近似地表示滿足 個(gè)點(diǎn) 的函數(shù)關(guān)系式 ——插值法建模。

       （1）計(jì)算方法課程中學(xué)習(xí)了兩種多項(xiàng)式插值：Lagrange插值和Newton均差插值：

已知n+1個(gè)數(shù)據(jù)點(diǎn)：

n次Lagrange插值公式：

特別地，當(dāng)n=1時(shí)， ————線性插值

當(dāng)n=2時(shí)，

                                        ———————拋物線插值或二次插值

Newton均差插值公式： ，其中 是k階均差，可由均差表方便計(jì)算得到。

Lagrange插值和Newton均差插值本質(zhì)上是一樣的，只是形式不同而已，因?yàn)椴逯刀囗?xiàng)式是唯一的。

（2）Runge現(xiàn)象和分段低次插值：

如 在[-5，5]上各階導(dǎo)數(shù)存在，但在此區(qū)間取n個(gè)節(jié)點(diǎn)構(gòu)造的Lagrange插值多項(xiàng)式在區(qū)間并非都收斂，而且分散得很厲害。（matlab\bin\ Lagrange.m是自己編寫的M文件）

[例]取 n=-10

hold off  

x=[-5:1:5];

y=1./(1+x.^2);

x0=[-5:0.1:5];

y0=lagrange(x,y,x0);

y1=1./(1+x0.^2);

plot(x0,y0)  

hold on

plot(x0,y1,'b:')   

legend('插值曲線','原數(shù)據(jù)曲線')  

因此插值多項(xiàng)式一般不要超過(guò)四次為宜。為避免Runge現(xiàn)象提出分段插值，所謂分段插值就是首先把插值點(diǎn)分開，在每一段上用低次插值，再連接起來(lái)。

如分段線性插值、分段二次插值、分段三次插值等等。

（3）三次樣條插值：

        三次樣條插值算法的插值精度較高，所構(gòu)造的曲線比較光滑（即在節(jié)點(diǎn)處連續(xù)可導(dǎo)）。因此，在許多工程設(shè)計(jì)或制造業(yè)中，例如飛機(jī)、導(dǎo)彈、高速機(jī)車、萬(wàn)噸級(jí)輪船等外形設(shè)計(jì)常常利用本算法進(jìn)行插值計(jì)算。

三次樣條插值函數(shù)的定義：

設(shè)已知n+1個(gè)數(shù)據(jù) ，如果函數(shù) 滿足如下性質(zhì)：

（1）在每個(gè)子區(qū)間 上， 是次數(shù)不超過(guò)三次的多項(xiàng)式；

（2）

（3） 在插值區(qū)間上二次連續(xù)可微，記為

則稱 為三次樣條插值函數(shù)，簡(jiǎn)稱三次樣條。

     三次樣條函數(shù) 在每個(gè)子區(qū)間 上可由4個(gè)系數(shù)唯一確定，因此， 在 上有4n個(gè)待定系數(shù)。由于 ，

故有

這里給出了3n-3個(gè)條件，加上插值條件（2），共有4n-2個(gè)條件。為了確定 ，通常還需要補(bǔ)充兩個(gè)邊界條件。

常用的邊界條件有三類：

第一種邊界條件：給定兩邊界節(jié)點(diǎn)處的一階導(dǎo)數(shù) ，并要求 滿足：

第二種邊界條件：給定兩邊界節(jié)點(diǎn)處的二階導(dǎo)數(shù) ，并要求 滿足：

特別地，若 ，則所得的樣條稱為自然樣條（MATLAB中即是自然樣條）。

第三種邊界條件：被插函數(shù) 是以 為周期的周期函數(shù)，要求 滿足：

定 理：三次樣條插值問(wèn)題的解存在且唯一。（證明略）                          

求三次樣條函數(shù) 有多種方法，這里介紹一種著名的三彎矩法。            

1、表達(dá)式：設(shè)三次樣條函數(shù)為 ，且 ，在 上， 由定義中的（2）及 處的二階導(dǎo)數(shù)為 ，則有

                                         (1)

                                         (2)

由于 是三次多項(xiàng)式，所以 是線性函數(shù)，于是線性插值，得

令 ，得

             （3）

對(duì)（3）式積分兩次，得

                    （4）

其中 為積分常數(shù)，利用（1），有

解得

代入（4）并化簡(jiǎn)，得

                                                                                 ——————————————————（5）                                                                                                                                                                   這就是系數(shù)用二階導(dǎo)數(shù)表示的三次樣條的表達(dá)式。因此，只要確定 ，即得三次樣條函數(shù) 。

2、   結(jié)點(diǎn)的 關(guān)系式

為了求 ，利用 的一階導(dǎo)數(shù)的連續(xù)性： ；

將（5）式對(duì) 求導(dǎo)數(shù)，得

類似地，可得

由 ，得

上式兩邊同乘以 ，得

令

則得到方程組

           （6）

稱為結(jié)點(diǎn)的M關(guān)系式（力學(xué)上稱為三彎矩方程組）

3、   邊界條件

   （6）式中是n+1個(gè)未知數(shù)的n-1個(gè)方程的方程組，不能唯一的確定 ，要唯一確定 ，必須附加兩個(gè)條件，可由實(shí)際情況在兩端點(diǎn) 處給出的條件，稱為邊界條件（端點(diǎn)條件）。常見兩種邊界條件：

（1）給定兩端點(diǎn)的二階導(dǎo)數(shù)，                               （7）

這時(shí)由結(jié)點(diǎn)的M關(guān)系式（6）可求出 ， 特別地當(dāng) 時(shí)為自然樣條。

（2）給出兩端點(diǎn)的一階導(dǎo)數(shù)（斜率）： ，此時(shí)由一階導(dǎo)數(shù)的連續(xù)性，可得到兩個(gè)方程：

                              （8）

邊界條件（7）和（8）可以寫成統(tǒng)一形式

                                                         （9）

其中

當(dāng) 時(shí)，即為（7）式；當(dāng) 時(shí)， 即為（8）式。

將（6）式與（9）式合在一起得到以下三對(duì)角方程組：

此方程組對(duì)角占優(yōu)，用追趕法一定可解出 。

4、   解題步驟

第一步：確定邊界條件；

第二步：建立結(jié)點(diǎn)的M關(guān)系式，求出各點(diǎn)的 ；

第三步：將所得的 代回樣條函數(shù)表達(dá)式（5），即得 在個(gè)子區(qū)間上的表達(dá)式 。

例題：某飛機(jī)頭部外形曲線數(shù)據(jù)如下：

設(shè)給定端點(diǎn)條件： 時(shí)， ； 時(shí)，

求插值點(diǎn) 處的函數(shù)值及一階、二階導(dǎo)數(shù)值。

答案：

2、    建模實(shí)例

例：丙烷導(dǎo)熱系數(shù)的估計(jì)

通過(guò)實(shí)驗(yàn)得到如下數(shù)據(jù)

表中T、P和K分別表示溫度、壓力和導(dǎo)熱系數(shù)，并假設(shè)在這個(gè)范圍內(nèi)導(dǎo)熱系數(shù)近似地隨壓力線性變化，求P=10130和T=99下的K。

對(duì)于此問(wèn)題進(jìn)行分析，不難發(fā)現(xiàn)這是一個(gè)兩個(gè)自變量的插值問(wèn)題，即要求同時(shí)對(duì)壓力和溫度進(jìn)行內(nèi)插。但是考慮到該問(wèn)題假設(shè)，內(nèi)插將分為二部分完成。首先對(duì)不同溫度求出在P=10130下的一組導(dǎo)熱系數(shù)，利用線性插值求出T=68 和P=10130下的導(dǎo)熱系數(shù)

類似可得到在T=87、106和140時(shí)的導(dǎo)熱系數(shù)：

其次對(duì)溫度的插值，采用三次Lagrange插值多項(xiàng)式計(jì)算得到：

P=10130和T=99下的K=0.0725

3、在MATLAB中的實(shí)現(xiàn)

     MATLAB提供一維、二維和N維插值函數(shù)interp1、 interp2和interpn，以及三次樣條插值函數(shù) spline等。

（1）           一維數(shù)據(jù)插值（一元函數(shù)）

MATLAB提供的函數(shù) interp1可以根據(jù)已知數(shù)據(jù)表[ X,Y]，用各種不同的算法計(jì)算XI各點(diǎn)上的函數(shù)近似值。該函數(shù)有三種調(diào)用格式。

（1.1）YI=interp1(X,Y,XI)

          根據(jù)數(shù)據(jù)表[X,Y]，用分段線性插值算法計(jì)算XI各點(diǎn)（或一個(gè)點(diǎn)）上的函數(shù)近似值YI 。當(dāng)Y是向量時(shí)，則對(duì)Y向量插值，得到結(jié)果YI是與XI同樣大小的向量；當(dāng)Y是矩陣時(shí)，則對(duì)Y的每一列向量插值，得到結(jié)果YI是一矩陣：它的列數(shù)與Y的列數(shù)相同，它的行數(shù)與向量XI的大小相同。

     （1.2）YI=interp1(Y,XI)

          格式調(diào)用方法與（1.1）相同，只是插值節(jié)點(diǎn)用其節(jié)點(diǎn)序號(hào)：X=1:N，這里的N=size(Y),即   X的大小與Y一樣。

      (1.3) YI=interp1(X,Y,Xi,method)

      函數(shù)調(diào)用與上述相同，只是需要指定輸入?yún)?shù) method的具體的算法：

'linear'——分段線性插值，可缺省

'cubic'——分段三次多項(xiàng)式插值

 'spline'——三次樣條插值：即在每個(gè)分段（子區(qū)間）內(nèi)構(gòu)造一個(gè)三次多項(xiàng)式，使其插值函數(shù)滿足插值條件外，還要求在各個(gè)節(jié)點(diǎn)處具有光滑的條件（導(dǎo)數(shù)存在）。

'nearest'——最鄰近區(qū)域插值：即在就近插值節(jié)點(diǎn)的區(qū)域上的函數(shù)取值該節(jié)點(diǎn)之函數(shù)值，該插值函數(shù)為一個(gè)階梯函數(shù)：

例1：已知數(shù)據(jù)表如下，試?yán)?/span>interp1的不同插值算法求xi=1,1.5,2,2.5,…,5各點(diǎn)的函數(shù)的近似值。

解：

  hold off

 xx=1:1:5;

 yx=[3.5,4.6,5.5,3.2,2];

 xxi=1:0.5:5;

f0=interp1(xx,yx,xxi);

f1=interp1(xx,yx,xxi,'linear');

f2=interp1(xx,yx,xxi,'cubic');

f3=interp1(xx,yx,xxi,'spline');

f4=interp1(xx,yx,xxi,'nearest');

f0,f1,f2,f3,f4

f0 =

  Columns 1 through 7

    3.5000    4.0500    4.6000    5.0500    5.5000    4.3500    3.2000

  Columns 8 through 9

    2.6000    2.0000

f1 =

  Columns 1 through 7

    3.5000    4.0500    4.6000    5.0500    5.5000    4.3500    3.2000

  Columns 8 through 9

    2.6000    2.0000

f2 =

  Columns 1 through 7

    3.5000    4.0750    4.6000    5.2625    5.5000    4.4813    3.2000

  Columns 8 through 9

    2.4625    2.0000

f3 =

  Columns 1 through 7

    3.5000    3.7734    4.6000    5.3766    5.5000    4.5953    3.2000

  Columns 8 through 9

    2.0797    2.0000

f4 =

  Columns 1 through 7

    3.5000    4.6000    4.6000    5.5000    5.5000    3.2000    3.2000

  Columns 8 through 9

    2.0000    2.0000  

plot(xxi,f1,xxi,f2,xxi,f3,xxi,f4) 

legend('線性插值','三次插值','樣條插值','最近區(qū)域插值')  

例2：某觀測(cè)站測(cè)得某日6：00~18：00之間每隔2小時(shí)的室內(nèi)外溫度列表如下，要求用三次樣條插值分別求得該日室內(nèi)外6：30~17：30之間每隔2小時(shí)各點(diǎn)的近似溫度。

解：設(shè)時(shí)間變量h為一行向量，溫度t為一個(gè)2列矩陣，第一列存儲(chǔ)室內(nèi)溫度，第二列存儲(chǔ)室外溫度。

h=6:2:18;

t=[18 20 22 25 30 28 24;15 19 24 28 34 32 30]';

xi=6.5:2:17.5;

yi=interp1(h,t,xi,'spline')  

yi =

   18.5020   15.6553

   20.4986   20.3355

   22.5193   24.9089

   26.3775   29.6383

   30.2051   34.2568

   26.8178   30.9594  

plot(h,t,xi,yi)  

例3：編制分段二次插值程序，即在每一個(gè)插值子區(qū)間 上用拋物線插值。

MATLAB\BIN\lag2.m

以例1數(shù)據(jù)為例：

xx=1:1:5;

yx=[3.5 4.6 5.5 3.2 2];

xxi=1:0.5:5; 

f0=interp1(xx,yx,xxi);

f1=interp1(xx,yx,xxi,'linear');

f2=interp1(xx,yx,xxi,'cubic');

f3=interp1(xx,yx,xxi,'spline');

f4=interp1(xx,yx,xxi,'nearest');

f5=lag2(xx,yx,xxi)

plot(xxi,f1,xxi,f2,xxi,f3,xxi,f4,xxi,f5,'b:') 

f5 =

  Columns 1 through 7

    3.5000    4.0750    4.6000    5.4500    5.5000    4.2125    3.2000

  Columns 8 through 9

    2.4625    2.0000

legend('線性插值','三次插值','樣條插值','最近區(qū)域插值','二次插值') 

（2）           三次樣條插值

        三次樣條插值算法的插值精度較高，所構(gòu)造的曲線比較光滑。因此，在許多工程設(shè)計(jì)或制造業(yè)中，例如飛機(jī)、導(dǎo)彈、高速機(jī)車、萬(wàn)噸級(jí)輪船等外形設(shè)計(jì)常常利用本算法進(jìn)行插值計(jì)算。它有兩種調(diào)用格式：

     （2.1）yi=spline(x,y,xi)

本格式根據(jù)數(shù)據(jù)表[x,y]用spline插值計(jì)算各節(jié)點(diǎn)處的函數(shù)近似值。這種格式與使用interp1函數(shù)做三次樣條插值作用一樣。

例：

 xx=1:1:5;

 yx=[3.5,4.6,5.5,3.2,2];

 xxi=1:0.5:5;

 f0=interp1(xx,yx,xxi,)

  f1=interp1(xx,yx,xxi,'spline')

  yi=spline(xx,yx,xxi)  

f0 =

  Columns 1 through 7

    3.5000    4.0500    4.6000    5.0500    5.5000    4.3500    3.2000

  Columns 8 through 9

2.6000                2.0000

f1 =

  Columns 1 through 7

    3.5000    3.7734    4.6000    5.3766    5.5000    4.5953    3.2000

  Columns 8 through 9

    2.0797    2.0000  

yi =

  Columns 1 through 7

    3.5000    3.7734    4.6000    5.3766    5.5000    4.5953    3.2000

  Columns 8 through 9

    2.0797    2.0000  

后兩者計(jì)算結(jié)果完全一樣。

      (2.2) pp=spline(x,y)

            本格式根據(jù)數(shù)據(jù)表[x,y]用spline插值算法獲取三次樣條插值pp形式，即在pp中存儲(chǔ)了各分段的三次樣條插值多項(xiàng)式的系數(shù)和其他有關(guān)必要的信息。欲求得xi的插值結(jié)果yi，應(yīng)利用yi=ppval(pp,xi)格式調(diào)用之。

例：

 xx=1:1:5;

 yx=[3.5,4.6,5.5,3.2,2];

 xxi=1:0.5:5;

 f0=interp1(xx,yx,xxi)

ps=spline(xx,yx);

yyi=ppval(ps,xxi)  

f0 =

  Columns 1 through 7

    3.5000    4.0500    4.6000    5.0500    5.5000    4.3500    3.2000

  Columns 8 through 9

    2.6000    2.0000

yyi =

  Columns 1 through 7

    3.5000    3.7734    4.6000    5.3766    5.5000    4.5953    3.2000

  Columns 8 through 9

    2.0797    2.0000  

注意：兩者計(jì)算結(jié)果不同的。

（3）           二維數(shù)據(jù)插值（二元函數(shù)）

a. 插值基點(diǎn)為網(wǎng)格節(jié)點(diǎn)

二維數(shù)據(jù)插值是構(gòu)造一個(gè)二元插值函數(shù) 去近似 ，即曲面插值。如測(cè)山高水深，畫出更為精確的等高線圖，就需要先插入更多的插值點(diǎn)。與interp1類似，指令interp2可以根據(jù)數(shù)據(jù)表[x,y,z]，用各種不同的算法計(jì)算[xi,yi]各點(diǎn)上的函數(shù)近似值zi。該函數(shù)有兩種調(diào)用格式：

(3.1) zi=interp2(x,y,z,xi,yi)

     本指令格式根據(jù)數(shù)據(jù)表[x,y,z]，用雙線性插值算法計(jì)算坐標(biāo)平面x-y上[xi,yi]各點(diǎn)（或一個(gè)點(diǎn)）的二元函數(shù)近似值zi。這里x可以是一行向量，它與z（矩陣）的各列向量相對(duì)應(yīng)；y可以為一個(gè)列向量，它與z的各行向量相對(duì)應(yīng)。對(duì)于[xi,yi]與zi間的對(duì)應(yīng)關(guān)系，則和[x,y]與z的對(duì)應(yīng)關(guān)系相同。

這里的雙線性插值，是指按x,y的雙向線性插值。

(3.2) i=interp2(z,xi,yi)

  本指令格式與 （3.1）的功能相同，只是插值節(jié)點(diǎn)用其節(jié)點(diǎn)序號(hào)：x=1:n,y=1:m，[m,n]=size(z).

(3.3) zi=interp2(x,y,z,xi,yi,method)

  本指令格式的調(diào)用方法與上述指令格式相同，只是規(guī)定在格式中指定具體的算法。這里輸入?yún)?shù)method有：

'linear' ——雙線性插值。默認(rèn)值，即格式(3.1)

'cubic'——雙三次插值

 'nearest'——最鄰近區(qū)域插值

例：某實(shí)驗(yàn)對(duì)一根長(zhǎng)10 m的鋼軌進(jìn)行熱源的溫度傳播測(cè)試。用x表示測(cè)試點(diǎn)0:2.5:10(m)，用h表示測(cè)試時(shí)間0:30:60(s)，用T表示測(cè)試所得各點(diǎn)的溫度。

試用雙線性插值在一分鐘內(nèi)每隔20s、鋼軌每隔1m處的溫度Ti。

解：

x=0:2.5:10

h=[0:30:60]'

T=[95 14 0 0 0 ;88 48 32 12 6;67 64 54 48 41]

mesh(x,h,T)

xi=[0:10];   %在x,y軸取小的步長(zhǎng)以平滑數(shù)據(jù)。

hi=[0:2:60]';

   Ti=interp2(x,h,T,xi,hi)  

x =

         0    2.5000    5.0000    7.5000   10.0000

h =

     0

    30

    60

T =

    95    14     0     0     0

    88    48    32    12     6

    67    64    54    48    41

Ti =

  Columns 1 through 7

   95.0000   62.6000   30.2000   11.2000    5.6000         0         0

   94.5333   63.2267   31.9200   13.4400    7.7867    2.1333    1.6000

   94.0667   63.8533   33.6400   15.6800    9.9733    4.2667    3.2000

   93.6000   64.4800   35.3600   17.9200   12.1600    6.4000    4.8000

   93.1333   65.1067   37.0800   20.1600   14.3467    8.5333    6.4000

   92.6667   65.7333   38.8000   22.4000   16.5333   10.6667    8.0000

   92.2000   66.3600   40.5200   24.6400   18.7200   12.8000    9.6000

   91.7333   66.9867   42.2400   26.8800   20.9067   14.9333   11.2000

   91.2667   67.6133   43.9600   29.1200   23.0933   17.0667   12.8000

   90.8000   68.2400   45.6800   31.3600   25.2800   19.2000   14.4000

   90.3333   68.8667   47.4000   33.6000   27.4667   21.3333   16.0000

   89.8667   69.4933   49.1200   35.8400   29.6533   23.4667   17.6000

   89.4000   70.1200   50.8400   38.0800   31.8400   25.6000   19.2000

   88.9333   70.7467   52.5600   40.3200   34.0267   27.7333   20.8000

   88.4667   71.3733   54.2800   42.5600   36.2133   29.8667   22.4000

   88.0000   72.0000   56.0000   44.8000   38.4000   32.0000   24.0000

   86.6000   71.5867   56.5733   45.9467   39.7067   33.4667   25.8400

   85.2000   71.1733   57.1467   47.0933   41.0133   34.9333   27.6800

   83.8000   70.7600   57.7200   48.2400   42.3200   36.4000   29.5200

   82.4000   70.3467   58.2933   49.3867   43.6267   37.8667   31.3600

   81.0000   69.9333   58.8667   50.5333   44.9333   39.3333   33.2000

   79.6000   69.5200   59.4400   51.6800   46.2400   40.8000   35.0400

   78.2000   69.1067   60.0133   52.8267   47.5467   42.2667   36.8800

   76.8000   68.6933   60.5867   53.9733   48.8533   43.7333   38.7200

   75.4000   68.2800   61.1600   55.1200   50.1600   45.2000   40.5600

   74.0000   67.8667   61.7333   56.2667   51.4667   46.6667   42.4000

   72.6000   67.4533   62.3067   57.4133   52.7733   48.1333   44.2400

   71.2000   67.0400   62.8800   58.5600   54.0800   49.6000   46.0800

   69.8000   66.6267   63.4533   59.7067   55.3867   51.0667   47.9200

   68.4000   66.2133   64.0267   60.8533   56.6933   52.5333   49.7600

   67.0000   65.8000   64.6000   62.0000   58.0000   54.0000   51.6000

  Columns 8 through 11

         0         0         0         0

    1.0667    0.7200    0.5600    0.4000

    2.1333    1.4400    1.1200    0.8000

    3.2000    2.1600    1.6800    1.2000

    4.2667    2.8800    2.2400    1.6000

    5.3333    3.6000    2.8000    2.0000

    6.4000    4.3200    3.3600    2.4000

    7.4667    5.0400    3.9200    2.8000

    8.5333    5.7600    4.4800    3.2000

    9.6000    6.4800    5.0400    3.6000

   10.6667    7.2000    5.6000    4.0000

   11.7333    7.9200    6.1600    4.4000

   12.8000    8.6400    6.7200    4.8000

   13.8667    9.3600    7.2800    5.2000

   14.9333   10.0800    7.8400    5.6000

   16.0000   10.8000    8.4000    6.0000

   18.2133   13.1867   10.7600    8.3333

   20.4267   15.5733   13.1200   10.6667

   22.6400   17.9600   15.4800   13.0000

   24.8533   20.3467   17.8400   15.3333

   27.0667   22.7333   20.2000   17.6667

   29.2800   25.1200   22.5600   20.0000

   31.4933   27.5067   24.9200   22.3333

   33.7067   29.8933   27.2800   24.6667

   35.9200   32.2800   29.6400   27.0000

   38.1333   34.6667   32.0000   29.3333

   40.3467   37.0533   34.3600   31.6667

   42.5600   39.4400   36.7200   34.0000

   44.7733   41.8267   39.0800   36.3333

   46.9867   44.2133   41.4400   38.6667

   49.2000   46.6000   43.8000   41.0000

mesh(xi,hi,Ti)

雙線性插值得到的鋼軌溫度立體圖

b. 插值基點(diǎn)為散亂節(jié)點(diǎn)

問(wèn)題：已知n個(gè)節(jié)點(diǎn) ，求點(diǎn) 處的插值 。

指令： cz=griddata(x,y,z,cx,cy,'method')

     這里x,y,z都是n維向量，表明數(shù)據(jù)點(diǎn)的橫坐標(biāo)、縱坐標(biāo)和豎坐標(biāo)。向量cx,cy是給定的網(wǎng)格點(diǎn)的橫坐標(biāo)和縱坐標(biāo)，指令cz=griddata(x,y,z,cx,cy,'method')返回在網(wǎng)格(cx,cy)處的函數(shù)值。cx與cy 應(yīng)是方向不同的向量，即一個(gè)是行向量，一個(gè)是列向量。

這里輸入?yún)?shù)method有：

'linear' ——雙線性插值。默認(rèn)值，即格式(3.1)

'cubic'——雙三次插值

'nearest'——最鄰近區(qū)域插值

例：在某海域測(cè)得一些點(diǎn)（x，y）處的水深z，在矩形區(qū)域（75，200）*（-50，150）內(nèi)畫出海底曲面圖形（即設(shè)計(jì)航線問(wèn)題——繪制等高線）

解：

clear

x=[129 140 103.5 88 185.5 195 105 158 108 77 81 162 162 118];

y=[7.5 142 23 147 23 138 86 –6.5 –81 3 57 –66 84 –34];

z=[-4 –8 –6 –8 –6 –8 –8 –9 –9 –8 –8 –9 –4 –9];

cx=75:5:200;

cy=-70:5:150;

cz=griddata(x,y,z,cx,cy','cubic')

meshz(cx,cy,cz)

contour3(cx,cy,cz,16)

contour3(cx,cy,cz,16)

4、    習(xí)題與上機(jī)實(shí)驗(yàn)

1）

h=1:12;

t=[5 8 9 15 25 29 31 30 22 25 27 24];

x1=1:0.2:12;

y=interp1(h,t,x1,'spline');

   plot(h,t,h,t,'+',x1,y,'b:')

   legend('線性插值','原數(shù)據(jù)點(diǎn)','樣條插值')  

線性插值與樣條插值的比較

2）對(duì) 在[-5，5]上用年n(=11)個(gè)節(jié)點(diǎn)（等分）分別作Lagrage插值、分段線性插值和三次樣條插值，用m(=21)個(gè)插值點(diǎn)（等分）作圖，比較結(jié)果。

解：

hold off

n=11;m=21;

x=-5:10/(m-1):5;

y=1./(1+x.^2);

z=0*x;

x0=-5:10/(n-1):5;

y0=1./(1+x0.^2);

y1=lagrange(x0,y0,x);

y2=interp1(x0,y0,x);

y3=interp1(x0,y0,x,'spline');

[x',y',y1',y2',y3']

plot(x,z,'k',x,y,'b',x,y1,'c',x,y2,'k',x,y3,'b:')

legend('','y=1/(1+x^2)','Lagrange','Piece-linear','Spline')

ans =

   -5.0000    0.0385    0.0385    0.0385    0.0385

   -4.5000    0.0471    1.5787    0.0486    0.0484

   -4.0000    0.0588    0.0588    0.0588    0.0588

   -3.5000    0.0755   -0.2262    0.0794    0.0745

   -3.0000    0.1000    0.1000    0.1000    0.1000

   -2.5000    0.1379    0.2538    0.1500    0.1401

   -2.0000    0.2000    0.2000    0.2000    0.2000

   -1.5000    0.3077    0.2353    0.3500    0.2973

   -1.0000    0.5000    0.5000    0.5000    0.5000

   -0.5000    0.8000    0.8434    0.7500    0.8205

         0    1.0000    1.0000    1.0000    1.0000

    0.5000    0.8000    0.8434    0.7500    0.8205

    1.0000    0.5000    0.5000    0.5000    0.5000

    1.5000    0.3077    0.2353    0.3500    0.2973

    2.0000    0.2000    0.2000    0.2000    0.2000

    2.5000    0.1379    0.2538    0.1500    0.1401

    3.0000    0.1000    0.1000    0.1000    0.1000

    3.5000    0.0755   -0.2262    0.0794    0.0745

    4.0000    0.0588    0.0588    0.0588    0.0588

    4.5000    0.0471    1.5787    0.0486    0.0484

    5.0000    0.0385    0.0385    0.0385    0.0385

3) 機(jī)床加工問(wèn)題

     待加工零件的外形根據(jù)工藝要求由一組數(shù)據(jù)(x,y)給出（在平面情況下），用程控銑床加工時(shí)每一刀只能沿x方向和y方向走非常小的一步，這就需要從已知數(shù)據(jù)得到加工所要求的步長(zhǎng)很小的(x,y)的坐標(biāo)。下表給出機(jī)翼截面下輪廓線上的部分?jǐn)?shù)據(jù)，假設(shè)需要得到x坐標(biāo)每改變0.1時(shí)的y的坐標(biāo)。試完成加工所需數(shù)據(jù)，畫出曲線，并求出x=0處的曲線的斜率和13<x<15范圍內(nèi)y的最小值。

機(jī)翼截面下輪廓線上的部分?jǐn)?shù)據(jù)

x0=[0 3 5 7 9 11 12 13 14 15];

y0=[0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.6];

x=0:0.1:15;

y1=lagrange(x0,y0,x);

y2=interp1(x0,y0,x);

y3=interp1(x0,y0,x,'spline');

[x',y1',y2',y3'];

subplot(3,1,1)

plot(x0,y0,'k+',x,y1,'r')

grid

title('Lagrange')

subplot(3,1,2)

plot(x0,y0,'k+',x,y2,'r')

grid

title('Piece-linear')

subplot(3,1,3)

plot(x0,y0,'k+',x,y3,'r')

grid

title('Spline')  

可見Lagrange插值的結(jié)果根本不能用，分段線性插值光滑性較差（特別當(dāng)x=14時(shí)），建立采用三次樣條插值?？梢灾苯訌臄?shù)據(jù)結(jié)果中得到：x=0時(shí)的曲線斜率為 0.0579/0.1=0.579；13<x<15范圍內(nèi)y的最小值為0.9826(x=13.8) 。

二、        曲線（面）擬合法建模

1、 數(shù)學(xué)原理

      已知一組（二維）數(shù)據(jù)，即平面上的 n個(gè)點(diǎn) 互不相同。尋求一個(gè)函數(shù)（曲線） ，使 在某種準(zhǔn)則下與所有數(shù)據(jù)點(diǎn)最為接近，即曲線擬合得最好。首先，確定所求曲線的形式（即經(jīng)驗(yàn)公式），線性最小二乘法是解決曲線擬合最常用的方法。其基本思路：令

其中 是事先選定的一組函數(shù)， 是待定系數(shù)（ 。擬合的準(zhǔn)則是使n個(gè)點(diǎn) 與 的距離 的平方和最小，稱為最小二乘準(zhǔn)則。

記                                

為求 使J達(dá)到最小，只需利用極值的必要條件 得到關(guān)于 的超定線性方程組。

當(dāng) 時(shí)，即為多項(xiàng)式插值。

比較常用的經(jīng)驗(yàn)公式：

多項(xiàng)式： （一般m=2,3，不宜過(guò)高）

特別地，一次多項(xiàng)式擬合又稱為“線性回歸”。

指數(shù)曲線：

雙曲線（一支）：

注意：指數(shù)曲線擬合首先需作變量代換，化成多項(xiàng)式擬合進(jìn)行。

2、 在MATLAB中的實(shí)現(xiàn)

在MATLAB中，實(shí)現(xiàn)最小二乘擬合通常采用兩種途徑：

（1）利用polyfit函數(shù)進(jìn)行多項(xiàng)式擬合。

P=polyfit(x,y,n) ——用 n階多項(xiàng)式擬合x,y向量給定的數(shù)據(jù)

PA=polyval(p,xi)——求xi點(diǎn)上的擬合函數(shù)的近似值。

（2）利用常用的矩陣除法解決復(fù)雜型函數(shù)的擬合。

例1：以一次、二次和三次多項(xiàng)式擬合下數(shù)據(jù)

解：

x=[0.5 1.0 1.5 2.0 2.5 3.0];

y=[1.75 2.45 3.81 4.80 7.00 8.60];

a1=polyfit(x,y,1)

a2=polyfit(x,y,2)

a3=polyfit(x,y,3)

x1=[0.5:0.05:3.0];

y1=a1(2)+a1(1)*x1;

y2=a2(3)+a2(2)*x1+a2(1)*x1.^2;

y3=a3(4)+a3(3)*x1+a3(2)*x1.^2+a3(1)*x1.^3;

  plot(x,y,'*')

  hold on

  plot(x1,y1,'b:',x1,y2,'k',x1,y3,'g')  

a1 =    2.7937   -0.1540

a2 =    0.5614    0.8287    1.1560

a3 =   -0.1156    1.1681   -0.0871    1.5200

  legend('原數(shù)據(jù)圖 ','  一次擬合',' 二次擬合  ','  三次擬合 ')  

p1=polyval(a1,x)  

p1 =

    1.2429    2.6397    4.0366    5.4334    6.8303    8.2271  

p2=polyval(a2,x)

p2 =

    1.7107    2.5461    3.6623    5.0591    6.7367    8.6950  

  p3=polyval(a3,x) 

p3 =

1.7540    2.4855    3.6276    5.0938    6.7974    8.6517

   v1=y-p1;

 v2=y-p2;

 v3=y-p3;

  s1=norm(v1,'fro')

  s2=norm(v2,'fro')

  s3=norm(v3,'fro')  

s1 =

    0.9558

s2 =

    0.4220

s3 =

    0.4057  

可見，三次擬合方差最小。

例2：用6階多項(xiàng)式對(duì)區(qū)間[0，2.5]上的誤差函數(shù) 進(jìn)行最小二乘擬合。

解：

x=0:0.1:2.5;

y=erf(x);   %計(jì)算“誤差函數(shù)”在[0,2.5]內(nèi)的數(shù)據(jù)點(diǎn)

p=polyfit(x,y,6)

px=poly2str(p,'x')

p =

    0.0084   -0.0983    0.4217   -0.7435    0.1471    1.1064    0.0004

px =

   0.0084194 x^6 - 0.0983 x^5 + 0.42174 x^4 - 0.74346 x^3 + 0.1471 x^2

   + 1.1064 x + 0.00044117  

例3      有效擬合的區(qū)間性圖示（用[0,2.5]區(qū)間數(shù)據(jù)擬合曲線擬合[0,5]區(qū)間數(shù)據(jù)）

x=0:0.1:5;

x1=0:0.1:2.5

y=erf(x);

y1=erf(x1);

p=polyfit(x1,y1,6)

f=polyval(p,x);

plot(x,y,'bo',x,f,'r-')

axis([0,5,0,2])

legend('擬合曲線','原數(shù)據(jù)線')  

x1 =

  Columns 1 through 7

         0    0.1000    0.2000    0.3000    0.4000    0.5000    0.6000

  Columns 8 through 14

    0.7000    0.8000    0.9000    1.0000    1.1000    1.2000    1.3000

  Columns 15 through 21

    1.4000    1.5000    1.6000    1.7000    1.8000    1.9000    2.0000

  Columns 22 through 26

    2.1000    2.2000    2.3000    2.4000    2.5000

p =

0.0084   -0.0983    0.4217   -0.7435    0.1471    1.1064    0.0004

說(shuō)明：在[0,2.5]區(qū)間之外，兩條曲線的表現(xiàn)完全不同。

例4：用最小二乘法求一個(gè)形如 的經(jīng)驗(yàn)公式，擬合下數(shù)據(jù)。

下面介紹另一種方法解此擬合問(wèn)題。用矩陣的方法求解，把a,b看成是為知量，已知 ，求解一超定方程：

在MATLAB中的實(shí)現(xiàn)：

x=[19 25 31 38 44];

y=[19.0 32.3 49.0 73.3 97.8];

x1=x.^2

x1=[ones(5,1),x1']

ab=x1\y'

x0=[19:0.2:44];

y0=ab(1)+ab(2)*x0.^2;

clf

plot(x,y,'bo')

hold on   %圖形合在同一坐標(biāo)下

plot(x0,y0,'r-')  

x1 =

         361         625         961        1444        1936

x1 =

           1         361

           1         625

           1         961

           1        1444

           1        1936

ab =

    0.9726

    0.0500

3、  習(xí)題與上機(jī)實(shí)驗(yàn)

1）已知數(shù)據(jù)表如下，試?yán)枚味囗?xiàng)式擬合下數(shù)據(jù)，求出xi=1,1.5,2,2.5,…,5各點(diǎn)的函數(shù)的近似值，并與用interp1的線性插值算法求xi=1,1.5,2,2.5,…,5各點(diǎn)的函數(shù)的近似值作比較。

解：%斷開與上例題“hold on”

 hold off 

xx=1:1:5;

 yx=[3.5,4.6,5.5,3.2,2];

 xxi=1:0.5:5;

 ff0=interp1(xx,yx,xxi)

 yy=polyfit(xx,yx,2);

 ff1=polyval(yy,xxi)

 plot(xx,yx,'*',xxi,ff0,'b:',xxi,ff1,'k')  

ff0 =

  Columns 1 through 7

    3.5000    4.0500    4.6000    5.0500    5.5000    4.3500    3.2000

  Columns 8 through 9

    2.6000    2.0000

ff1 =

  Columns 1 through 7

    3.5257    4.2807    4.7571    4.9550    4.8743    4.5150    3.8771

  Columns 8 through 9

    2.9607    1.7657

legend('線性插值','二次擬合') 

2）分別用線性回歸、二次多項(xiàng)式和十次多項(xiàng)式擬合下數(shù)據(jù)：

 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];

p1=polyfit(x,y,1)

p2=polyfit(x,y,2)

p3=polyfit(x,y,10)

xi=linspace(0,1,100);

z1=polyval(p1,xi);

z2=polyval(p2,xi);

z3=polyval(p3,xi);

plot(x,y,'*',xi,z1,'b:',xi,z2,'g',xi,z3,'k')  

p1 =

   10.3185    1.4400

p2 =

   -9.8108   20.1293   -0.0317

p3 =

  1.0e+006 *

  Columns 1 through 7

   -0.4644    2.2965   -4.8773    5.8233   -4.2948    2.0211   -0.6032

  Columns 8 through 11

    0.1090   -0.0106    0.0004   -0.0000

legend('原數(shù)據(jù)圖','線性擬合','二次擬合','高階擬合')  

可見，高次擬合曲線在數(shù)據(jù)點(diǎn)附近更加接近數(shù)據(jù)點(diǎn)的測(cè)量值，但從整體上來(lái)說(shuō)，曲線波動(dòng)比較大，并不一定適合實(shí)際應(yīng)用的需要。所以，在進(jìn)行高次曲線擬合時(shí)，“越高越好”的觀點(diǎn)是不一定對(duì)的。

                                      此文完成于2001年暑假。

                                           王 正 盛