第二次課程-綜合實驗

6合實，基礎性實驗后，本章將進入綜合實驗部分。在綜合實驗里通過編寫短程序或者調用工具箱來解決實驗問題，同時進一步熟悉 展等幾個部分。本章所編實驗問題均來自日常數學教學和科學研究，重點在于展現使，實驗問題：25開開NYNYYN輸出結果矩陣結圖6-1 fori=1:mforif

fori=1:tifif

20120131312315例1比較矩陣A

和

的對 6 6應位置，找出相同元素。

45 5

63 63解：分別調用程序一和程序二。res111223440516555res ifm==n&s==tforifif

ifres==[]開開NNYNY輸出刪除為'Z'的6-2functionfrep(x)fori=1:length(x)-1forj=i+1:length(x)%if

x11，2，2，3，3，-1，5，-1，4，5，2，2，去掉數列中的重復項，并輸解：x functionfrep(x)fori=1:length(x)-forj=i+1:length(x)ifx(j)==x(i)

xifdisp(['共刪除num2str(b)項重復項

x 6實驗練習中的函數sort可以實現數列排序功能本實驗是在不使用sort函數的情況下，開開輸入數輸入排序結冒泡排6-3NNNYNY變換x[j]與x[j+1]的j<length(x)-6-4functionpaixu(x)forforj=1:N- x(j)x(j+1)如果前項大于后項，交換兩項

110 產生一個m行n列、均值為M、標準差為s的隨機矩陣。 xx(j)x(j+1)改為x(j)WW6-5。functionpaixu(x)fori=1:Nforj=1:N- x(j)>

ifxNNNYYNY變換x[j]與x[j+1]的值j<length(x)-6-5要把這N個數從小到大排序。首先x[1]xN的范圍內，選出最小值與x[1]交換；然x[2xN范圍內，選出最小值與x[2交換；依次是x[3xN的范圍，這樣依實驗問題：x5x10的近似根（誤差<10-二分法的理論依據是根的存在性定理，定理的內容是：若函數fx在a,b上連續(xù)，fafb0fx在a,b內至少有一個實數根。二分法的基本思想是：首先選取區(qū)間a,bab，保留有根的半個區(qū)間aab或者abb

2

a

ab

a具體步驟（1）取區(qū)間a,b中2

f

0，

（2）fab0fabfb0，則

a bb；否aa2 2

ab1 ；得到新的有根區(qū)間a1,b1，將它看作區(qū)間a,b；（3）重復執(zhí)行（1a2x5x10的近似根（誤差<10-5）holdon;whileabs(a-b)>jingduiff(c)==0elseiff(c)*f(b)<0plot(c,0,'*'畫出每次二分區(qū)間的中點c%顯示每次二分區(qū)間的中點k=k+1 近方程的根，也就是曲線f

1x6-11234567890 6-6 法。設方程fx0在a,b內有一個單x*，且曲線為凹型，如圖6-7，則可用弦截法求根的近似值，計算步驟如下：取定一個端點b,fb，連接它與另一端點作弦交x軸x1b

bfbf

fyya obx6-7以x1,fx1為新端點，并連接端點b,fb作弦交x軸x2b

bfbfx1

fn1xnn

弦截法得到的迭代序列xn也是收斂的，且極限是方程fx0的根f

3x21fx0有根區(qū)間，并求根（精度106x2y22yex實驗問題：問題1.使用截痕法觀察馬鞍面zx2y2；問題2.繪制旋轉拋物面zx2y2問題1.對于空間曲面Fx,y,z0，通常可以采取截痕法來認識該曲面的特性，即用Fx,y,z0，通過研究截得交線的性質來認識曲面的性質。本實驗中，通過截痕法來研究馬鞍面zx2y2。問題2.通會在長方形區(qū)域內繪制函數圖像這可能導致某些函數圖像失真針這個問題，可以利用令來精細繪制空間特殊曲面的圖像。此處將通zx2y2的圖像，來了解怎樣繪制出理想的圖像。holdonplot3(px,py',pz','r*')

6-8zx2y2x=26-8x=2的交線，右圖是單獨畫出了交線。2，首先輸入如下程序：

6-9zx2y2的圖像（一圖6-9所示的旋轉拋物面在xOy坐標平面上的投影是一個正方形區(qū)域，而不是圓域。利用提供的find命令可以對其繪圖過程進行精選處理，也就是將橫縱坐標的平4的那些點去掉。

6-10zx2y2的圖像（二zx2y2y=3，z=222a2

y

zc2

1x2 y z精細繪 1的圖像a2 c2重積分計算是多元函數微分學中較為復雜的計算。本節(jié)將講解利用求解二重、三重積分的方法。求解分為計算數值解和精確解，計算數值解時，可以使用件的內置函數；在精確計算時，則通常首先利用繪圖，將重積分轉化為累次積矩形區(qū)域 x,yaxb,cyd上的二重積分1Ixy2dxdyDxy0x11y3D解：dblquadf=inline('x.*I=dblquad(Iinlinem文件定義函數的情況下，直接定義函數的表達式；'x.*長方體閉區(qū)域xyaxbcyd,ezftriplequad2Icosxzdxdydz，其中x解：f=@(x,y,z)cos(x)+I=triplequad(f,0,pi,0,2,-I

0很多重積分是轉化為累次積分進行計算的。中計算累次積分可以調用int函數y3Iy3I2dyxydx 解：symsxI=4Iarctanydxdy，其中Dx2y216x2y24，及直線yx y0x=-5:0.01:5;y=x;holdontext(2,2.6,'y=x');a=0:0.01:2*x1=4*cos(a);y1=4*sin(a);plot(x1,y1,'-r')x2=2*cos(a);y2=2*sin(a);text(sqrt(15),1,'x^2+y^2=16');y3=0;x3=0;plot(x,y3,'-k');plot(x3,x,'k');text(sqrt(3),1,'x^2+y^2=4');text(3,-0.2,'y=0');text(3,1.5,'D');hold543y=2D x2y2= y=

x2 6-11D6-11，將二重積分轉化為：IarctancosddDDD,024int symsrtheta;f=r*theta;Ix2例5利用柱面坐標計算Izex2y2dxdydz，其中為z 及x2（h0）解：利 繪制積分區(qū)域的程序如下x=-2:0.08:2;y=-[X,Y]=meshgrid(x,y)Z1=sqrt(X.^2+Y.^2);Z2=2.*ones(size(X));mesh(X,Y,Z1);text(1,2,2.8,'z=sqrt(x^2+y^2)');holdontext(2,-hiddenoff

6-12xoys=abs(Z1-xx=s.*X;yy=s.*Y;zzplot3(xx(s~=0),yy(s~=0),zz(s~=0),'b*');t=zeros(size(yy(s~=0)));plot3(xx(s~=0),yy(s6-13

6-13zD,020zz有zex2y2dxdydzze2dd

2dd，計算該三重積分令為 symsrhzf=z*r*exp(r*r);I1=int(f,r,0,z);

0DI2=int(I1,theta,0,2*pi);I=int(I2,z,0,h)I-(pi*(h^2-exp(h^2)+sinydxdyDyxxy2 x2y2z24a2x2y22ax0（a0）所截且含在圓柱面內的立zdxdydz，其中體Vx2y2z24（z0）Vx2y23z 元線性回歸的實現。 N N

。E(Y)01x1pxp01,p假定對這一組變量(

p;Y作了n T01,pT11

(T

x1p

0

y1

y其中X

2p,

1y2，XTX

1x1x

xnp

p

yny

y?0?x?x1 p通常2

?

)2ninni

1

p當n較小時，通常取2

?

)2np1

1

pY

21

間是否存性關系也需要進行檢驗。即要檢驗假設H012p0SRSQSRSQ(np 且當H0成立

(yy)2，稱其為總離差平方和， (yy?)2，稱其為剩余平方和， (y?y)2，稱其為回歸平方和。Fnp1

~Fpnp1 Fpnp1)，并根據樣本求出F SR

的觀測值，若FF(p,np1)，H，認為Y與自變0SQ(np0

21 之間 性相關關系。反之，則接受H0，認為回歸效果不顯著命令為：[bbintrrint其中，YyY2X

yny

x1pX

x2p x np間，stats表示用于檢驗回歸模型的統(tǒng)計量。

6-2所示。6-2123456789y9y7關系。程序為：y=[33.240.338.746.841.437.539.040.730.152.938.231.843.344.142.533.34.248.035.940.436.845.2x1=[3.55.35.15.84.26.06.85.53.17.24.54.98.06.56.63.76.27.04.04.55.95.64.8x2=[9201833 11233539217403523332734x3=[6.16.47.46.77.55.96.04.05.88.35.06.47.67.05.04.45.57.06.03.54.94.38.0subplot(1,3,1plot(x1,y'*title('yx1的散點圖subplot(1,3,2),plot(x2,y'*title('yx2的散點圖'),subplot(1,3,3),plot(x3,y,'*'),title('yx3的散點圖')

6-146-14可以看出這些點大致分布在一條直線附近，故因變量與各自變量之間有著良好的y01x12x2Y=[33.240.338.746.841.437.539.040.730.152.938.231.843.344.142.533.34.248.035.940.436.845.2[bbintrrintstats]=regress(Y,X,alpha);bbintstats ResidualCaseOrder6420

Case6-150.05；統(tǒng)計FFpR2越接近1FF1(1n H0，F值越大，說明回歸方程越顯著；與F對應的概率p

rcoplt(,rint)新進行回歸。本例題中結果：統(tǒng)計量stats顯示：R20.9126，F69.5702 0p0.0000.05y6-15顯示大部分數據的殘差都離零點比較近，故可得回歸方程

2其中，beta表示估計出的回歸系數；r表示殘差J表示Jacobian矩陣xy表示輸入數據x為矩陣（對于一元非線性回歸為列向量，y為n為列向量；model表示事先用m-文件定義的非線性函數；beta0表示回歸系數的初值。alpha為顯著性水平。和誤差估計命令為：其中nlinfit或nintool所得的回歸函數在x處的值Y及值的顯著性為1-alpha的置信Y±DELTA。2（1）functionyhat=volum(beta,x)y=[6.428.209.589.59.7109.939.9910.4910.5910.6010.8010.6010.90beta0=[81即獲得回歸模型：y 6-169876 圖6-161.1（x1（x2yaxb實驗問題：無理數、數的近似計算，以值計算為例以和E為代表的數在自然科學中扮演著重要的角色很長一段時間里對于熟，本實驗將以π值計算為例，介紹數、無理數的近似計算。圓周率π的值精確到任意精度的方法。他提出，圓周長大于內接正多邊形周長而小于外9622322 ，2000多年前的《周髀算經》中稱“周三徑一”π的第一個近似值。魏晉24576邊形，而得到一個結論：3.1415926＜81593π2

2 2

2 2 27之后多種表達式陸續(xù)給出。1706JohnMachin7arctanxx

計算到了圓周率的100位。1914年，數學家SrinivasaRamanujan了下面的公式 2

n044n

在1985年，Gosper用這個公式計算到了圓周率的17500000位。，David和GregoryChudnovsky了下面的公式：1

并在1994年計算到 . 426880

(6n)!(545140134n13 (n!)3(3n)!(64032（BBP公式： 1

16n8n n1位數。5萬億位。π值進行近似計算，分別為：劉徽割圓術，冪級數計算方法，n邊形的邊長為an，圓內接正nSn。根據勾股定理，邊

an12 422 42 42SSnS2nS2S2nSna(1)=sqrt(3);b(1)=3*sqrt(3)/2;fori=2:6a(i)=sqrt(2-sqrt(4-a(i-1)^2));result利用與arctanx

2n

.x1

.4arctan14(111 n=1000;選擇展開式的次數digits(22);定義計算過程中的精度fork=1:nvpa(s,20201 dx

0 將積分區(qū)間[0,1]分成n等份，在每一個小區(qū)間上，選取中點i，可得近似計算程序為n=50;定義等分積分區(qū)間數fork=1:length(i)-1蒙特卡洛(MonteCarlo)模擬，也稱統(tǒng)計模擬方法，是二十世紀四十年代中期由于科學技使用蒙特卡洛方法計算6-17nmm、n4 分大時 6-17n=500隨機取點數fori=1:nifa(1)^2+a(2)^2<=1pi=pi+1;4*蒙特卡洛方法簡單易行，但是該算法收斂的速度比較慢。比較表6-3與表6-4的結果可E的近似值，計算公式為：E11111forfori=1:k

利用韋達公式給出依據下述四個公式：2 ，2

(1)n1，3

，(

n1(2n

n1(2n8

(2n

π 基于關系 dx ，利用蒙特卡洛方法近似計算值（提示：隨機地向正方01x MN充分大時，MN就近似等22

轉化為微分方程，從而建立抽象系統(tǒng)發(fā)展變化的動態(tài)模型，簡記為GM(hn)。該模型中hax

式（6.9.1）dxx及控制量udt(0)(2),,x(0)

，對X

k(1(2),x(1nx(1k)x(0i)k1,2,ndt

ax(1)

其中，a,u為待定參數，a稱作發(fā)展系數，反映的發(fā)展態(tài)勢；u稱作灰作用量，反映數 (1)(x(1)(k1))aZ(1)(k1)

(1x(1k1x(1k1x(1kx(0k1X(1在(k1序列；

(t)

12

(k)

(k1)]

dt

在(k

x(0)(k1)a[1(x(1)(k)x(1)(k1))]2

1(x(1)(1)x(1) 1x(0) x(0)(3) (x(1)(2)x(1) 1

u

1(x(1)(n1)

(n))1(x(1)(1)x(1) x(0) x(0)(3)

(x(1)(2)x(1) 令Y

，B

，[a,

1(x(1)(n1)x(1) （6.9.5）

Y

其中

(

還原到原始數據得

?(0)?(1)?(1)（6.9.8（6.9.9）

GM(1,1aaa2的條件下，這時建立的GM(1,1)模型才有意義。一般來說有如下結論：當0.3a0.5時，GM(1,1)模型可用于短期，中長期慎用當0.5a0.8時，用GM(1,1)模型作短期應十分謹慎當0.8a1.0時，應采用殘差修正GM(1,1驗法。一般地將GM(1,1模型的精度檢驗等級分為四級：相對誤差檢驗法、后驗差檢驗法6-5。

cp1級（好c0.952級（合格0.80p3級（勉強0.70p4級（不合格p1某公司組織促銷活動，1號~86-61號~5號的數據建立灰色模型，對6號至8號的銷售量進行，并與實際值進行對比。表6-6銷售數 63748解：x00=[56985703570757195724572957385751x0=x00-5696*ones(1,8);%提高精forx1對原始數據進行累加fori=1:4forfori=1:4BBu=c(2,1)a,ufort=1:8x11;建立時間響應函數for *ones(1,8)%forq1(i)=x00(i)-q11(i)=abs((x00(i)-x13(i))/x00(i));x1 Baux131.0e+003 q1 q11= ?(1)(k1)(x(0)(1)u)eaku15.91e042k 其中0.3a0.420.5，故所建立的GM(1,1)模型只適用于短期表6-76號~8號銷售量及誤差檢序 實際數據 數據?(0)殘差(k相對誤差23456781.某地區(qū)十年的農業(yè)產值如表6-8，用灰色模型未來5年的農業(yè)產值表6-810年的農業(yè)產 年 產 的加密、（plaintext（ciphertext簡單的加密過程模型加密發(fā)送者明 密加密器道器接收 明圖6-18加密過程模

、干在學中，凱撒是一種最為古老的加密體制。在古羅馬的時候，凱撒曾用此方法3AD，BE，以此類推?？傻妹魑拿魑淖帜副恚好芪淖帜副恚合枺℉illCipher）是基于矩性變換的替換，由LesterS.Hill在1929年發(fā)明。希爾相對于凱撒的好處是很好地隱藏了字母的頻率信息，使得傳統(tǒng)的通過字母頻率來的方式失效本實驗以Hill2為例講碼生成和的程序同。一般地，Hilln，明文每n個分為一組，按組轉換成密文。Hill2的加密與流程如下260-251-1ABCDE IJKLM12345 9NOPQR VWXYZ 0

A 0 將明文字母依次按每兩個字母一組依據表6-9aibibiAaimod26bi

aA1biiii關于模運算mod模ma，babkm，kZa模m等價于bab abmodm 模m逆：設aZm，若存在bZm使得ab1modm，稱a有模mba1modm命題：整數a有模m逆元a與m無