Matlab实验综合实验

Matlab实验综合实验内容提要π的计算与数值积分Hill密码饮酒驾车模型圆周率π的计算历程所谓“圆周率”是指一个圆的周长与其直径的比值。古今中外，许多人致力于圆周率的研究与计算。为了计算出圆周率的越来越好的近似值，一代代的数学家为这个神秘的数贡献了无数的时间与心血。回顾历史，人类对π的认识过程，反映了数学和计算技术发展情形的一个侧面。π的研究，在一定程度上反映这个地区或时代的数学水平。德国数学家康托说：“历史上一个国家所算得的圆周率的准确程度，可以作为衡量这个国家当时数学发展水平的指标。”直到19世纪初，求圆周率的值应该说是数学中的头号难题。为求得圆周率的值，人类走过了漫长而曲折的道路。实验时期基于对一个圆的周长和直径的实际测量而得出的。在古代世界，实际上长期使用π＝3这个数值。最早见于文字记载的有基督教《圣经》中的章节，其上取圆周率为3。这一段描述的事大约发生在公元前950年前后。几何法时期阿基米德：科学地研究这一常数的第一个人；是他首先提出了一种能够借助数学过程而不是通过测量，能够把π的值精确到任意精度的方法。由此，开创了圆周率计算的第二阶段。圆周长大于内接正多边形周长而小于外切正多边形周长．据说阿基米德用到了正96边形才算出他的值域。

在中国刘徽：公元263年前后，刘徽提出著名的“割圆术”求出了比较精确的圆周率。他发现：当圆内接正多边形的边数不断增加后，多边形的周长会越来越逼近圆周长，而多边形的面积也会越来越逼近圆面积。于是，刘徽利用正多边形面积和圆面积之间的关系，从正六边形开始，逐步把边数加倍：正十二边形、正二十四边形，正四十八边形……，一直到正三○七二边形，算出圆周率等于三点一四一六，将圆周率的精度提高到小数点后第四位。在中国祖冲之:在刘徽研究的基础上，进一步地发展，经过既漫长又烦琐的计算，一直算到圆内接正24576边形，而得到一个结论：

3.1415926＜π＜3.1415927

同时得到π的两个近似分数：约率为22／7；密率为355／113。他算出的π的8位可靠数字，不但在当时是最精密的圆周率，而且保持世界记录九百多年。以致于有数学史家提议将这一结果命名为“祖率”。分析法时期这一时期人们开始摆脱求多边形周长的繁难计算，利用无穷级数或无穷连乘积来算π。1593年，韦达给出

1）π的最早分析表达式；2）公式优美，令人赞叹；3）借助数字2，通过一系列的加、乘、除和开平方就可算出π值。接着有多种表达式出现。如沃利斯1650年给出：发现了下面的公式1706年，英国天文学教授JohnMachin利用并利用这个公式计算到了圆周率的100位.1914年，印度数学家SrinivasaRamanujan发表了下面的公式：在1985年，Gosper用这个公式计算到了圆周率的17500000位．

1989年，David和GregoryChudnovsky发表了下面的公式

并在1994年计算到了4044000000位．它的另一种形式是1995年，由DavidBailey,PeterBorwein和SimonPlouffe共同发表了下面的圆周率计算公式（简称BBP公式）

该公式的最大优点在于：经后来人将该公式变形后打破了传统的计算方法，可以直接计算圆周率的任意第n位数，而不是先计算前面的n-1位数。1997年，FabriceBellard发表了一个比BBP算法更快的公式从而，大大降低了圆周率近似值的计算量.电子计算机的出现使π值计算有了突飞猛进的发展。1949年美国马里兰州阿伯丁的军队弹道研究实验室首次用计算机（ENIAC）计算π值，一下子就算到2037位小数，突破了千位数。1989年美国哥伦比亚大学研究人员用克雷－2型和IBM－VF型巨型电子计算机计算出π值小数点后4.8亿位数，后又继续算到小数点后10.1亿位数，创下最新的纪录。

2010年1月7日——法国一工程师将圆周率算到小数点后27000亿位。

2010年8月30日——日本计算机奇才近藤茂利用家用计算机和云计算相结合，计算出圆周率到小数点后5万亿位。π的计算方法古典方法（略）级数逼近法数值积分法随机模拟法π的幂级数计算方法思想：利用一些特殊函数的幂级数展开式计算圆周率.例如编写下面的程序:n=10;%选择展开式的次数s=0;digits(22);%定义计算过程中的精度fork=1:ns=s+4*(-1)^(k+1)/(2*k-1);endvpa(s,20)%定义显示精度为20位上机实验:

请按照上述公式分别编程计算，并与就精度与叠加次数进行比较，能得出怎样的结论？函数的泰勒展开式的matlab指令

MATLAB软件求一元函数泰勒展式的命令是taylor,具体使用格式如下：taylor(f)求函数f的5阶Maclaurin展开式；taylor(f,n)求函数f的n-1阶Maclaurin展开式；taylor(f,n,a)求函数f在x=a处的n-1阶泰勒展开式；

taylor(f,a)求函数f在x=a处的5阶泰勒展开式；对Taylor级数，|x|越小，级数收敛速度越快使用这种算法收敛速度要快的多。Machin公式英国天文学家麦琴（Machin）于1706年给出，他利用这个公式求得了π的第100位小数。利用类似方法，到1948年英国的弗格森和美国的伦奇共同发表了π的第808位小数，成为人工计算圆周率值的最高纪录。π的数值积分计算方法根据上述微积分公式，将求解圆周率转换为数学上经常遇到的一类问题——数值积分问题。

实际中常遇到这样的数值积分，被积函数不是由解析式给出，或者虽然是解析式但求其原函数十分困难。这时就要考虑定积分的近似计算了。定积分的近似计算梯形法

分割区间，用梯形面积代替曲边梯形面积抛物线法（Simpson法）

分割区间，在每个区间添加其中点，用过两端点及中点的抛物线面积代替曲边梯形面积。n=50;%定义等分积分区间数，可以更改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梯形法等分区间数n103.14242598500110203.14180098689309503.141625986923001003.141600986923125003.1415929869231210003.1415927369231350003.14159265692313下表给出的是不同等分区间数计算出的

圆周率π的近似值抛物线法Matlab提供了积分函数quad()clear;f=inline('4./(1+x.*x)');a=0;b=1;n=1;z=quad(f,a,b);fprintf('z=%.10f\n',z)按照实际问题所遵循的概率统计规律，用电子计算机进行直接的抽样试验，然后计算其统计参数。蒙特卡洛方法也可以人为地构造出一个合适的概率模型，依照该模型进行大量的统计实验，使它的某些统计参量正好是待求问题的解。在物理学的各种各样问题中得到广泛的应用。1）蒙特卡洛方法圆周率π的随机模拟计算方法cs=0n=500%随机取点数fori=1:na=rand(1,2);ifa(1)^2+a(2)^2<=1cs=cs+1endend4*cs/n１１０正方形中任取一点，落入四分之一圆中的概率为正方形中任取n个点，其中m个点落入四分之一圆中，则当m、n充分大时，π的蒙特卡洛计算方法

算法优点：简单易行；缺点：算法收敛的速度比较慢．

在精度要求不是很高的情况下，这种取随机数进行数据模拟的方法还是有一定的实用价值，通过简单编程我们可以模拟出许多数学现象．

例:计算定积分1）随机地向正方形内投点，最后统计落在曲线下的点数M；2）当总的掷点数N充分大时，N/M就近似等于积分值I。上机实验：基于关系式利用蒙特卡洛方法近似计算保密通讯无论在军事、政治、经济还是日常生活中都起着非常重要的作用。为了将信息传递给己方的接受者，同时又要防止他人（特别是敌人）知道信息的内容，必须将要传递的信息（明文）加密，变成密文后发送出去，这样，即使敌方得到密文也看不懂，而己方的接受者收到密文后却可以按照预先定好的方法加以解密。Hill密码密码可分为古典密码和现代密码本实验主要介绍古典密码的加密与破译原理，同时介绍如何用Matlab编程来实现加密、解密和破译过程。问题背景和实验目的

古典密码：以字符为基本加密单元；

现代密码：以信息块为基本加密单元。加密信息传递过程明文（信息）加密器密文密文明文（信息）解密器普通信道发送敌方截获破译发送方接收方Hill2密码的加密过程Hill2密码中所用的数学手段是矩阵运算。

加密过程：①

将汉语拼音的26个字母

与0到25

之间的整数建立一一对应关系，称为字母的表值，然后根据明文字母的表值，将明文信息用数字表示。ABCDEFGHIJKLM12345678910111213NOPQRSTUVWXYZ1415161718192021222324250设通讯双方所给出的26个字母的表值如下：注：这里假定明文中只使用26个大写字母Hill2密码的加密过程②

选择一个二阶可逆整数方阵A，称为Hill2密码的加密矩阵，它是加密体制的“密钥”，是加密的关键，仅通讯双方掌握。③

将明文字母分组。Hill2

使用的是二阶矩阵，所以将明文字母每2个一组（可以推广至Hilln密码）。查出每个字母的表值，这样，每组字母构成一个二维列向量。若最后仅剩一个字母，则补充一个没有实际意义的哑字母（哑元）。这样使得每组都有2个字母，④令=

A

，由的两个分量反查字母表值表，得到相应的两个字母，即为密文字母。Hill2加密举例例：设明文为“KMXYSXX”（昆明学院数学系），试给出这段明文的Hill2

密文。其中加密矩阵为

将明文字母分组：

KMXYSXXX最后的一个字母X为哑字母，无实际意义。解:ABCDEFGHIJKLM12345678910111213NOPQRSTUVWXYZ1415161718192021222324250查表得每组字母的表值，得到4个二维列向量：将上述4个二维向量左乘密钥矩阵A

得：作模26运算，将所有的数都化为0到25之间的整数：Hill2

加密举例反查字母表值得每个向量对应的字母组为：KMXYSXXKMVWOTTTHill2加密Hill2

加密举例ABCDEFGHIJKLM12345678910111213NOPQRSTUVWXYZ1415161718192021222324250KMVWOTTT问题：怎样解密？明文字母查表值分组一组向量加密矩阵左乘一组新的向量反查表值密文Hill2

加密过程模运算先查出密文字母“KMVWOTTT”所对应的向量：在模运算下解方程组：

A

=

解密：加密的逆过程，将加密过程逆转回去即可。Hill2

解密过程上面的向量是由经过模26运算得来的，现在的问题是怎样逆转回去？例：怎么得到密文“KMVWOTTT

”的原文模m可逆记定义1：设A为定义在集合Zm上的n阶方阵，若存在一个定义在Zm上的方阵B，使得

则称A模m可逆，B为A的模m逆矩阵，记为定义2：设a

Zm，若存在b

Zm使得ab=1(modm)，则称b

为

a

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

Zm中的数命题：定义在集合Zm

上的n

阶方阵A模m

可逆的充要条件是：m和

det(A)

无公共素数因子，即m

与det(A)

互素。Hill2密码的加密矩阵必须满足上述条件。m=26m的素数因子只有2和13定义在Z26上的方阵A模26可逆的充要条件是：模m可逆det(A)不能被2和13整除

问题：是否

Zm中所有的数都存在模m倒数？a

存在唯一的模m倒数a

与m无公共素数因子Z26中具有模26

倒数的整数及其模26倒数表a1357911151719212325a-11921153197231151725m=26时的Matlab程序见教材第121页fuluA.m

可以用来判断一个2阶矩阵在模26运算下是否可逆，并在可逆的前提下给出其模26逆矩阵。模26可逆

思考：如何用Matlab编程来找出所有模m

倒数的整数及其模m倒数？在模运算下解方程组：A

=

Hill2

解密过程？

问题：如何计算

？模m逆矩阵的计算设B=kA*为A

的模26逆，其中k

为待定系数A*为A

的伴随矩阵本计算方法可推广到求矩阵A

的模m逆矩阵Hill2

解密过程设加密矩阵用B左乘密文对应的向量得：模26运算后得：查表后得明文分别为：KMXYSXXXMatlab的Hill2加密程序见附录B，相应解密程序见附录C。？Hill2加密过程总结①

通讯双方确定加密矩阵

(密钥)和字母的表值对应表②

将明文字母分组，通过查表列出每组字母对应的向量

③

令=

Amod(m)

，由

的分量反查字母表值表，

得到相应的密文字母若明文只含奇数个字母，则补充一个哑元Hill2解密过程总结①将密文字母分组，通过查表列出每组字母对应的向量②

求出加密矩阵A

的模m逆矩阵

B③

令=

B

mod(m)

，由

的分量反查字母表值表，

得到相应的明文字母甲方收到乙方（己方）的一个密文信息，内容为：Hill2解密举例WKVACPEAOCIXGWIZUROQWABALOHDKCEAFCLWWCVLEMIMCC按照甲方与乙方的约定，他们之间采用Hill2密码，密钥为，字母表值见下表，问这段密文的原文是什么？ABCDEFGHIJKLM12345678910111213NOPQRSTUVWXYZ1415161718192021222324250Hill2解密举例①将密文字母分组，通过查表列出每组字母对应的向量②

求出加密矩阵A

的模26逆矩阵③

用B

左乘每组密文字母组成的向量，然后再反查字母表值表，得到相应的明文字母序号分组密文密文表值明文表值分组明文1WK2311721GU2VA22149DI3CP316114AN4EA51139MI5OC153131MA6IX924198SH序号分组密文密文表值明文表值分组明文7GW723925IY8IZ9090IZ9UR211896IF10OQ15172123UW11WA23159EI12BA21109JIHill2解密举例序号分组密文密文表值明文表值分组明文13LO121525BE14HD841410NJ15KC11391IA16EA51139MI17FC6341DA18LW12231425NY序号分组密文密文表值明文表值分组明文19WC233211UA20VL2212144ND21EM513513EM22IM913913IM23CC3311AAHill2解密举例即：“古典密码是以字符为基本加密单元的密码”GUDIANMIMASHIYIZIFUWEIJIBENJIAMIDANYUANDEMIMAAWKVACPEAOCIXGWIZUROQWABALOHDKCEAFCLWWCVLEMIMCC原文Hill2解密举例密文ABCDEFGHIJKLM12345678910111213NOPQRSTUVWXYZ1415161718192021222324250经分析该密文是用Hill2密码加密，且密文(U,C)和(R,S)分别对应明文(T,A)和(C,O)，问能否破译这段密文？Hill2密码破译举例MOFAXJEABAUCRSXJLUYHQATCZHWBCSCP我方截获一段密文猜测密文是由26个字母组成，即m=26，

经破译部门通过大量的统计分析和语言分析确定表值破译这段密文的关键是找到“密钥”和字母对应的表值密文(U,C)和(R,S)分别对应明文(T,A)和(C,O)Hill2密码破译举例查字母表值PCP、C

模26可逆可唯一确定加密矩阵AHill2密码破译举例注：这里的运算都是在模运算意义下进行HEWILLVISITACOLLEGETHISAFTERNOON相应的Matlab程序见附录D

得到加密矩阵的模26逆矩阵后，根据前面的解密方法即可得密文的原文Hill2密码破译举例饮酒驾车模型在2003年全国道路交通事故死亡人数中，饮酒驾车造成的占有相当比例，为此，国家发布了新的《车辆驾驶人员血液、呼气酒精含量阈值与检验》国家标准。在新标准下，大李在中午12点喝了一瓶啤酒，下午6点检查时符合标准，接着晚上又喝了一瓶，但凌晨2点检查时却被定为饮酒驾车，为什么喝同样多的酒，两次检查结果不一样？药物在体内的分布与排除

药物进入机体形成血药浓度(单位体积血液的药物量)

血药浓度需保持在一定范围内——给药方案设计

药物在体内吸收、分布和排除过程——药物动力学

建立房室模型——药物动力学的基本步骤

房室——