Matlab實(shí)驗(yàn)綜合實(shí)驗(yàn)內(nèi)容提要π的計(jì)算與數(shù)值積分Hill密碼飲酒駕車模型圓周率π的計(jì)算歷程所謂“圓周率”是指一個(gè)圓的周長(zhǎng)與其直徑的比值。古今中外，許多人致力于圓周率的研究與計(jì)算。為了計(jì)算出圓周率的越來(lái)越好的近似值，一代代的數(shù)學(xué)家為這個(gè)神秘的數(shù)貢獻(xiàn)了無(wú)數(shù)的時(shí)間與心血。回顧歷史，人類對(duì)π的認(rèn)識(shí)過(guò)程，反映了數(shù)學(xué)和計(jì)算技術(shù)發(fā)展情形的一個(gè)側(cè)面。π的研究，在一定程度上反映這個(gè)地區(qū)或時(shí)代的數(shù)學(xué)水平。德國(guó)數(shù)學(xué)家康托說(shuō)：“歷史上一個(gè)國(guó)家所算得的圓周率的準(zhǔn)確程度，可以作為衡量這個(gè)國(guó)家當(dāng)時(shí)數(shù)學(xué)發(fā)展水平的指標(biāo)。”直到19世紀(jì)初，求圓周率的值應(yīng)該說(shuō)是數(shù)學(xué)中的頭號(hào)難題。為求得圓周率的值，人類走過(guò)了漫長(zhǎng)而曲折的道路。實(shí)驗(yàn)時(shí)期基于對(duì)一個(gè)圓的周長(zhǎng)和直徑的實(shí)際測(cè)量而得出的。在古代世界，實(shí)際上長(zhǎng)期使用π＝3這個(gè)數(shù)值。最早見(jiàn)于文字記載的有基督教《圣經(jīng)》中的章節(jié)，其上取圓周率為3。這一段描述的事大約發(fā)生在公元前950年前后。幾何法時(shí)期阿基米德：科學(xué)地研究這一常數(shù)的第一個(gè)人；是他首先提出了一種能夠借助數(shù)學(xué)過(guò)程而不是通過(guò)測(cè)量，能夠把π的值精確到任意精度的方法。由此，開(kāi)創(chuàng)了圓周率計(jì)算的第二階段。圓周長(zhǎng)大于內(nèi)接正多邊形周長(zhǎng)而小于外切正多邊形周長(zhǎng)．據(jù)說(shuō)阿基米德用到了正96邊形才算出他的值域。

在中國(guó)劉徽：公元263年前后，劉徽提出著名的“割圓術(shù)”求出了比較精確的圓周率。他發(fā)現(xiàn)：當(dāng)圓內(nèi)接正多邊形的邊數(shù)不斷增加后，多邊形的周長(zhǎng)會(huì)越來(lái)越逼近圓周長(zhǎng)，而多邊形的面積也會(huì)越來(lái)越逼近圓面積。于是，劉徽利用正多邊形面積和圓面積之間的關(guān)系，從正六邊形開(kāi)始，逐步把邊數(shù)加倍：正十二邊形、正二十四邊形，正四十八邊形……，一直到正三○七二邊形，算出圓周率等于三點(diǎn)一四一六，將圓周率的精度提高到小數(shù)點(diǎn)后第四位。在中國(guó)祖沖之:在劉徽研究的基礎(chǔ)上，進(jìn)一步地發(fā)展，經(jīng)過(guò)既漫長(zhǎng)又煩瑣的計(jì)算，一直算到圓內(nèi)接正24576邊形，而得到一個(gè)結(jié)論：

3.1415926＜π＜3.1415927

同時(shí)得到π的兩個(gè)近似分?jǐn)?shù)：約率為22／7；密率為355／113。他算出的π的8位可靠數(shù)字，不但在當(dāng)時(shí)是最精密的圓周率，而且保持世界記錄九百多年。以致于有數(shù)學(xué)史家提議將這一結(jié)果命名為“祖率”。分析法時(shí)期這一時(shí)期人們開(kāi)始擺脫求多邊形周長(zhǎng)的繁難計(jì)算，利用無(wú)窮級(jí)數(shù)或無(wú)窮連乘積來(lái)算π。1593年，韋達(dá)給出

1）π的最早分析表達(dá)式；2）公式優(yōu)美，令人贊嘆；3）借助數(shù)字2，通過(guò)一系列的加、乘、除和開(kāi)平方就可算出π值。接著有多種表達(dá)式出現(xiàn)。如沃利斯1650年給出：發(fā)現(xiàn)了下面的公式1706年，英國(guó)天文學(xué)教授JohnMachin利用并利用這個(gè)公式計(jì)算到了圓周率的100位.1914年，印度數(shù)學(xué)家SrinivasaRamanujan發(fā)表了下面的公式：在1985年，Gosper用這個(gè)公式計(jì)算到了圓周率的17500000位．

1989年，David和GregoryChudnovsky發(fā)表了下面的公式

并在1994年計(jì)算到了4044000000位．它的另一種形式是1995年，由DavidBailey,PeterBorwein和SimonPlouffe共同發(fā)表了下面的圓周率計(jì)算公式（簡(jiǎn)稱BBP公式）

該公式的最大優(yōu)點(diǎn)在于：經(jīng)后來(lái)人將該公式變形后打破了傳統(tǒng)的計(jì)算方法，可以直接計(jì)算圓周率的任意第n位數(shù)，而不是先計(jì)算前面的n-1位數(shù)。1997年，F(xiàn)abriceBellard發(fā)表了一個(gè)比BBP算法更快的公式從而，大大降低了圓周率近似值的計(jì)算量.電子計(jì)算機(jī)的出現(xiàn)使π值計(jì)算有了突飛猛進(jìn)的發(fā)展。1949年美國(guó)馬里蘭州阿伯丁的軍隊(duì)彈道研究實(shí)驗(yàn)室首次用計(jì)算機(jī)（ENIAC）計(jì)算π值，一下子就算到2037位小數(shù)，突破了千位數(shù)。1989年美國(guó)哥倫比亞大學(xué)研究人員用克雷－2型和IBM－VF型巨型電子計(jì)算機(jī)計(jì)算出π值小數(shù)點(diǎn)后4.8億位數(shù)，后又繼續(xù)算到小數(shù)點(diǎn)后10.1億位數(shù)，創(chuàng)下最新的紀(jì)錄。

2010年1月7日——法國(guó)一工程師將圓周率算到小數(shù)點(diǎn)后27000億位。

2010年8月30日——日本計(jì)算機(jī)奇才近藤茂利用家用計(jì)算機(jī)和云計(jì)算相結(jié)合，計(jì)算出圓周率到小數(shù)點(diǎn)后5萬(wàn)億位。π的計(jì)算方法古典方法（略）級(jí)數(shù)逼近法數(shù)值積分法隨機(jī)模擬法π的冪級(jí)數(shù)計(jì)算方法思想：利用一些特殊函數(shù)的冪級(jí)數(shù)展開(kāi)式計(jì)算圓周率.例如編寫(xiě)下面的程序:n=10;%選擇展開(kāi)式的次數(shù)s=0;digits(22);%定義計(jì)算過(guò)程中的精度f(wàn)ork=1:ns=s+4*(-1)^(k+1)/(2*k-1);endvpa(s,20)%定義顯示精度為20位上機(jī)實(shí)驗(yàn):

請(qǐng)按照上述公式分別編程計(jì)算，并與就精度與疊加次數(shù)進(jìn)行比較，能得出怎樣的結(jié)論？函數(shù)的泰勒展開(kāi)式的matlab指令

MATLAB軟件求一元函數(shù)泰勒展式的命令是taylor,具體使用格式如下：taylor(f)求函數(shù)f的5階Maclaurin展開(kāi)式；taylor(f,n)求函數(shù)f的n-1階Maclaurin展開(kāi)式；taylor(f,n,a)求函數(shù)f在x=a處的n-1階泰勒展開(kāi)式；

taylor(f,a)求函數(shù)f在x=a處的5階泰勒展開(kāi)式；對(duì)Taylor級(jí)數(shù)，|x|越小，級(jí)數(shù)收斂速度越快使用這種算法收斂速度要快的多。Machin公式英國(guó)天文學(xué)家麥琴（Machin）于1706年給出，他利用這個(gè)公式求得了π的第100位小數(shù)。利用類似方法，到1948年英國(guó)的弗格森和美國(guó)的倫奇共同發(fā)表了π的第808位小數(shù)，成為人工計(jì)算圓周率值的最高紀(jì)錄。π的數(shù)值積分計(jì)算方法根據(jù)上述微積分公式，將求解圓周率轉(zhuǎn)換為數(shù)學(xué)上經(jīng)常遇到的一類問(wèn)題——數(shù)值積分問(wèn)題。

實(shí)際中常遇到這樣的數(shù)值積分，被積函數(shù)不是由解析式給出，或者雖然是解析式但求其原函數(shù)十分困難。這時(shí)就要考慮定積分的近似計(jì)算了。定積分的近似計(jì)算梯形法

分割區(qū)間，用梯形面積代替曲邊梯形面積拋物線法（Simpson法）

分割區(qū)間，在每個(gè)區(qū)間添加其中點(diǎn)，用過(guò)兩端點(diǎn)及中點(diǎn)的拋物線面積代替曲邊梯形面積。n=50;%定義等分積分區(qū)間數(shù)，可以更改i=0:1/n:1;s=0;fork=1:length(i)-1s=s+(1/(1+((i(k)+i(k+1))/2)^2))*1/n;end4*s梯形法等分區(qū)間數(shù)n103.14242598500110203.14180098689309503.141625986923001003.141600986923125003.1415929869231210003.1415927369231350003.14159265692313下表給出的是不同等分區(qū)間數(shù)計(jì)算出的

圓周率π的近似值拋物線法Matlab提供了積分函數(shù)quad()clear;f=inline('4./(1+x.*x)');a=0;b=1;n=1;z=quad(f,a,b);fprintf('z=%.10f\n',z)按照實(shí)際問(wèn)題所遵循的概率統(tǒng)計(jì)規(guī)律，用電子計(jì)算機(jī)進(jìn)行直接的抽樣試驗(yàn)，然后計(jì)算其統(tǒng)計(jì)參數(shù)。蒙特卡洛方法也可以人為地構(gòu)造出一個(gè)合適的概率模型，依照該模型進(jìn)行大量的統(tǒng)計(jì)實(shí)驗(yàn)，使它的某些統(tǒng)計(jì)參量正好是待求問(wèn)題的解。在物理學(xué)的各種各樣問(wèn)題中得到廣泛的應(yīng)用。1）蒙特卡洛方法圓周率π的隨機(jī)模擬計(jì)算方法cs=0n=500%隨機(jī)取點(diǎn)數(shù)fori=1:na=rand(1,2);ifa(1)^2+a(2)^2<=1cs=cs+1endend4*cs/n１１０正方形中任取一點(diǎn)，落入四分之一圓中的概率為正方形中任取n個(gè)點(diǎn)，其中m個(gè)點(diǎn)落入四分之一圓中，則當(dāng)m、n充分大時(shí)，π的蒙特卡洛計(jì)算方法

算法優(yōu)點(diǎn)：簡(jiǎn)單易行；缺點(diǎn)：算法收斂的速度比較慢．

在精度要求不是很高的情況下，這種取隨機(jī)數(shù)進(jìn)行數(shù)據(jù)模擬的方法還是有一定的實(shí)用價(jià)值，通過(guò)簡(jiǎn)單編程我們可以模擬出許多數(shù)學(xué)現(xiàn)象．

例:計(jì)算定積分1）隨機(jī)地向正方形內(nèi)投點(diǎn)，最后統(tǒng)計(jì)落在曲線下的點(diǎn)數(shù)M；2）當(dāng)總的擲點(diǎn)數(shù)N充分大時(shí)，N/M就近似等于積分值I。上機(jī)實(shí)驗(yàn)：基于關(guān)系式利用蒙特卡洛方法近似計(jì)算保密通訊無(wú)論在軍事、政治、經(jīng)濟(jì)還是日常生活中都起著非常重要的作用。為了將信息傳遞給己方的接受者，同時(shí)又要防止他人（特別是敵人）知道信息的內(nèi)容，必須將要傳遞的信息（明文）加密，變成密文后發(fā)送出去，這樣，即使敵方得到密文也看不懂，而己方的接受者收到密文后卻可以按照預(yù)先定好的方法加以解密。Hill密碼密碼可分為古典密碼和現(xiàn)代密碼本實(shí)驗(yàn)主要介紹古典密碼的加密與破譯原理，同時(shí)介紹如何用Matlab編程來(lái)實(shí)現(xiàn)加密、解密和破譯過(guò)程。問(wèn)題背景和實(shí)驗(yàn)?zāi)康?/p>

古典密碼：以字符為基本加密單元；

現(xiàn)代密碼：以信息塊為基本加密單元。加密信息傳遞過(guò)程明文（信息）加密器密文密文明文（信息）解密器普通信道發(fā)送敵方截獲破譯發(fā)送方接收方Hill2密碼的加密過(guò)程Hill2密碼中所用的數(shù)學(xué)手段是矩陣運(yùn)算。

加密過(guò)程：①

將漢語(yǔ)拼音的26個(gè)字母

與0到25

之間的整數(shù)建立一一對(duì)應(yīng)關(guān)系，稱為字母的表值，然后根據(jù)明文字母的表值，將明文信息用數(shù)字表示。ABCDEFGHIJKLM12345678910111213NOPQRSTUVWXYZ1415161718192021222324250設(shè)通訊雙方所給出的26個(gè)字母的表值如下：注：這里假定明文中只使用26個(gè)大寫(xiě)字母Hill2密碼的加密過(guò)程②

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

將明文字母分組。Hill2

使用的是二階矩陣，所以將明文字母每2個(gè)一組（可以推廣至Hilln密碼）。查出每個(gè)字母的表值，這樣，每組字母構(gòu)成一個(gè)二維列向量。若最后僅剩一個(gè)字母，則補(bǔ)充一個(gè)沒(méi)有實(shí)際意義的啞字母（啞元）。這樣使得每組都有2個(gè)字母，④令=

A

，由的兩個(gè)分量反查字母表值表，得到相應(yīng)的兩個(gè)字母，即為密文字母。Hill2加密舉例例：設(shè)明文為“KMXYSXX”（昆明學(xué)院數(shù)學(xué)系），試給出這段明文的Hill2

密文。其中加密矩陣為

將明文字母分組：

KMXYSXXX最后的一個(gè)字母X為啞字母，無(wú)實(shí)際意義。解:ABCDEFGHIJKLM12345678910111213NOPQRSTUVWXYZ1415161718192021222324250查表得每組字母的表值，得到4個(gè)二維列向量：將上述4個(gè)二維向量左乘密鑰矩陣A

得：作模26運(yùn)算，將所有的數(shù)都化為0到25之間的整數(shù)：Hill2

加密舉例反查字母表值得每個(gè)向量對(duì)應(yīng)的字母組為：KMXYSXXKMVWOTTTHill2加密Hill2

加密舉例ABCDEFGHIJKLM12345678910111213NOPQRSTUVWXYZ1415161718192021222324250KMVWOTTT問(wèn)題：怎樣解密？明文字母查表值分組一組向量加密矩陣左乘一組新的向量反查表值密文Hill2

加密過(guò)程模運(yùn)算先查出密文字母“KMVWOTTT”所對(duì)應(yīng)的向量：在模運(yùn)算下解方程組：

A

=

解密：加密的逆過(guò)程，將加密過(guò)程逆轉(zhuǎn)回去即可。Hill2

解密過(guò)程上面的向量是由經(jīng)過(guò)模26運(yùn)算得來(lái)的，現(xiàn)在的問(wèn)題是怎樣逆轉(zhuǎn)回去？例：怎么得到密文“KMVWOTTT

”的原文模m可逆記定義1：設(shè)A為定義在集合Zm上的n階方陣，若存在一個(gè)定義在Zm上的方陣B，使得

則稱A模m可逆，B為A的模m逆矩陣，記為定義2：設(shè)a

Zm，若存在b

Zm使得ab=1(modm)，則稱b

為

a

的模m倒數(shù)或乘法逆，記作b=a-1(modm)。注:a,b都是

Zm中的數(shù)命題：定義在集合Zm

上的n

階方陣A模m

可逆的充要條件是：m和

det(A)

無(wú)公共素?cái)?shù)因子，即m

與det(A)

互素。Hill2密碼的加密矩陣必須滿足上述條件。m=26m的素?cái)?shù)因子只有2和13定義在Z26上的方陣A模26可逆的充要條件是：模m可逆det(A)不能被2和13整除

問(wèn)題：是否

Zm中所有的數(shù)都存在模m倒數(shù)？a

存在唯一的模m倒數(shù)a

與m無(wú)公共素?cái)?shù)因子Z26中具有模26

倒數(shù)的整數(shù)及其模26倒數(shù)表a1357911151719212325a-11921153197231151725m=26時(shí)的Matlab程序見(jiàn)教材第121頁(yè)fuluA.m

可以用來(lái)判斷一個(gè)2階矩陣在模26運(yùn)算下是否可逆，并在可逆的前提下給出其模26逆矩陣。模26可逆

思考：如何用Matlab編程來(lái)找出所有模m

倒數(shù)的整數(shù)及其模m倒數(shù)？在模運(yùn)算下解方程組：A

=

Hill2

解密過(guò)程？

問(wèn)題：如何計(jì)算

？模m逆矩陣的計(jì)算設(shè)B=kA*為A

的模26逆，其中k

為待定系數(shù)A*為A

的伴隨矩陣本計(jì)算方法可推廣到求矩陣A

的模m逆矩陣Hill2

解密過(guò)程設(shè)加密矩陣用B左乘密文對(duì)應(yīng)的向量得：模26運(yùn)算后得：查表后得明文分別為：KMXYSXXXMatlab的Hill2加密程序見(jiàn)附錄B，相應(yīng)解密程序見(jiàn)附錄C。？Hill2加密過(guò)程總結(jié)①

通訊雙方確定加密矩陣

(密鑰)和字母的表值對(duì)應(yīng)表②

將明文字母分組，通過(guò)查表列出每組字母對(duì)應(yīng)的向量

③

令=

Amod(m)

，由

的分量反查字母表值表，

得到相應(yīng)的密文字母若明文只含奇數(shù)個(gè)字母，則補(bǔ)充一個(gè)啞元Hill2解密過(guò)程總結(jié)①將密文字母分組，通過(guò)查表列出每組字母對(duì)應(yīng)的向量②

求出加密矩陣A

的模m逆矩陣

B③

令=

B

mod(m)

，由

的分量反查字母表值表，

得到相應(yīng)的明文字母甲方收到乙方（己方）的一個(gè)密文信息，內(nèi)容為：Hill2解密舉例WKVACPEAOCIXGWIZUROQWABALOHDKCEAFCLWWCVLEMIMCC按照甲方與乙方的約定，他們之間采用Hill2密碼，密鑰為，字母表值見(jiàn)下表，問(wèn)這段密文的原文是什么？ABCDEFGHIJKLM12345678910111213NOPQRSTUVWXYZ1415161718192021222324250Hill2解密舉例①將密文字母分組，通過(guò)查表列出每組字母對(duì)應(yīng)的向量②

求出加密矩陣A

的模26逆矩陣③

用B

左乘每組密文字母組成的向量，然后再反查字母表值表，得到相應(yīng)的明文字母序號(hào)分組密文密文表值明文表值分組明文1WK2311721GU2VA22149DI3CP316114AN4EA51139MI5OC153131MA6IX924198SH序號(hào)分組密文密文表值明文表值分組明文7GW723925IY8IZ9090IZ9UR211896IF10OQ15172123UW11WA23159EI12BA21109JIHill2解密舉例序號(hào)分組密文密文表值明文表值分組明文13LO121525BE14HD841410NJ15KC11391IA16EA51139MI17FC6341DA18LW12231425NY序號(hào)分組密文密文表值明文表值分組明文19WC233211UA20VL2212144ND21EM513513EM22IM913913IM23CC3311AAHill2解密舉例即：“古典密碼是以字符為基本加密單元的密碼”GUDIANMIMASHIYIZIFUWEIJIBENJIAMIDANYUANDEMIMAAWKVACPEAOCIXGWIZUROQWABALOHDKCEAFCLWWCVLEMIMCC原文Hill2解密舉例密文ABCDEFGHIJKLM12345678910111213NOPQRSTUVWXYZ1415161718192021222324250經(jīng)分析該密文是用Hill2密碼加密，且密文(U,C)和(R,S)分別對(duì)應(yīng)明文(T,A)和(C,O)，問(wèn)能否破譯這段密文？Hill2密碼破譯舉例MOFAXJEABAUCRSXJLUYHQATCZHWBCSCP我方截獲一段密文猜測(cè)密文是由26個(gè)字母組成，即m=26，

經(jīng)破譯部門通過(guò)大量的統(tǒng)計(jì)分析和語(yǔ)言分析確定表值破譯這段密文的關(guān)鍵是找到“密鑰”和字母對(duì)應(yīng)的表值密文(U,C)和(R,S)分別對(duì)應(yīng)明文(T,A)和(C,O)Hill2密碼破譯舉例查字母表值PCP、C

模26可逆可唯一確定加密矩陣AHill2密碼破譯舉例注：這里的運(yùn)算都是在模運(yùn)算意義下進(jìn)行HEWILLVISITACOLLEGETHISAFTERNOON相應(yīng)的Matlab程序見(jiàn)附錄D

得到加密矩陣的模26逆矩陣后，根據(jù)前面的解密方法即可得密文的原文Hill2密碼破譯舉例飲酒駕車模型在2003年全國(guó)道路交通事故死亡人數(shù)中，飲酒駕車造成的占有相當(dāng)比例，為此，國(guó)家發(fā)布了新的《車輛駕駛?cè)藛T血液、呼氣酒精含量閾值與檢驗(yàn)》國(guó)家標(biāo)準(zhǔn)。在新標(biāo)準(zhǔn)下，大李在中午12點(diǎn)喝了一瓶啤酒，下午6點(diǎn)檢查時(shí)符合標(biāo)準(zhǔn)，接著晚上又喝了一瓶，但凌晨2點(diǎn)檢查時(shí)卻被定為飲酒駕車，為什么喝同樣多的酒，兩次檢查結(jié)果不一樣？藥物在體內(nèi)的分布與排除

藥物進(jìn)入機(jī)體形成血藥濃度(單位體積血液的藥物量)

血藥濃度需保持在一定范圍內(nèi)——給藥方案設(shè)計(jì)

藥物在體內(nèi)吸收、分布和排除過(guò)程——藥物動(dòng)力學(xué)

建立房室模型——藥物動(dòng)力學(xué)的基本步驟

房室——