计算机仿真技术实验指导书.doc

计算机仿真技术实验指导书黄剑航 编 莆田学院电子信息工程学系2011年8月 前言计算机仿真技术是利用计算机科学和技术的成果建立被仿真对象的模型，并在某些实验条件下对模型进行动态实验的一门综合性技术。它具有高效、安全、受环境条件的约束较少、可改变时间比例尺等优点，已成为分析、设计、运行、评价、培训系统（尤其是复杂系统）的重要工具，是电子信息本科专业的专业必修课。计算机仿真离不开基本的专业理论知识和计算机编程，因此学习本门课程可以使同学们巩固相关专业领域的基本知识、原理和方法，训练编程能力，掌握仿真的基本概念，并将仿真技术应用在专业领域的分析，问题解决上面，为今后进一步的学习、科研打下基础。本课程注重实践和理论结合，在实验室上机完成全部课程。结合课程和专业特点，实验指导书选取了MATLAB语言基础、基本算法仿真、趣味算法仿真、信号与系统仿真和Simulink仿真这几部分主要内容，作为学习仿真技术的主要内容。希望同学们在实验课前要做好预习工作，掌握基本方法原理，实验课堂上做好每个实验，对实验结果要做好必要的记录和保存，课后认真完成实验报告。只要认真对待每一次实验练习，相信同学们在专业理论知识和编程仿真实践上真正会有所收获。目录第1部分 MATLAB基础. 1实验1.1 MATLAB入门 . 1实验1.2 MATLAB数组及矩阵运算 . . . .3实验1.3 MATLAB编程应用.6实验1.4 MATLAB绘图. .10实验1.5 MATLAB匿名函数. 11第2部分 基本算法实现.16实验2.1 数值积分算法.16实验2.2 插值与拟合.18第3部分 趣味算法实现.23实验3.1 Hill密码.23实验3.2 随机模拟 .28第4部分 信号与系统仿真.34实验4.1 连续时间信号在MATLAB中的表示. 34实验4.2 连续时间信号在MATLAB中的运算.40实验4.3 傅里叶变换及其性质.45实验4.4 连续时间LTI系统的时域分析.51第5部分 Simulink建模仿真.57实验5.1 Simulink仿真基础.57实验5.2 连续系统建模仿真.61第1部分 MATLAB基础实验1.1 MATLAB入门1. 实验目的熟悉MATLAB工作环境和基本操作。2. 实验内容2.1熟悉MATLAB工作环境，掌握clc，clear，cd这几个常见命令。 图1 MATLAB桌面，具体布局可能因机器和版本的不同而会有变化2.2 假设x=3，y=4。利用MTLAB 计算下列表达式： ， 要求：（1）在命令窗口直接输入，得到运算结果；（2）编写M文档并在命令窗口执行，然后用工作空间平台查看有哪些变量在当前工作区。2.3 一小球从空中下落的位移公式为： 利用MATLAB 计算小球在t=5s 时的位置，已知：. 2.4 通过以下两种方式得到关于exp 函数的帮助 ：（1）在命令窗口中输入help exp 命令。（2）运用帮助空间窗口。2.5计算函数在为如下取值时候的函数值，取4.9 3.2, 100, 1.5, 9.75, 2.56，16， 4.9， 10。3. 思考题3.1 MATLAB 的命令窗口的作用是什么？编辑/调试窗口？图像窗口？ 3.2 列出几种不同的得到MATLAB 帮助的方法。 3.3什么是工作区？在同一工作区内，你如何决定它里面存储了什么？ 3.4 你怎样清空MATLAB 工作区内的内容？ 实验1.2 MATLAB数组及矩阵运算1. 实验目的熟悉MATLAB数组及矩阵运算。 2. 实验内容2.1 确定下列数组的大小，通过whos 或工作空间窗口（The workspace browser）检查你 的答案。注意在本练习中后面的数组可能要用到前面数组的定义。 (1) u=10 20*i 10+20(2) v=-1;20;3(3) w=1 0 -9;2 -2 0;1 2 3(4) x=u v(5) y(3,3)=-7(6) z=zeros(4,1) ones(4,1) zeros(1,4)(7) v(4)=x(2,1)2.2 w(2,1)的值是多少？ 2.3 x(2,1)的值是多少？ 2.4 y(2,1)的值是多少?2.5 当语句3(7)执行后，v(3)的值是多少？ 2.6 c 数组的定义如下，写出下面子数组的内容。 c =1.1000 -3.2000 3.4000 0.60000.6000 1.1000 -0.6000 3.10001.3000 0.6000 5.5000 0(1) c(2,:) (2) c(:,end) (3) c(1:2,2:end) (4) c(6)(5) c(4:end) (6) c(1:2,2:4) (7) c(1 4,2) (8) c(2 2,3 3)2.7 当赋值语句执行后，下列数组的内容是多少？ (1) a=1 2 3; 4 5 6; 7 8 9;a(3 1,:)=a(1 3,:);(2) a=1 2 3; 4 5 6; 7 8 9;a(1 3,:)=a(2 2,:);(3) a=1 2 3; 4 5 6; 7 8 9;a=a(2 2,:);2.8 当数组执行后，下列数组a 的内容是多少？ (1) a=eye(3,3);b=1 2 3;a(2,:)=b;(2) a=eye(3,3);b=4 5 6 ;a(:,3)=b;(3) a=eye(3,3);b=7 8 9;a(3,:)=b(3 1 2);(4) a=eye(3,3);b=7 8 9;a(3,:)=b(3 1 2); 2.9 假设a,b,c 和d 的定义如下： 分别运行出下列表达式的运算结果，并思考。 (1) a + b (2) a .* c (3) a * b (4) a * c(5) a + c (6) a + d (7) a .* d (8) a * d3. 思考题3.1 数组,矩阵,向量有什么区别?3.2 回答关于下列矩阵的有关问题 (1) C 的大小是多少？ (2) C(2,3)的值是多少？ (3) 列出值为0.6 的元素的下标实验1.3 MATLAB编程应用1. 实验目的学习函数的编制，掌握MATLAB的实际应用。通过不同的程序结构和不同的实际应用问题，掌握MATLAB的编程方法。2. 实验内容2.1 画出函数 在区间0,10的波形。 用MATLAB 编辑器创建一个新的M 文件，把上面的语句写入这个文件并命名为 ： test1.m。然后在命令窗口中输入test1 执行这个文件。看得到什么结果？2.2 编写一个MATLAB 程序，要求输入圆柱体的半径和高，然后输出体积。2.3 一个程序实例学习： （温度转换）设计一个MATLAB 程序，读取一个华氏温度的输入，输出开尔文温度。 华氏温度和开尔文温度的转换关系式可在物理学课本中找到。其关系式为： 在物理学参考书中举了一些例子，我们可以用来检验我们程序是否正确。例如 我们设计程序的步骤如下 ：（1）提示用户键入华氏温度值 （2）读取输入值 （3）通过关系式转换为开氏温度 （4）输出结果，结束 我们将会用input 函数输入华氏温度，用fprintf 函数输出结果。 % Script file:temp_conversion.m% Purpose:% To convert an input temperature from degrees Fahrenheit to% an output temperature in kelvins.% Record of revisions:% Date Programmer Description of change% = = =% 12/01/97 S.J.Chapman Original code%Define variables:% temp_f -Temperature in degrees Fahrenheit% temp_k -Temperature in kelvins%Prompt the user for the input temperature.temp_f=input(Enter the temperature in degrees Fahrenheit:);%Converttokelvins.temp_k=(5/9)*(temp_f-32)+273.15;%Writeouttheresult.fprintf(%6.2f degrees Fahrenheit = %6.2f kelvins.n,.temp_f,temp_k);我们输入上面的例子中的华氏温度值，以检测程序的正确性。注意用户的输入值已用黑 体字标出。 temp_conversionEnter the temperature in degrees Fahrenheit:212212.00 degrees Fahrenheit = 373.15 kelvins. temp_conversionEnter the temperature in degrees Fahrenheit:-110-110.00 degrees Fahrenheit = 194.26 kelvins.这个结果和物理教科书的结果相同。 2.4 编写一个函数文件，计算出坐标系中用户指定两点(X1,Y1)和(X2,Y2)之间的距离。要求有输入、输出及其相关提示。2.5 双曲余弦的定义如下： ，编写一个程序，计算出用户指定的x 的值对应的双曲余弦值。用这个程序计算双曲余弦值的若干值，并和MATLAB 中的内建函数cosh(x)得到的值比较看看是否完全相同。并用MATLAB 打印出这个函数的图象。2.6 电子工程：负载的最大输出功率一个内阻Rs50，电动势V120V 的电源驱动一个负载RL。当RL 为多少时，RL 的功率最大？在这种情况下，功率为多少？画以RL 为自变量的RL 功率图。2.7 利用公式求的近似值，直到最后一项的绝对值小于为止。2.8 Fibonacci(斐波纳契）序列的元素满足Fibonacci 规则： 且；现要求该序列中第一个大于20000 的元素，并指明该元素是序列的第几项。2.9 在田径比赛中，一个身高为一米八零的铅球运动员，大概以多大的角度（和水平方向夹角）推铅球，才能使铅球推得最远，并求出最远距离。不计空气阻力，假设铅球出手点和运动员高度相等，且铅球出手瞬间初始速度大小为14m/s，重力加速度取g=10m/s2。2.10 打印出所有的“水仙花数”，所谓“水仙花数”是指一个三位数，其各位数字立方和等于该数本身。例如153是一个水仙花数，满足15313+53+33。2.11 有一个函数：，写一程序，输入x，输出y值。2.12 输入3个整数，要求按大小输出。3. 思考题3.1 MATLAB 启动时，当前工作路径是什么？在任意路径下新建一个名字为mynewdir的新文件夹， 把当前工作路径改为mynewdir，然后打开一个编辑窗口，增加以下语句： % create an input array from -2*pi to 2*pit = -2*pi:pi/10:2*pi;% calculate|sin(t)|x=abs(sin(t);%plot resultplot(t,x);把此文件以test2.m 为文件名保存，然后在命令窗口中输入test2 执行此文件，观察结果。关闭画图窗口，然后返回到原始工作路径，再在命令窗口中输入test2。看有何情况发生，为什么？3.2 MATLAB中的函数文件和脚本文件有什么区别？实验1.4 MATLAB绘图1. 实验目的掌握MATLAB的几种绘图方法和相关命令，包括plot，meshgird等命令。掌握结果图的一些参数设置（如坐标范围，坐标轴标示，曲线类型等等）。2. 实验内容2.1 绘制正切函数 在0x1 上的图象，并为x 轴和y 轴添加标签。增量设为0.1。2.2 显示同样的图象，把sinx添加到第二个图形中作为第二条曲线。2.3 使用plot3 函数绘制曲线、和的图象，不要添加坐标轴标签，但要打开网格。2.4 绘制以参数方程表示的空间曲线。2.5 绘制以极坐标方程表示的心脏线。2.6 分别在，的条件下使用命令fplot做出函数的图形。2.7 试在矩形区域，上分别绘制函数与对应的三维网格表面图和三维曲面图。3. 思考题3.1 执行完 x,y=meshgrid(-10:0.5:10)之后的x、y值是多少？3.2 说明在计算诸如 z=x2+y2 这样的f(x,y)函数时，必须先使用meshgrid函数的原因。实验1.5 MATLAB匿名函数1. 实验目的MATLAB7.0以后引入了一种新的函数类型-匿名函数(anonymous function)。匿名函数不仅完全实现了内联函数(inline function)的所有功能，而且一些用内联函数实现起来不大方便或者不简洁的功能，用匿名函数也非常容易实现。除此之外，匿名函数的调用效率要远高于内联函数。本实验的主要目的就是熟悉匿名函数的用法。2. 匿名函数的定义2.1 匿名函数的基本定义在程序编写过程中，经常需要计算诸如，等这样一些数学表达式的值，如何方便地构造这些函数？写成M文件再调用肯定可以，但是从简洁以及管理的方便性考虑这不是好办法。这里引入了匿名函数，基本定义如下：Fhandle=(arglist) expr其中，expr是具体的函数表达式，arglist是指定的函数的自变量。具体见下面示例：f=(x) x.2;fx=f(1:10) fx = 1 4 9 16 25 36 49 64 81 100g=(x,y) x.2+y.2;gxy=g(1:10,2:11) gxy = 5 13 25 41 61 85 113 145 181 2212.2 匿名函数的种类这里按照自变量的个数以及层数可以分为以下几种：单变量匿名函数、多变量匿名、单重匿名函数以及多重匿名函数。1.单变量匿名函数单变量匿名函数是最简单的匿名函数，只含有一个变量，如：f=(x) x.2就是单变量匿名函数，除此之外，含有参数，参数值已知的单个自变量的匿名函数也是单变量匿名函数，比如：a=10;b=15;f=(x) a*x+b;f(1:5) ans = 25 35 45 55 652.多变量匿名函数多变量匿名函数含有两个或者两个以上的自变量，如g=(x,y) x.2+y.2。同单变量匿名函数一样，多变量匿名函数也可以有自己的参数值。比如：a=1;b=2;g=(x,y) a*x+y.b;g(1:5,1:5) ans = 2 6 12 20 303.单重匿名函数到目前为止，上面列的匿名函数，无论单变量还是多变量匿名函数，都属于单重匿名函数。这类匿名函数的特点是：只有一个符号引导，符号之后就是具体的函数表达式。自变量输入单重匿名函数后，得到的就是具体的数值结果。除了单重匿名函数外，还有二重、多重匿名函数。这些多重匿名函数在参数传递方面非常方便。4.多重匿名函数下面以二重匿名函数为例来介绍多重匿名函数。例子如下：f=(a,b) (x) a*x+b f = (a,b)(x)a*x+b其中a，b是外层变量，x是内层变量。可以这样理解这个表达式：每个符号后面括号里的变量作用域一直到表达式结尾。这样a，b的作用域就是(x) a*x+b，而x的作用域就是 a*x+b。因此，若有个表达式为fab=f(a,b)，意思是这是一个在给定的参数a，b情况下，以x为变量的单层的单变量匿名函数。可以利用functions函数观察建立的匿名函数的信息，如下：f=(a,b) (x) a*x+b f = (a,b)(x)a*x+b f23=f(2,3) f23 = (x)a*x+b3. 应用实例3.1 匿名函数可以非常方便地表示所求方程，并供fzero等求解函数调用。求下列方程的根：求解代码如下： f=(x) exp(x)+x2+x(sqrt(x)-100f = (x)exp(x)+x2+x(sqrt(x)-100 x0=fzero(f,3)x0 = 4.1635f(x0) ans = 2.8422e-0143.2 显式表示下列y关于x的隐函数：，利用匿名函数，可以在MATLAB中显式写出y和x的关系式如下：y=(x) fzero(y) (exp(y)+xy)(1/y)-x2*y,1)这样，对于任意的x，只需调用，就能求得相应的y值，如：y1=y(1)y1 = 2.7779y2=y(2)y2 = 1.10553.3 显式表示下列z关于x，y的隐函数：借助匿名函数，可以在MATLAB中显式写出z关于x，y的关系式如下： z=(x,y) fzero(z) z-sin(z*x-0.5)2+x*2*y2-z/10)*.exp(-(x-0.5-exp(-y+z)2+y2-z/5+3),rand);其中，fzero函数求解该隐函数的初值为随机值rand。利用上述匿名函数，可以画出z关于x，y的图形，如图1所示：X,Y=meshgrid(-1:0.1:7,-2:0.1:2);Z=arrayfun(x,y) z(x,y),X,Y);surf(X,Y,Z)xlabel(fontsize15fontnametimes new romanx,color,b)ylabel(fontsize15fontnametimes new romany,color,b)zlabel(fontsize15fontnametimes new romanz,color,b)title(fontsize15fontnametimes new romanz(x,y),color,b) 图1 实例3.3结果图3.4 求，k取值0,5区间不同值时的函数值，并画出函数图像。分别用inline内联函数类型和匿名函数类型实现。该例对于一些复杂的、无显式积分表达式的带参数积分问题具有通用性。解决代码如下：tic;k=linspace(0,5);y1=zeros(size(k);for i=1:length(k) kk=k(i); fun=inline(sin(,num2str(kk),*x).*x.2); y1(i)=quadl(fun,0,5);endtime=toc;disp(Inline method costs time:,num2str(time),seconds!); tic;f=(k) quadl(x) sin(k.*x).*x.2,0,5);kk=linspace(0,5);y2=zeros(size(kk);for ii=1:length(kk) y2(ii)=f(kk(ii);endtime=toc;disp(Anonymous function method costs time:,num2str(time),seconds!); plot(kk,y2);xlabel(k);ylabel(f(k);grid on;得到的运行结果图如下图2所示：图2 实例3.4结果图第2部分 基本算法实现实验2.1 数值积分算法 1. 实验目的为了研究连续系统的性质，需要建立连续系统的数学模型，然后进行仿真。连续系统通常可以用常微分方程来描述。对于连续系统的仿真，就是用计算机求解这些数学模型，也就是对微分方程求数值解。各种数值积分算法在常微分方程的求解中起到了关键作用。通过算法仿真实验，达到对数值积分算法本质更好的理解和应用的目的。2. 算法原理欧拉法，改进的欧拉法，龙格-库塔法，线性多步法等，详见课本。3. 实验内容3.1 设一微分方程为，初始条件。试编写一个程序用欧拉法求其数值解。3.2 已知方程的解析解为，且时，。取计算步长为，试用欧拉法，改进欧拉法和四阶龙格-库塔法求取值是05之间的数值结果，并和解析解做一比较。3.3 对于线性微分方程，分别用改进欧拉法和龙格-库塔法进行求解，然后将结果与用MATLAB中现成函数求得的结果进行比较。3.4 设连续系统的数学模型用下列微分方程描述：试计算当时，的时间响应。4. 思考题4.1 对于高阶微分方程的数值积分求解，如何进行预处理，转换成为适合欧拉法或者龙格-库塔法求解的一阶形式？4.2 已知一阶系统的微分方程为，初始条件，取仿真步长，分别用欧拉法，改进欧拉法和龙格-库塔法计算该系统仿真第一步的值。实验2.2 插值与拟合 1. 实验目的理解插值与拟合各自的基本定义，算法原理。实现算法在实际中的具体应用。掌握相关的MATLAB函数用法2. 算法原理2.1 插值方法原理在一些实际应用中，会得到一些数据，通常将这些数据列成数据表格的形式。这种用数据表格形式给出的函数称为列表函数，即，其中点成为结点。根据函数已有的数据表格来计算函数在一些新的点处的函数值，这就是插值所要解决的问题。一种最简单的插值方法就是拉格朗日插值，利用拉格朗日多项式即可计算相应区间任意点的插值结果。候拉格朗日插值多项式也被称为拉格朗日公式，如下： 拉格朗日公式常用于以下两种简单的情形。（1）线性插值设函数在处的函数值分别为，则线性插值公式为利用线性插值公式就可以计算在之间任意一个新点x处的近似值。（2）二次（抛物线）插值设函数在处的函数值分别为，则二次插值公式为如果这三点不在一条直线上，则上述插值公式得到的是一个二次函数，通过这三点的曲线是抛物线。利用三点的二次插值公式就可以计算在所在区间的任意一个新点x处的近似值。2.2 拟合方法原理给定平面上的一组数据，要求确定一个初等函数，使给定的数据点（在某种距离下）与函数表示的曲线整体上最为接近，这类问题称为数据拟合。设有经验公式（这里x和c均可为向量），要求根据数据，确定参数c，这样的问题称为曲线拟合，其基本原理是最小二乘法，即求使误差平方和达到最小的c。MATLAB中有关拟合的命令如下：p=polyfit(x,y,k) 用k次多项式拟合向量数据(x,y)，返回多项式的降幂系数向量。yi=polyval(p,xi) 根据多项式p计算出在xi处插值函数值yi.c=lsqnonlin(fun,c0) 使用迭代法搜索最优参数c，其中fun是以参数c（可以是向量）为自变量的表示误差 (x,y是数据向量)的函数，c0是c的近似值，作为迭代初始值。c=lsqcurvefit(fun,c0,x,y) 根据数据向量x,y 确定经验公式（c可以是向量）中的参数c。这里fun为函数，c0是c的初始值。与lsqnonlin不同，使用lsqcurvefit时数据x,y从外部输入。3. 实例分析3.1 以下是美国人口两个世纪以来的统计数据，试以此预测2010年，2020年的美国人口。年份180018101820183018401850186018701880人口/百万5.37.29.612.917.123.231.438.650.2年份189019001910192019301940195019601970人口/百万62.976.092.0106.5123.2131.7150.7179.3204.0年份198019902000人口/百万226.5251.4275.0解：先做散点图，为了方便，也为了避免数据相差太大造成误差，年份均减去1800.t=0:10:200;y=5.3 7.2 9.6 12.9 17.1 23.2 31.4 38.6 50.2 62.9 76.0 92.0. 106.5 123.2 131.7 150.7 179.3 204.0 226.5 251.4 275.0;subplot(2,2,1);plot(t,y,bo);axis(0,300,0,400);text(20,350,Population);从数值上看给定数据接近抛物线或者指数曲线。图1 美国人口散点图及预测结果图（1）考虑用二次多项式拟合，代码如下：p=polyfit(t,y,2)Q1=sqrt(sum(y-polyval(p,t).2)ti=0:2:220;yi=polyval(p,ti);yj=polyval(p,210,220)subplot(2,2,2)plot(t,y,bo,ti,yi,r,210,220,yj,ro)text(20,350,Bipoly);结果：拟合曲线见图1右上角，均方误差=12.5425，预测人口=303.0567 331.8160.（2）采用马尔萨斯模型，人口按指数增长，因此可设人口曲线为，其中N0, r是待拟合的参数，相关代码如下：fun1=(c,t) c(1)*exp(c(2)*t);c= lsqcurvefit(fun1,6,0.1,t,y)Q2=sqrt(sum(y-c(1)*exp(c(2)*t).2)ti=0:2:220;yi=c(1)*exp(c(2)*ti);yj=c(1)*exp(c(2)*210,220)subplot(2,2,3);plot(t,y,bo,ti,yi,r,210,220,yj,ro)text(20,350,Exp curve);结果：拟合曲线见图1左下角，均方误差=48.5696，预测人口=338.3433 389.2984.（3）采用Logistic模型，人口增长规律为 ，其中Nm,N0,r是待拟合的参数，相关代码如下：fun2=(c,t) c(1)./(1+(c(1)/c(2)-1)*exp(-c(3)*t);c=lsqcurvefit(fun2,500,6,0.2,t,y)Q3=sqrt(sum(y-c(1)./(1+(c(1)/c(2)-1)*exp(-c(3)*t).2)ti=0:2:220; yi=c(1)./(1+(c(1)/c(2)-1)*exp(-c(3)*ti);yj=c(1)./(1+(c(1)/c(2)-1)*exp(-c(3)*210,220)subplot(2,2,4)plot(t,y,bo,ti,yi,r,210,220,yj,ro)text(20,350,Logistic);结果：拟合曲线见图1右下角，均方误差=19.8133，预测人口=292.9388 311.8921从以上结果看，二次多项式拟合均方误差最小，Logistic模型次之，马尔萨斯模型最差，但从机理上考虑，二次多项式和马尔萨斯模型确定的人口最终会趋于无穷大，不合理，它们都只能用作近期预测。相对来说，Logistic模型较为合理。.4. 实验内容4.1 设函数在直角坐标平面内四个点(0,0), (0,1), (1,0), (1,1)的所对应的函数值分别为1，3，2，5，试用两次线性插值（双线性插值）求坐标为(0.2,0.4)该点处的函数值。4.2 设x0=0.4,x2=0.8,而x1可取0.5,0.6,0.7中任何一点，根据下表使用二次拉格朗日插值求sin0.6389的近似值时，问如何选节点x1，可使其近似值的误差较小？试计算说明。x0.40.50.60.70.8sinx0.38940.47940.56460.64420.71744.3 用电压U=10v的电池给电容充电，电容上时刻t的电压为：，其中V0是电容的初始电压，是充电常数，试由下面一组t，U数据确定U0和。t/s0.51234579U/V6.366.487.268.228.668.999.439.635. 思考题5.1 说明双线性插值的原理。5.2 结合例子说明参数拟合时候如何选择参数的初始值。第3部分 趣味算法实现实验3.1 Hill密码 1. 实验目的了解加密和解密原理，掌握模运算下的矩阵运算以及数字与字符串之间的转换。2. 算法原理保密通信具有悠久的历史，现在更是被广泛应用于军事、经济、商业等各行各业。在保密通信中奖原信息称为明码，加密后的信息称为密码。如果不知道加密方法，一般人无法知道明文内容，这样就起到了保密的作用。加密和解密的过程可以抽象为一个数学模型。移位加密法 移位加密法是一种简单的加密方法。它通过将明文的字母按字母表中的次序平移若干位实现加密。如加密方法是字母平移5位，则明文字母和密文字母对应关系如下：明文字母 ABCDEFGHIJKLMNOPQRSTUVWXYZ密文字母 FGHIJKLMNOPQRSTUVWXYZABCDE明文THEOLYMPICGAMES 加密后就成了 YMJTQDRUNHLFRJX. 一般人不知道它是什么意思，这就起到了加密的作用。如果知道它是由字母表平移5位得来的，那就很容易获得原文。数字5是解开密码的一把钥匙，称为密钥。在移位加密法中，明文字母和密文字母之间的对应关系是固定的，这种加密方法可以通过分析词频，利用统计方法破译。下面要介绍的Hill密码是用矩阵运算实现加密的，它不保持明文字母和密文字母之间固定的对应关系，破译起来比移位法困难。模n运算 在模n运算下，参与运算的只有0，1，2，n-1这n个元素。规定两个元素的和（积）是按普通加法（乘法）运算的结果减去或加上n的某个倍数后得到的在0和n-1之间的数。例如在模15运算下4+13=2，4*13=7。通常记做4+13=2(mod 15),4*13=7(mod 15)。在模n运算下，元素m的加法负元是n-m。如果m*k=1(mod n)，则称k是m的乘法逆元。如在模15运算下，13的逆元是7，4的逆元还是4。值得注意的是在模15运算下，除了0没有逆元外，3，5及3，5的倍数6，9，10，12也没用逆元。可以证明：若m和n的最大公因子等于1，即gcd(m,n)=1，则在模n运算下元素m存在唯一的逆元。特别，当n是素数时，元素1,2,n-1都有逆元，此时数0,1,2,n-1在模n下可以进行加减乘除(0除外)四则运算，它们构成一个数域，称为有限域。模n下矩阵加法、乘法以及数乘矩阵运算与模n下数字的加法、乘法运算相似，先对矩阵做普通加法、乘法和数乘矩阵运算，然后所有数字以n取模即可。定义 两个方阵A、B如满足AB=BA=E(mod n)，则称矩阵A模n可逆，B称为A的逆矩阵，记做(mod n)。矩阵A的逆矩阵可以利用伴随矩阵，按求得。Hill加密法（1） 先将英文字母变换成数字；A B C D E F G H I J K L M N O P Q R S T U V W X Y Z 1 2 3 4 5 6 7 8 9 10 24 25 0称其为字母表的表值。（2） 将明文按n个字母分组，并用对应的数字取代字母，构成一个个n维向量。（3） 取一个在模n运算下可逆的n阶矩阵A，用A左乘（2）中的向量得到新向量；把新向量中的数字按（1）中的对应关系换回成字母即得到密文。（4） 按例中的方法求出可逆矩阵A在模n下的逆矩阵B，用B左乘（3）中得到的新向量则重新得到（2）的向量，从而得到原始明文。矩阵A是解密的关键，称为密钥。3. 实例分析3.1 在模26运算下，求矩阵的逆矩阵。解：由于det(A)=3，且(mod 26)，故3.2 采用Hill2加密法（即明文2个2个分组），加密矩阵取，对 “THE OLYMPIC GAMES” 进行加密，再将结果解密。解：（1）为了保留英文单词之间的空格，可以增加逗号“，”，句号“.”和空格“ ”当做字母，连同26个英文字母一起与数字1,2,28,0之间建立对应关系，其中126分别对应AZ，“，”对应27，“.”对应28，“ ”对应0。（2）将“THE OLYMPIC GAMES”两两分组，最后不足两个字母时用“空格”补足，并用数字代替，得（3）取矩阵，输入A=2 7;4 5;x=20 8; 5 0;15 12;25 13; 16 9; 3 0; 7 1;13 5; 19 0;y=mod(A*x,29)结果为：y = 9 10 27 25 8 6 21 3 9 4 20 4 20 22 12 4 19 18这个换成字母就是“IDJT,DYTHVFLUDCSIR”,这就是密文。注意，明文中的两个“E”，在密文中分别变成了“J”和“S”；而密文中的三个“D”对应明文中的“H”，“L”和“A”，它们并不一一对应。（4）如果知道密钥A，解密是一件很容易的事情。A是模29下的可逆矩阵，其逆矩阵，用B左乘（3）中的向量重新得到20 5 15 25 16 3 7 13 198 0 12 13 9 0 1 5 0换成字母就是“THE OLYMPIC GAMES”，即原始明文。以下是本题加密和解密的一个完整的MATLAB程序，程序名为Hill2.mclear;clc;close all;A=2 7;4 5;x=THE OLYMPIC GAMES % encryption x2=double(x);x3=(x2-64).*(x264)+(x2-17).*(x2=44)+(x2-18).*(x2=46)+(x2-32).*(x2=32);if rem(length(x2),2) =1 x3=x3,0;endx4=reshape(x3,2,length(x3)/2);x5=mod(A*x4,29);x6=x5(:);x7=(x6+64).*(x664)+(y2-17).*(y2=44)+(y2-18).*(y2=46)+(y2-32).*(y2=32);y4=reshape(y3,2,length(y3)/2);y5=mod(B*y4,29);y6=y5(:);if y6(end)=0 y6(end)=;endy7=(y6+64).*(y60)+(y6+17).*(y6=27)+(y6+18).*(y6=28)+(y6+32).*(y6=0);z=char(y7)4. 实验内容4.1 根据移位加密法编写一段加密和解密的程序，并将所编程序应用于明文 “HELLOWORLD” ，密钥取3.4.2 甲方收到与之有秘密通信往来的另一方的一个密文信息，密文内容为：JCFOAMOOLABUOHFAXOIKVLYWQSTPKPYABHROARXAMO按照甲方与乙方的约定，他们之间的密文通信采用Hill2密码，密钥为矩阵， 26个字母的表值如正文所示，问这段密文的原始明文是什么？实验3.2 随机模拟 1. 实验目的了解随机模拟的基本概念和原理，通过实验内容，掌握利用计算机进行随机模拟的方法。2. 背景介绍随机模拟是一种随机试验的方法，也称蒙特卡洛(Monte Carlo)方法。这种方法源于美国第二次世界大战期间研制原子弹的“曼哈顿计划”，该计划的主持人之一，冯诺依曼用驰名世界的赌城-摩纳哥的蒙特卡洛来命名这种方法，使它蒙上了一层神秘的色彩。设计一个随机试验，只要使一个事件的概率与某个未知数有关，然后通过重复试验，以频率近似表示概率，即可求出该未知数的近似解。现在，随着计算机的发展，已按照上述思路建立起一类新的方法-随机模拟方法。计算机产生的随机数都是按照某种确定的算法产生的，它遵循一定的规律，一旦初值确定，所有随机数也就随之确定，这显然不满足真正随机数的要求，因此我们称这种随机数为“伪随机数”。在仿真中，通常，它们能当做真正的随机数那样使用。符合各种统计规律的随机数可以通过MATLAB的随机数生成函数很方便地调用。3. 实例分析3.1 的模拟计算建立一个概率模型，它和有关，然后设计适当的随机试验，并通过这个试验的结果来确定。通过这个实验理解模拟计算的基本思想。大家知道圆周率的值本身没有解析解。现在用随机模拟的方法设计一种求的近似值的方法，并计算它的近似值。在平面上画一组距离为a的平行线，把长度为l(la)的针投到平面上，则针与平行线相交的概率为。这是著名的蒲丰投针问题，如图1所示。图1 蒲丰投针问题用x表示针的中点与最近一条平行线的距离，用表示针与此线的交角，显然样本空间为：以R表示边长为a/2和的长方形，针与平行线相交当且仅当，如图2所示。设在R中满足这个