數(shù)學(xué)實(shí)驗(yàn) 課件 第11、12章 插值與擬合實(shí)驗(yàn)、Hill密碼與解密

11.1插值

測量數(shù)據(jù)的數(shù)據(jù)量較小并且數(shù)據(jù)值是準(zhǔn)確的，或者基本沒有誤差，這時(shí)我們一般用插值的方法來解決問題．11.1.1一維插值

已知離散點(diǎn)上的數(shù)據(jù)集

，即已知在點(diǎn)集

上的函數(shù)值

，

構(gòu)造一個解析函數(shù)（其圖形為一曲線）．通過這些點(diǎn)，并能夠求出這些點(diǎn)之間的值，這一過程稱為一維插值．完成這一過程可以有多種方法，MATLAB提供interp1函數(shù)，這個函數(shù)的調(diào)用格式為：

vq=interp1(x,v,xq,method)

該命令用指定的算法找出一個一元函數(shù)v=f(x)，然后以f(x)給出xq處的值.

xq可以是一個標(biāo)量，也可以是一個向量，是向量時(shí)，必須單調(diào)，method可以是下列方法之一：

'linear'：分段線性插值（默認(rèn)方式）；

'nearest'：最鄰近插值；

'cubic'：分段三次Hermite函數(shù)插值;

'spline'：三次樣條插值.例11.1已知某產(chǎn)品從1900年到2010年每隔10年的產(chǎn)量為75.995,91.972,105.711,123.203,131.699,150.697,179.323,203.212,226.505,249.633,256.344,267.893，計(jì)算出1995年的產(chǎn)量，用三次樣條插值的方法，畫出每隔一年的插值曲線圖形，同時(shí)將原始的數(shù)據(jù)畫在同一圖上．解程序如下year=1900:10:2010;product=[75.995,91.972,105.711,123.203,131.699,150.697,179.323,203.212,226.505,249.633,256.344,267.893]p1995=interp1(year,product,1995,'spline')x=1900:2010;y=interp1(year,product,x,'spline');plot(year,product,'o',x,y);計(jì)算結(jié)果為p1995=254.4043，圖形如圖11-1a所示.圖11-1a

三次樣條插值如果用線性插值，則程序的后四行改為：p1995=interp1(year,product,1995,'linear')x=1900:2010;y=interp1(year,product,x,'linear');plot(year,product,'o',x,y);計(jì)算結(jié)果為p1995=252.9885，圖形如圖11-1b所示.圖11-1b

線性插值

這兩種計(jì)算方法得到的數(shù)據(jù)有微小的差異，這種差異我們從兩個圖形上也能夠看到，主要表現(xiàn)在節(jié)點(diǎn)（那些繪制成圓點(diǎn)的點(diǎn)）的附近.前者是光滑的，后者有角點(diǎn)出現(xiàn).例11.2某日測得從零點(diǎn)開始每隔2小時(shí)的環(huán)境溫度數(shù)據(jù)如下（℃）：9，10，10，12，20，26，30，29，27，22，20，14，8請推測早上11點(diǎn)的溫度，并畫出這一天的溫度曲線.解>>x=0:2:24;>>y=[910101220263029272220148];>>y1=interp1(x,y,11,'spline')y1=28.4833可推測出早上11點(diǎn)的溫度是28.48℃.>>xi=linspace(0,24,100);>>yi=interp1(x,y,xi,'spline');>>plot(x,y,'o',xi,yi,'-')得到這一天的溫度曲線，見圖11-2，圓圈為插值點(diǎn)，折線為三次樣條曲線.圖11-2溫度曲線11.1.2二維插值

已知離散點(diǎn)上的數(shù)據(jù)集

，即已知在點(diǎn)集

,上的函數(shù)值

，構(gòu)造一個解析函數(shù)（其圖形為一曲面）.通過這些點(diǎn)，并能夠求出這些已知點(diǎn)以外的點(diǎn)的函數(shù)值，這一過程稱為二維插值.MATLAB利用interp2和griddata函數(shù)進(jìn)行插值.1.二維網(wǎng)格數(shù)據(jù)的插值問題

Zi=interp2(X,Y,Z,Xi,Yi,method)

該命令用指定的算法找出一個二元函數(shù)z=f(x,y)，然后以f(x,y)給出(x,y)處的值，返回?cái)?shù)據(jù)矩陣Zi.Xi，Yi是向量，且必須單調(diào)，Zi和meshgrid(Xi,Yi)是同類型的.method可以用下列方法之一：

'linear'：線性插值（缺省方式）；

'nearest'：最近鄰點(diǎn)插值；

'cubic‘：三次插值；

'spline'：三次樣條函數(shù)插值；例11.3利用插值方法繪制二元函數(shù)

圖像.解假設(shè)僅知其中較少的數(shù)據(jù)，則可以由下面的命令繪制出已知數(shù)據(jù)的網(wǎng)格圖，如圖11-3所示.[x,y]=meshgrid(-3:.6:3,-2:.4:2);z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y);%生成樣本點(diǎn)surf(x,y,z)axis([-33,-22,-1,1.5])[x1,y1]=meshgrid(-3:.2:3,-2:.2:2);%選擇更密集的插值點(diǎn)z1=interp2(x,y,z,x1,y1);%線性插值figure(2)surf(x1,y1,z1)axis([-33,-22,-1,1.5])z2=interp2(x,y,z,x1,y1,'nearest');%最近鄰點(diǎn)插值figure(3)surf(x1,y1,z2)axis([-33,-22,-1,1.5])z3=interp2(x,y,z,x1,y1,'cubic');figure(4)surf(x1,y1,z3)axis([-33,-22,-1,1.5])z4=interp2(x,y,z,x1,y1,'spline');figure(5)surf(x1,y1,z4)axis([-33,-22,-1,1.5])adcb圖11-3二元函數(shù)插值圖像

從圖11-3a和11-3b可以看出，由這些數(shù)據(jù)繪制的圖形很粗糙的.圖11-3c和11-3d比較理想.例11.4已知1950年到1990年每隔10年，服務(wù)年限從10年到30年每隔10年的勞動報(bào)酬表如表11-1.服務(wù)年限年份1020301950150.697169.592187.6521960179.323195.072250.2871970203.212239.092322.7671980226.505273.706426.7301990249.633370.281598.243表11-1某企業(yè)工作人員的月平均工資（元）試計(jì)算1975年時(shí)，15年工齡的工作人員平均工資.解程序如下：years=1950:10:1990;service=10:10:30;wage=[150.697169.592187.652179.323195.072250.287203.212239.092322.767226.505273.706426.730249.633370.281598.243]w=interp2(service,years,wage,15,1975)計(jì)算結(jié)果為235.6287.可知1975年時(shí)，15年工齡的工作人員平均工資為235.6287元.例11.5設(shè)有數(shù)據(jù)x=1,2,3,4,5,6，y=1,2,3,4，在由x，y構(gòu)成的網(wǎng)格上，數(shù)據(jù)為12,10,11,11,13,1516,22,28,35,27,2018,21,26,32,28,2520,25,30,33,32,20畫出原始網(wǎng)格圖和將網(wǎng)格細(xì)化為間隔為0.1后的插值網(wǎng)格圖.解程序如下：x=1:6;y=1:4;t=[12,10,11,11,13,15;16,22,28,35,27,20;18,21,26,32,28,25;20,25,30,33,32,20];subplot(1,2,1)mesh(x,y,t)x1=1:0.1:6;y1=1:0.1:4;[x2,y2]=meshgrid(x1,y1);t1=interp2(x,y,t,x2,y2,'cubic');subplot(1,2,2)mesh(x1,y1,t1);結(jié)果如圖11-4所示.11-4a是給定的網(wǎng)格處的數(shù)據(jù)，11-4b是插值后的數(shù)據(jù).ab圖11-4插值網(wǎng)格圖2.二維散點(diǎn)分布數(shù)據(jù)的插值問題

通過上面的例子可以看出，MATLAB提供的二維插值函數(shù)還是能較好地進(jìn)行二維插值運(yùn)算的.但該函數(shù)有一個重要的缺陷，就是它只能處理以網(wǎng)格形式給出的數(shù)據(jù)，如果已知數(shù)據(jù)不是以網(wǎng)格形式給出的，則用該函數(shù)是無能為力的.在實(shí)際應(yīng)用中，大部分問題都是以實(shí)測的

散點(diǎn)給出的，所以不能直接使用函數(shù)interp2進(jìn)行二維插值.

MATLAB語言提供了一個更一般的griddata函數(shù)，用來專門解決這樣的問題，調(diào)用格式如下：

vq

=griddata(x,y,v,xq,yq,method)

該命令用以處理插值基點(diǎn)為散亂節(jié)點(diǎn)的插值問題，變量含義與interp2相同，但不要求x、y分量數(shù)值單調(diào)，所用插值方法也有所不同，method可以用下列方法'之一：

‘linear’：基于三角剖分的線性插值（默認(rèn)）；

‘nearest’：基于三角剖分的最鄰近插值；

‘natural’：基于三角剖分的自然鄰近插值；

‘cubic’：基于三角剖分的三次插值；

‘v4‘：雙調(diào)和樣條插值.例11.6已知某處山區(qū)地形選點(diǎn)測量坐標(biāo)數(shù)據(jù)為：x=0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5y=0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

5.5

6海拔高度數(shù)據(jù)(單位：米)為：z=8990878592919693908782929698999591898684828496989592908885848381858081828995969392898686828587989996978885828382858994959392918684888892939495898786838192929697989693958482818485858182808081859093958486819899989796958487808185828384879095868880828184858683828180828788899899979698949287請畫出地形地貌圖和等高線圖.解程序如下：x=0:.5:5;y=0:.5:6;[xx,yy]=meshgrid(x,y);z=[8990878592919693908782;9296989995918986848284;9698959290888584838185;8081828995969392898686;8285879899969788858283;8285899495939291868488;8892939495898786838192;9296979896939584828184;8585818280808185909395;8486819899989796958487;8081858283848790958688;8082818485868382818082;8788899899979698949287];mesh(xx,yy,z)%原始地貌圖（見圖11-5a）xi=linspace(0,5,50);%加密橫坐標(biāo)數(shù)據(jù)到50個yi=linspace(0,6,80);%加密縱坐標(biāo)數(shù)據(jù)到80個[xii,yii]=meshgrid(xi,yi);%生成網(wǎng)格數(shù)據(jù)zii=interp2(x,y,z,xii,yii,'cubic');%插值figure;mesh(xii,yii,zii)%加密后的地貌圖（見圖11-5b）ab圖11-5山區(qū)地形地貌圖例11.7（海底曲面圖）某海域測得一些點(diǎn)（x,y)處的水深z由表11-2給出，在矩形區(qū)域(75,200)×(-50,150)內(nèi)畫出海底曲面的圖形.x129140103.588185.5195105y7.5141.52314722.5137.585.5z4868688x157.5107.57781162162117.5y-6.5-81356.5-66.584-33.5z9988949表11-2海底測量數(shù)據(jù)解程序如下：x=[129 140

103.588185.5

195 105

157.5107.5 77

81

162

162 117.5];y=[7.5 141.52314722.5

137.585.5

-6.5

-81

3 56.5-66.584 -33.5];z=-[4 8

6

8 6

8

8 99

8

8

9

4

9];plot3(x,y,z,'o');holdon%原始數(shù)據(jù)點(diǎn)cx=75:0.5:200;cy=-70:0.5:150;cz=griddata(x,y,z,cx,cy','v4');%插值mesh(cx,cy,cz)結(jié)果如圖11-6所示.圖11-6海底曲面圖11.2擬合

測量數(shù)據(jù)的數(shù)據(jù)量較大或者測量值與真實(shí)值有誤差，這時(shí)一般用曲線擬合的方法來解決問題.

已知離散點(diǎn)上的數(shù)據(jù)集

，即已知在點(diǎn)集

上的函數(shù)值

，構(gòu)造一個解析函數(shù)f(x)（其圖形為一曲線），使f(x)在原離散點(diǎn)xi上的值盡可能接近給定的yi值，這一構(gòu)造函數(shù)f(x)的過程稱為曲線擬合.最常用的曲線擬合方法是最小二乘法，該方法是尋找函數(shù)f(x)使得

最小.

在MATLAB中，利用polyfit和lsqcurvefit函數(shù)進(jìn)行擬合.11.2.1多項(xiàng)式擬合polyfit函數(shù)調(diào)用格式如下：p=polyfit(x,y,n)其中x和y是節(jié)點(diǎn)數(shù)據(jù)向量，n是多項(xiàng)式的次數(shù)；p是返回的多項(xiàng)式p(x)=p1xn+p2xn?1+...+pnx+pn+1系數(shù)向量，系數(shù)按降冪排列，p的長度為n+1.例11.8求如下給定數(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].解MATLAB命令如下：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];p=polyfit(x,y,2)%選取次數(shù)的時(shí)候，次數(shù)長度不得超過length(x)-1x1=0.5:0.05:3.0;y1=polyval(p,x1);%求以p為系數(shù)的多項(xiàng)式的值plot(x,y,'*r',x1,y1,'-b')計(jì)算結(jié)果為p=0.56140.82871.1560此結(jié)果表示擬合函數(shù)為：用此函數(shù)擬合數(shù)據(jù)的效果如圖11-7所示.圖11-7多項(xiàng)式擬合例11.9某鄉(xiāng)鎮(zhèn)企業(yè)2015-2021年的生產(chǎn)利潤如表11-3.

試預(yù)測2022年的利潤.

解作已知數(shù)據(jù)的散點(diǎn)圖，>>x0=[2014201520162017201820192020];>>y0=[70122144152174196202];>>plot(x0,y0,'*')發(fā)現(xiàn)該鄉(xiāng)鎮(zhèn)企業(yè)的年生產(chǎn)利潤幾乎直線上升（見圖11-8）.因此，我們可以y=a1x+a0用作為擬合函數(shù)來預(yù)測該鄉(xiāng)鎮(zhèn)企業(yè)未來的利潤.年份2015201620172018201920202021利潤（萬元）70122144152174196202表11-3生產(chǎn)利潤

圖11-8數(shù)據(jù)的散點(diǎn)圖11.2.2非線性擬合

lsqcurvefit函數(shù)是MATLAB提供的用于求解非線性最小二乘擬合問題的專用函數(shù)，用作各種類型曲線的擬合，調(diào)用格式如下：

x=lsqcurvefit(fun,x0,xdata,ydata)

其中，fun是非線性函數(shù)模型；x0是求解初始值；xdata和ydata是原始輸入和輸出數(shù)據(jù)向量；返回值x為擬合輸出的未知參數(shù).例11.10用最小二乘法求形如

的函數(shù)擬合下列數(shù)據(jù).x0.10.20.150-0.20.3y0.950.840.861.061.500.72表11-4擬合數(shù)據(jù)解先將參數(shù)，合寫為，編寫如下程序：fun=@(c,x)c(1)*exp(c(2)*x)x=[0.1,0.2,0.15,0,-0.2,0.3];y=[0.95,0.84,0.86,1.06,1.50,0.72];c=lsqcurvefit(fun,[0,0],x,y)norm(feval(fun,c,x)-y)^2其中，[0,0]是初始值，最后一句是計(jì)算殘差的平方和，也就是擬合函數(shù)在給定點(diǎn)的值和原始數(shù)據(jù)的差的平方和，運(yùn)行結(jié)果為：Optimizationterminated:relativefunctionvaluechangingbylessthanOPTIONS.TolFun.c=1.0997-1.4923ans=0.0031可得a=1.0997，b=-1.4923，擬合函數(shù)為

.ans=0.0031說明殘差很小.例11.11我國從1954年到2005年全國總?cè)丝诘臄?shù)據(jù)如表11-5，請預(yù)測2019年和2030年的全國人口數(shù).表11-5各年份全國總?cè)藬?shù)（單位：千萬）解因?yàn)長ogistic人口模型的基本形式為其中t0=1954，N0=60.2.當(dāng)t→+∞時(shí)，人口數(shù)量達(dá)到最大值Nm.首先定義函數(shù)文件：functionf=flsqcurvefit(x,t)f=x(1)./(1+(x(1)/60.2-1).*exp(-x(2).*(t-1954)));MATLAB程序如下：cleart0=1954:2005;y0=[60.261.562.864.666.067.266.265.967.369.170.472.5...74.576.3 78.580.783.0 85.287.1 89.2 90.992.4 93.7 95.0...96.259 97.598.705100.1 101.654 103.008 104.357 105.851 107.5...109.3 111.026 112.704114.333 115.823 117.171 118.517 119.850...121.121122.389123.626124.761125.786 126.743 127.627...128.453129.277129.988130.756];x=[0.03,250];a=lsqcurvefit('flsqcurvefit',x,t0,y0)plot(t0,y0,'o')holdonlt=1954:0.1:2040;ly=flsqcurvefit(a,lt);plot(lt,ly)t1=[2019,2030];y1=flsqcurvefit(a,t1)holdoff運(yùn)行后可得：Localminimumpossible.lsqcurvefitstoppedbecausethefinalchangeinthesumofsquaresrelativetoitsinitialvalueislessthanthevalueofthefunctiontolerance.<stoppingcriteriadetails>a=0.0336181.0280y1=147.6920156.6069得到圖形如圖11-9所示.

從圖11-9可以看出，擬合效果較好，2019年全國總?cè)藬?shù)為147.6920千萬人，接近實(shí)際人口數(shù)140.005千萬人，并預(yù)計(jì)2030年全國總?cè)藬?shù)為156.6069千萬人.圖11-9人口增長預(yù)測圖Matlab工具箱提供了命令cftool，該命令給出了一維數(shù)據(jù)擬合的交互式環(huán)境.具體執(zhí)行步驟如下：

（1）把數(shù)據(jù)導(dǎo)入到工作空間；

（2）運(yùn)行cftool，打開用戶圖形界面窗口；

（3）對數(shù)據(jù)進(jìn)行預(yù)處理；

（4）選擇適當(dāng)?shù)哪Ｐ瓦M(jìn)行擬合；

（5）生成一些相關(guān)的統(tǒng)計(jì)量，并進(jìn)行預(yù)測.

可以通過幫助（運(yùn)行doccftool）熟悉該命令的使用細(xì)節(jié).

12.1Hill2密碼加密

一般的加密過程是這樣的：

明文→加密器→密文→普通信道→解密器→明文

其中的“→普通信道→解密器”這個環(huán)節(jié)容易被敵方截獲并加以分析.

在這個過程中，運(yùn)用的數(shù)學(xué)手段是矩陣運(yùn)算，加密過程的具體步驟如下：ABCDEFGHIJKLM12345678910111213NOPQRSTUVWXYZ1415161718192021222324250表12-1明文字母的表值

1)根據(jù)明文字母的表值，將明文信息用數(shù)字表示，設(shè)明文信息只需要26個拼音大寫字母A—Z(也可以不止26個，如還有小寫字母、數(shù)字、標(biāo)點(diǎn)符號等)，通信雙方給出這26個字母表值(見表12-1).

2)選擇一個二階可逆整數(shù)方陣A，稱為Hill2密碼的加密矩陣，它是這個加密體制的“密鑰”(是加密的關(guān)鍵，僅通信雙方掌握).

3)將明文字母依次逐對分組.Hill2密碼的加密矩陣為二階矩陣，則明文字母每2個一組(可以推廣至Hilln密碼，則每n個明文字母為一組).若最后一組僅有一個字母，則補(bǔ)充一個沒有實(shí)際意義的啞字母，這樣使每一組都由2個明文字母組成.查出每個明文字母的表值，構(gòu)成一個二維列向量α.

4)A乘以α，得一新的2維列向量β=Aα，由β的兩個分量反查字母表值得到的兩個字母即為密文字母.

以上4步即為Hill2密碼的加密過程.解密過程，即為上述過程的逆過程.例12.1明文為WAWDZG(“我愛我的祖國”的拼音縮寫)，

，求這段明文的Hill2密文.解將明文相鄰字母每2個分為一組：

WAWDZG(12-1)查表12-1得到每對字母的表值，并構(gòu)造2維列向量：

(12-2)將上述3個向量左乘矩陣，得到3個2維列向量：

(12-3)作模26運(yùn)算(每個元素都加減26的整數(shù)倍，使其化為0～25之間的一個整數(shù))得到：反查表12-1得到每對表值對應(yīng)的字母為：

YCELNU(12-4)這就得到了密文“YCELNU”.12.2Hill2密碼破譯

要將例12.1密文解密，只要將上述加密過程逆轉(zhuǎn)回去，即將密文按同樣方式分組，查它們的表值即得：

(12-5)

(12-5)是前面的(12-3)經(jīng)模26運(yùn)算的結(jié)果.但如何由式(12-5)中的向量求得(12-2)中的向量呢？這是在模運(yùn)算意義下，如何解方程組：

Aα=β

(12-6)的問題.一個一般的n階方陣可逆的充要條件為det(A)≠0.但在模26意義下矩陣可逆與一般的矩陣可逆有所不同.記整數(shù)集合

，m為一正整數(shù)，模m可逆定義如下：定義12.1

對于一個元素屬于集合zm的n階方陣A，若存在一個元素屬于集合zm的方陣B，使得稱A為模m可逆，B為A的模m逆矩陣，記為B=A-1(modm).

的意義是，每一個元素減去m的整數(shù)倍后，可以化成單位矩陣.

例如：

定義12.2

對于Zm的一個整數(shù)a，若存在Zm中的中一個整數(shù)

b，使得

，則稱b為

a的模倒數(shù)或乘法逆，記作

.

可以證明，如果

a與

m無公共素?cái)?shù)因子，則

a有唯一的模

m倒數(shù)(素?cái)?shù)是指除了1與自身外，不能被其他非零整數(shù)整除的正整數(shù))，反之亦然.

例如，

.

命題

元素屬于Zm的方陣

A模

m

可逆的充要條件是，m和det(A)沒有公共素?cái)?shù)因子，即m和det(A)互素.

顯然，所選加密矩陣必須符合該命題的條件.

a1357911151719212325a-11921153197231151725表12-2模26倒數(shù)表表12-2可用下列程序求得：m=26;fora=1:m

fori=1:m

ifmod(a*i,m)==1

fprintf('TheINVERSE(mod%d)ofnumber:%dis:%d\n',m,a,i)end;end;end若

則的模26逆矩陣：例12.3甲方收到與之有秘密通信往來的乙方的一個密文信息，密文內(nèi)容：WKVACPEAOCIXGWIZUROQWABALOHDKCEAFCLWWCVLEMIMCC按照甲方與乙方的約定，他們之間的密文通信采用Hill2密碼，密鑰為二階矩陣問這段密文的原文是什么？

解已知密鑰

的逆矩陣是

，根據(jù)解密方法，得到分組明文見表12-3.

于是，實(shí)際問題(甲)的解為：

GUDIANMIMASHIYIZIFUWEIJIBENJIAMIDANYUANDEMIMA

即為：“古典密碼是以字符為基本加密單元的密碼”.

12.3MATLAB求解C

=char(A)

將數(shù)組

A

轉(zhuǎn)換為字符數(shù)組

從32到127的整數(shù)對應(yīng)于可打印的ASCII字符.ASCII字符表中65到90的整數(shù)代表大寫字母A到Z（見表12-3）.序號分組密文密文表值明文表值分組明文1WK2311721GU2VA22149DI3CP316114AN4E513M

A19I5O1513M

C31A6I919S

X248H7G79I

W2325Y8I99I

Z00Z9U219I

R186F10O1521U

Q1723W11W235E

A19I表12-3分組明文12BA210J

19I13L122B

O155E14H814N

D410J15K119I

C31A16E513M

A19I17F64D

C31A18L1214N

W