MATLAB语言与数学实验综合实训报告13统计班20134046030张淑娟

《MATLAB语言与数学实验》综合实训报告

班 级： 13级统计 姓 名： 张淑娟 学 号： 20134046030 指导老师： 莫宏敏老师

吉首大学数学与统计学院

2016-6-27

《MATLAB语言与数学试验》实训大纲

一、课程基本信息

1、课程名称：MATLAB语言与数学实验 2、课程类别：选修

3、适用专业：信息与计算科学，统计学，金融工程 4、课程地位：专业选修课 5、总学时：一周

6、先修课程：数学基础课，计算机基础 二、实习目的和任务

开设数学实验实训的目的是使学生掌握数学实验的基本思想和方法,提高学生的动手能力。从实际问题出发，借助计算机，通过学生亲自设计和动手，体验解决问题的全过程，从实验中去探索、学习和发现数学规律，充分调动学生学习的主动性,培养学生的创新意识.培养学生运用所学知识,借助计算机并利用数学软件解决实际问题的能力，最终达到提高学生数学素质和综合能力的目的。

三、教学基本要求

1、掌握MATLAB语言概况及基本知识；

2、掌握符号计算和数值计算的使用，能进行交互状态下的程序编制； 3、掌握程序的设计与编写，能编程解决一些简单的数学问题；

4、学习matlab在线性代数，概率统计、微分方程、数值分析等方面的应用； 5、掌握用计算机语言作函数图形的方法及对图形的修饰、修改； 6、掌握MATLAB的优化工具箱的用法,能求解简单的数学规划问题。

四、课程内容与学时分配

本课程主要分为三部分：第一部分就是数学软件MATLAB的学习、掌握和应用，是利用软件进行数学实验的基础；第二部分是实验部分，围绕大学数学的基本内容，让学生充分利用计算机及软件的数值、统计计算和图形功能展示基本概念和理论，去体验如何发现、总结和应用数学规律；第三部分是综合部分，以学数学用数学为中心培养学生综合应用数学的能力（具体见实验内容）。

五、有关说明

1、《MATLAB语言与数学实验》实训教学是课堂理论讲授的补充。 2、考核形式：写出实训报告作为评分依据。

3、要求写出源程序及运行结果，用A4纸双面打印并按班级统一打包交电子版（个人电子版命名方式：班级—学号—姓名）。

实验内容

1（10分）、分别就x?[?10,10]和x?[?1,1]，作出函数z?cos(4x2?9y2)的图形并观察图形的区别。

解：x1=-10:0.1:10; x2=-1:0.1:1; y1=x1; y2=x2;

subplot(1,2,1)

z1=cos(4*x1.^2+9*y1.^2) plot3(x1,y1,z1); subplot(1,2,2)

z2=cos(4*x2.^2+9*y2.^2) plot3(x2,y2,z2)

110.50.500-0.5-0.5-110100-10-100-1110-1-10

2（10分）、编制函数文件hillt.m其调用格式为y=hillt（m,n），其中m,n为正整数，返回值为m*n的矩阵，且y(i,j)?1/(i?j)，且满足要求：

（1）y=hillt(n)等价于调用y=hillt(n,n).

（2）当输入参数大于2个或没有参数时，给出报错的信息“输入参数个数不对！”，当输出参数不是一个时给出报错信息“输出参数只有一个！”。 解：function y=hillt(m,n) if nargin==1 n=m; else

if nargin~=2

error('输入参数个数不对') end end y=[]; for i=1:m for j=1:n

y(i,j)=1/(i+j) end end

>> y=hillt()

错误使用 hillt (line 6) 输入参数不对

>> y=hillt(2,3) y =

0.5000 0.3333 0.2500

0.3333 0.2500 0.2000

>> y=hillt(3) y =

0.5000 0.3333 0.2500 0.3333 0.2500 0.2000 0.2500 0.2000 0.1667

3（10分）、试用二分法编程求方程x?sinx?1?0在[0, 1]小数6位。

解：f=inline('x-sin(x)-1','x'); a=2;b=3;e=1;k=0; while e>1e-6 c=(a+b)/2; if f(a)*f(c)<0 b=c;

elseif f(a)*f(c)>0 a=c; else a=c;b=c; end e=b-a k=k+1; end x=(a+b)/2; [k,x]

之间的根，精确到

0.510973

?12?A?4（15分）、调试课本126页附录B，模26意义下密钥矩阵为??时，对

03??任意输入的大写字符串的加密程序及附录C的相应解密程序，并写出

str=‘MYFRIENDWILLCOMETODAY’的密文，和str=‘ROOAJJKSLSWOXSFAWXH’的明文。

解： m=26;enmat=[1 2;0 3];demat=[1 8;0 9];ZERO=64;c=[];en=[]; astr=input('输入要加密的明文文字（全部为大写字母):') astr=input('输入要解密的密文密码（全部为大写字母):') while any(double(astr)>90|double(astr)<65) astr=input('输入错误，应该全部为大写字母:') end

an=double(astr);lh=length(an); if mod(length(an),2)==1 an=[an,an(length(an))]; end

an=an-ZERO;

for i=1:length(an) if an(i)==26 an(i)=0; end end

c=reshape(an,2,length(an)/2); dn1=mod(enmat*c,m); dn2=mod(demat*c,m);

en1=reshape(dn1,1,length(an)); en2=reshape(dn2,1,length(an)); en1=en1+ZERO; en2=en2+ZERO; for i=1:length(en1) if en1(i)==64 en1(i)=90; end end

for i=1:length(en2)

if en2(i)==64 en2(i)=90; end end

en1=en1(1:lh); en2=en2(1:lh); char(en1) char(en2)

输入要加密的明文文字（全部为大写字母):'MYFRIENDWILLCOMETODAY' astr =MYFRIENDWILLCOMETODAY

输入要解密的密文密码（全部为大写字母):'ROOAJJKSLSWOXSFAWXH' astr =ROOAJJKSLSWOXSFAWXH ans =VSQCDDWEXEASJEHCSTX ans =HEWILLGOHOMETONIGHT

5（15分）、设A0(0,0)为一导弹发射点，发现位于B(0,100)处一架敌机沿水平方向逃离，随即发射一枚导弹予以打击，现已知导弹时刻对准敌机，且速率为飞机速率的两倍（设飞机速度为1）。试编程模拟导弹打击敌机的动态过程，并实时给出飞机和导弹的位置坐标。如果敌机飞行60单位距离之外即逃出我方空域，那么，要想在我方空域内击落敌机，则导弹的速度至少应提高到敌机速度的多少倍？ 解：k=0;

A=[0,0];%导弹初始位置 B=[0,100];%飞机初始位置 v=1; dt=1; d=100; while d>0.5

plot(A(1),A(2),'r');%画导弹位置 hold on

plot(B(1),B(2),'g+');%画飞机位置 pause(0.2); k=k+1; B=B+[v*dt,0]; e=B-A;

d=norm(e);e0=e/d; A=A+2.0*v*dt*e0; fprintf('k=%.0f

d=%.f\\n',k,A(1),A(2),B(1),d); end

k=1 A(0.02,2) B(1.00,100) d=100 k=2 A(0.06,4) B(2.00,100) d=98 k=3 A(0.12,6) B(3.00,100) d=96 k=4 A(0.20,8) B(4.00,100) d=94 k=5 A(0.31,10) B(5.00,100) d=92 k=6 A(0.43,12) B(6.00,100) d=90 ??

k=64 A(61.67,100) B(64.00,100) d=4 k=65 A(63.66,100) B(65.00,100) d=3 k=66 A(65.66,100) B(66.00,100) d=2 k=67 A(67.66,100) B(67.00,100) d=1 k=68 A(69.66,100) B(68.00,100) d=0

120A(%.2f,%.f) B(%.2f,100)

100806040200010203040506070

所以导弹的速度至少提高为敌机的68倍。

（15分）6、设某城市有某种物品的10个需求点，第i个需求点Pi的坐标为(ai,bi),道路网与坐标轴平行，彼此正交。现打算建一个该物品的供应中心，且由于受到城市某些条件的限制，该供应中心只能设在x界于[5,8]，y界于[5,8]的范围之内。问该中心应建在何处为好? Pi点的坐标为：

ai 1 bi 2 4 3 5 9 12 6 4 5 20 17 8 10 8 9 10 8 解：function f = ff15(x) a=[1 4 3 5 9 12 6 20 17 8]; b=[2 10 8 18 1 4 5 10 8 9];

18 1 f(1) = abs(x(1)-a(1))+abs(x(2)-b(1)); f(2) = abs(x(1)-a(2))+abs(x(2)-b(2)); f(3) = abs(x(1)-a(3))+abs(x(2)-b(3)); f(4) = abs(x(1)-a(4))+abs(x(2)-b(4)); f(5) = abs(x(1)-a(5))+abs(x(2)-b(5)); f(6) = abs(x(1)-a(6))+abs(x(2)-b(6)); f(7) = abs(x(1)-a(7))+abs(x(2)-b(7)); f(8) = abs(x(1)-a(8))+abs(x(2)-b(8)); f(9) = abs(x(1)-a(9))+abs(x(2)-b(9)); f(10) = abs(x(1)-a(10))+abs(x(2)-b(10)); x0 = [6; 6]; AA=[-1 0 1 0 0 -1 0 1];

bb=[5;8;5;8];

[x,fval] = fminimax(@ff15,x0,AA,bb)

求解得到中心应建在（8,8）位置最好。

7（25分）、Logistic 模型

N(t)?K/(1?ea?rt)

式中N(t)为 t时刻种群密度，K为环境最大容纳量，r为种群内增长率，a是由 K和种群初始密度决定的常数。在总的增长趋势上很好地描述了种群的增长规律，但有一定局限性，如对于季节性生物来说，其生长过程主要受季节性生态因子周期性变化的影响，但该模型没有考虑生物种群数量季节性的周期性变化。已知棉蚜种群为季节性生物种群，通过实验，有如表1所示的一组实验数据，试建立模型，并拟合试验数据：

表1 棉蚜种群数量变化的实验数据

t 1 2 3 4 5 6 7 8 9 10 11 y 3940 5186 6208 5420 3074 3384 1706 1100 710 316 20 提示：各种生态因子的影响最终都可归结为对环境容纳量的影响，从而影响种群的增长过程。实验要求（1）写出数学模型与实验步骤；（2）画出拟合曲线；（3）进行结果分析。 解：（1）种群的实际增长过程可表示为普适公式：

N（t）=[K(1+f(t)]/（1+e(a-rt)）

影响因数f（t）的具体形式可近似地表示为单一的正弦或余弦函数：

F（t）=h*cos(α+w*t)

季节性生物种群总的增长规律可表示为：

N（t）=[K(1+ h*cos(α+w*t)]/（1+e(a-rt)） 代码：

t=[1 2 3 4 5 6 7 8 9 10 11];

y=[3940 5186 6208 5420 3074 3384 1706 1100 710 316 20]; scatter(t,y,5,'r','filled'),hold on;

ji=inline('b(1).*(1+b(2).*cos(b(3)+b(4).*t))./(1+exp(b(5)-b(6).*t))','b','t');

b=[4200 1 1 0.1 0.1 0.20];

[beta,Res,Re]=lsqcurvefit(ji,b,t,y); parameters=beta,predictions=y+Re ss=sum((y-mean(y)).^2); rs=sum(Re.^2); R=(ss-rs)/ss syms b t;b=beta;

y=subs(b(1).*(1+b(21).*cos(b(3)+b(4).*t))./(1+exp(b(5)-b(6).*t))); ezplot(y,[0,11])

title('棉蚜种群消长实验数据拟合曲线'),xlabel('时间'),ylabel('数量') 运行 M 文件，结果显示为： parameters = 1.0e+003 *

4.2067 0.0010 0.0002 0.0003 0.0013 0.0013

Predictions =

1.0e+003 * 3.8323

Columns 1 through 9

5.5295 5.7515 5.0468 4.0153

2.9312 1.9286 1.0961 0.4995

Columns 10 through 11 0.1846 0.9647 （2）图形显示

0.1755

R =

（3）结果分析

由上述拟合结果可以看出模型很好的描述了季节性生物种群的增长规律，拟合数据与实验数据的相关系数大于0.90。模型中参数h,w分别反映了生态因子对种群增长影响的幅度和角频率，而且h和w可以为时间t的函数，因此对于周期性生态因子的影响，可用正弦或余弦函数来描述。

本文来源：https://www.bwwdw.com/article/1ytv.html