数值建模与仿真实验报告

1、1、启动和退出matlabo在命令窗口认识help、demo命令，并査找cell、floor、fix、round、 rem、sign等函数的用法。2a=l,0,-l;2,4，l;-2,0,5,b=0,-l,0;2,l,3;l,l,2】求 2a+b、a2-3b、a*b、b*a、a.*b、a/b、ab、a./b、a.b. » b=0,-l,0;2, 1,3;1, 1,20-10 213112» a=l,0,-l;2,4, l:-2,0,5a =10-12 41-205» c=2*a+bc =2-1-2695-3112» d=a*a-3*bd =3 3-621

2、3-215-321» 5=a*be =-1_2-293145710?=b*a-2-2_1-444-11410g=a.*b000443-2010h=a/b2. 00003. 00005. 00005. 00003. 0000-4.00007. 0000-9.000016. 0000» i=abi =0.3333-1.33330. 66670.25001. 00000. 25000. 33330.33330. 6667» j=a. /bwarning: divide by zero.j =inf1.0000 -2.000004.00000inf 0. 3333 2.

3、5000» k=a. bwarning: divide by zero.k =01.0000 -0.5000inf0.2500inf03. 0000 0.40003、利用函数产生3x4阶单位矩阵和全部元素都是4.5的4x4阶常数矩阵» l=eye(3,4)10 0 0 0 10 0 0 0 10» =4.5*ones(4)4.50004.50004.50004.50004.50004.50004.50004.50004.50004.50004.50004.50004.50004.50004.50004.50004、利用函数产生5x5阶随机分布的矩阵和5x5阶正态分

4、布的随机矩阵» n=rand(5, 5)n =0. 95010. 76210.61540.40570. 05790. 23110.45650. 79190.93550. 35290. 60680.01850. 92180.91690.81320.48600. 82140. 73820.41030. 00990. 89130.44470.17630. 89360. 1389» 0=randn(5j 5)0 =-0.43261. 1909-0.18670. 11390. 29441.66561. 18920.72581.0668-1.33620.1253-0. 0376-0.5

5、8830. 05930.71430.28770. 32732. 1832-0.09561.6236-1.14650. 1746-0.1364-0.8323-0.69185、练习文件读写操作1、gcd函数用于求两个整数的最大公约数。先用help命令查看该函数的用法，然后利用该 函数求15和35的最大公约数。» :g, c, d】=ccd<15, 352、己知矩阵 a=8,9,5;36,7,11;2v8,5,b=-1,3，-2;2,o,3;-3,1,9求下列表达式: (1)a+5*b 和 a-b(2) a*b 和 a.*b(3) a/b 和 ba(4) aa3 和 a.a3(5)

6、a,b72-633. 134 1 -1.2b01 -1.分8781o. 7 i9g-1. 26s39. 47s61. os544. 8o9r-1. <a«rs3. 0244 3. 7 195 i.31t1627216714 1423342-85619002m4g20o2990&1272040060-343261*51212b1331129分6-23、分别用for和while循环结构编写程序，计算下式的值:2t -io o 33-s 40< oe34-o.<102412. 0244-3.23176. <1612-2.0366» s=0;

7、7; for k=l:100»=(2*k-l)*2;s=s+m;end» s13333004、求下而分数序列前20项之和: 2 3 5 8 t 2 3 5 » s s=0；for n=0:10- _s s=s+l/ga»»a(n):end32.66031、1行100列的fibonacc数组a，a=a(2)=l,a(i)=a(卜l)+a(i-2),川for循环语句来寻求该数组屮第一个大于10000的元素，并指出其位置h要求编写m文件。n=.00;a=oncs(l,n); for i=3:na(i)=a(i-l)+a(i-2); if a(i)&g

8、t;10000 a(i), break;ans =10946end;212、创建一个4x3阶的服从0-1均匀分布的随机矩阵a,求岀a中各行的平均伉b(列句量)， 将b补在a的右方构成4x4阶的矩阵c,并提取c的下三角矩阵。 mrand(4> 3); b=a(:,1)+ag,2)+ag,3> oio.b.0.95010000.23110.7621000.60680.45650.615400.48600.01 於o. r»i»0.4»l f«u,(ch； c»trxl(f)3、根据麦克劳林公式可以得到：el+l+l/2kl/3!+-+l

9、/n!,编写一段程序，求当n=10时e 的近似值。s=0;for n-0:10s=s+l/ganuha(n);s =ends2.71834、编写一个函数，使其能够产生如下的分段函数2x-l,x<0< jca2-1,0<x<1function g=f (x) if x=0 g= (-2)*x-l; elseif x>0ax<l (=x2-l; elseif x>=l 家=x*1/2-1: end1、选择合适的步长绘制出下列函数的图形(1) in 岩，x e (-1,1)(2)yjcosx) ,x g71 71 2'2»丨龙鼴竈on丨麝m

10、f，启泊 clearclcx=*1:0. 01:1: t=：l-x)./(l4x); y=og(t); plrt (x,y/-); ax,s(h.2,1.2,-2,2】) axw square")jitor c:matlab7workuntwed.m edit text cell tools debug desktop多画 x 龜on函典丄_cl»arcl:x=-pi/2:0. ul:pi/2; t=:os(x); y=iqrt(t); plrt gc，y，): axi«(-2,2,0,1.2)3xls( square )view insert tools des

11、ktop window help 及叹 ® 斐 q i 2、在同一坐标下绘制函数x5x 5sin(x) 在0，力的曲线法bs k 81d0 idxb «n caor - cmatlabaworkuntwed.mdit 了ext cell tools debug desktop window helpclearclcx-0:0.01:px;plot(x, (x), - ,x,x. .2, : ,x,-(x). 2, sin(x)j legend(*x /x 2*/ (-x) .2 / sm(x)-l);3、在极坐标系绘制下列函数的曲线3 (1) cos (t)-12(2) 2

12、t +1insert tools desktop window help; x. a m o. 14敦 m a i" qclearclct=0:0.01:2tpi; r=cos(t). 3-1; polar(t, r);:ait i ext ceil ioois ueoug desktop winacb吸叹吠os n210120150240270|龙黾龜m m曾9a aclearclcfe戏叹)®甓ul b t=0:0.01:2*pi; r=2. *t/2+l; polar (t,r);4、绘制二维正态分布密度函/(u)1+y2)的三维图形dit text cell too

13、ls debug desktop windowclearclcx,y=>eshjrid(-2:0. 1:2); z=l/2*pi. »exp(-(x. *2+y. a2)/2); plot3(x,y, z：dit view insert tools desktop window help1、对表达式进行化简ii/command windowclearclcsyms x ff=sqrt(1/x. 3+4/x.-2+6/x+8). a(1/3) simplify (f)(1/x*3+4/x*2+6/x+8)*(1/6) ans(l+4*x+6*x*2+8*x *3)/x*3)a(1

14、/6)2、求表达式的极限lim y/x: - 2x -h 5 ; limx->c.vl + -v2-lx->0command ?vindowclear(xk*x+5* (1/2clcsynis x sartss=sqrt(x. a2-2*x+5) limit (s)5* (1/2)射 command windowclearclcs =sy>s x ss= (sqrt (1+x. a2)-l)./x)(l+x-2)*(l/2)-l)/x limit (s)ans =03 + 2xdx多 x囈亀oclearclcsyms x s msqrt (lx*2) int (s)|comm

15、and windowans =-l/3*(l-xa2)"(3/2)clearclcs =syjns x ss-l/(3+2*x)1/(3+2*x)int (s，x，0, 3) |ans =l/2*log(3)3、求罕只分4、求下列函数的极限问题:(1)= ea,x 分别求当 x 0+时的极限clearclcsyms x ss=l/ (1+exp (-1/x) limit (s,x, 0, right')|1/(1+exp (-1/x)ans =clearclcsyms x ss=1/(1+exp(-1/x) limit (s, x, inf) i1/(1+exp (-1/x)

16、ans =1/25、微分m题求y"并简化结果表达式clcsy»s x a s m nw=log(exp(x)+sqrt(x*2+a"2) n=diff(») s=diff (n) simplify(s)/.1、求下列多项式f (x) =0时的根。(1) f (x) =x3-2x2-5(2) f (x) =x3+2x2+10x-20j *g rl command windowclearclcp=l,-2,0,-5 x=roots(p)p =1-20-52.6906-0.3453 + 1.3187i 0.3453 - 1. 3187icommand wind

17、ow loc(kp(i)+(x-2n-2)#ot)r2(«p(:w(f2ht(l/2w/(«p(xh(i»r(l/2)s =舫：(expfc)n apwaz+apwh'hili+a's'lwl+expwhwi'jl'iwlia'htttpwn'ktexpwntal/tx'zl'q/a/feipw+lx'jhl'lwd'z»多 i x也亀1n dclearclcp= 1, 2, 10, -20x=roots(p)2 10 -201.6844 + 3.4313i

18、1.6844 3.4313i 1. 36882、求函数f (x) =2x2-6在x=-4 3之间的极小值和x=-2附近的零点。command window典iclearx,y=fjninbnd(*2.*x. 2-6',-4,3)x,y=fzeroc 2. *x. a2-6，-2)6 7321-8. 8818e-0163、求下列微分方程在13】区间内的数值解：(1) (2)2xdxdvdxxo) = 2.v(0) = 3edit text £ell tols debug运丨x "b亀n a |鼉|funct ion f=f (x, y) f=2*x+2*x/y;

19、87; b,t)=0de45cf»rtfu鬈colons 1 throu(h 151.(0001.03351.06701.1005i.1m01.1m01.23401.2m01.33401.38401.0101.48401.53401.58401.6340coluuu 16 through 30ij8401.7m01.78401.8m01.8940i.9m01.98402.03402.08402.13402.18402.23402.28402.33402.3840column 31 throu(h 452.<m02侧2.53402.59402.6m02.6m02.73402.fm

20、02.sm02.8m02.9»02.95052. w02.983510000coluvu 1 through 142.(0002.10132.20442.30912.415$2.57802.?<432.91473.08923.26793.45)73.63783.82914.0247coluuu 15 throu(h 284.：2464.42894.63754.85045.05785.28965.51595.74655.98176.22136.46556.71426.96737.2251coluiru 2$、throu(h 427.<8747.75428.02578.301

21、78.w238.86759.15m9.45189.751010.05v10.363110.676210.994011.0999coluvu 43 throu(h 45ll.：06311.313311.4208£dit text £cll tfiols defcug qei篡囈竃m oafunction f=f (x,y) f=l/(5*y);» clear.tt,tl 吻wscfjl,3】,3);»»rcolws 1 throo(h 161.0)001.05001.10001.15001.20001.25001.30001.35001.4000

22、1.45001.50001.55001.60001.65001.7000colians if through 32l.卿u8soo1,90001.95002. oogo2.05002.10002.15002.20002.2500130002.35002.40002.45002.5000colws 33 thiou(h 412.6)00165002.f0002.f5002.80002.85002.90002.95001 0000»ran:colwu 1 thiou(h 16x(d001 0033xooer3.0100x01333.016$10199x02321 02653.0299x0

23、332x03643.0»?xoxm63colws 17 through 32xc529105613.05541062?10659106923.072510757x0z901 08221 08541 088f3.09191 09521 0984cohans 33 throu(h 41xbm811081x 11133.1145iu713.120911241x12z3113051、已知多项式 pl(x)=3x+2,p2(x)=5x2-x+2,p3(x>=x2-0.5，求：(1) p(x)=pl(x)p2(x)p3(x)；1.75002.5500xw963.1016clearclcpl

24、=3, 2; p2=5,-l, 2: p3=l, 0,-0. 5; p4=conv(pl, p2) p5=conv(p3, p4)15.00007.0000-3.50000.5000-2.0000-2.0000ans =0. 70710. 1000 + 0. 6245i 0. 1000 - 0. 6245i -0. 7071 -0.6667(2) p(x)=o的全部根。ans =2/5clearclcsyns x ss=sin(2*x)./sin(5*x): linit(s, x, 0)2、求极限值：(1) lim sin 2%(2) lim(l + )2sin 5xxclearclcarts

25、 =syms x sexps=(l+l./x)' (2*x); limit (s, x, inf)2*x/(l+it2)158z1<9y=dsolvec dy= (y"2-x*y)/x 2', y(l= l，x3、求微分方程心娜=人刈=1的特解。 程序：y=dsolve('dy-(ya2-x*y)/(xa2) * r，y(l)=l，'x')结果： y =2*x/(l+xa2)4、解线性方程组3x + 4 y - 2 z = 12 、45x+5y+4z=23 6x + 2 y - 3 z = 4syjns x y zsl=3*x+4*y-2

26、*z-12; s2=45*x+5*y+4*z-23;s3=6*x+2*y-3*z-4 ;|x, y, z=solve(sl, s2, s3)实验3人口预测与数据拟合一、实验目的：通过对人口预测问题的分析求解，了解利用最小二乘法进行数据拟合的基本想，熟悉 寻找最佳拟合曲线的方法，掌握建立人口增长数学模型的思想方法。二、实验器材和环境matlab2014 版本，windows7 系统 malthus 模型、logistic 模型三、实验内容和步骤实验问题：1981-2016年各年我国人口数的统计数据如下表所示（单位：亿）：19811982198319841985198619871988198919

27、9010.00710.16510.30110.43610.58510.75110.93011.10311.2711.433199119921993199419951996199719981999200011.58211.71711.85211.98512.11212.2391236312.47612.57912.674200120022003200420052006200720082009201012.76312.84512.92312.99913.07613.14513.21313.28013.34513.40920112012201320142015201613.47413.54013.6

28、0713.67813.74613.827根据上述数据，建立我国人口增长的近似曲线，并预测2020年、2025年、2030年我国 的人口数量。1.最小二乘法程序如下：clearclcx0=l:36;yo=10.007 10.165 10.301 10.436 10.58510.751 10.930 11.103 11.27 11.433.11.582 11.717 11.852 11.985 12.112.12.239 12.363 12.476 12.579 12.674.12.763 12.845 12.923 12.999 13.076.13.145 13.213 13.280 13.34

29、5 13.409.13.474 13.540 13.607 13.678 13.746.13.827;a=polyf it (xo, yoz 1)y2020=polyval (a，40)y2025=polyval (a, 45)y2o3o=polyval (a, 50)求得：a =0.1081 10.2622y2020 =14.5860y2025 =15.1264y2o3o =15.66692. malthus 模型clearclct=1:36;x(t)=10.007 10.165 10.301 10.436 10.58510.75110.93011.10311.2711.433.11.582

30、11.71711.85211.98512.112.12.23912.36312.47612.57912.67412.76312.84512.92312.99913.07613.14513.21313.28013.34513.409.13.47413.54013.60713.67813.746.13.827 ; y=log (x (t); a=polyfit (t,y,1) r=a (1), x0=exp (a (2) xi=x0 *exp(r*t); plot(t,x(t), r，t,xl, b)人口预测程序：clearclct = 40;%t是变量，此时t对应2020年r=0.0093x0=

31、10.291x(t)=x0*exp (r*t)求得：x(40)=14.9285malthus 模型首先计算参数增长率人口最大值n程序如下.clearclct=18;x0=10.007;xl=12.476;x2=13.827;r=(l/t)*log(1/xo-l/xl)/(l/xl-l/x2)n= (1-exp(-r*t)/(1/xl-(1/xo)*exp(-r*t) r =0.0515n =14.8838 然后进行曲线拟合 程序如下： clear clct=l:36;x(t)=10-007 10.165 10.301 10.436 10.585.10.75110.93011.10311.271

32、1.43311.58211.71711.85211.98512.11212.23912.36312.47612.57912.67412.76312.84512.92312.99913.07613.145 13.213 13.280 13.345 13.409.13.474 13.540 13.607 13.678 13.746.13.827 ;y(t)=14.8838./(1+(14.8838/10.007-1)*exp(-0.0515*t); plot(t,x (t), 'o:1,t,y (t), r- )拟合阁形如下:最后根据malthus模型进行人口预测程序如下：clearclc

33、t=40 %t是变量，此时t对应2020年y (t)=14.8838./ (1+(14.8838/10.007-1)*exp(-0.0515*t) y(40) =14.0134 y(45) =14.2019 y(50) =14.3512模型一：用malthus模型的基本假设是：假设人口的增长率为常数r，记每年时刻t的人口为x(t), 即x(t)为模型的状态变量，且初始时刻的人口为xo,于是得到如下函数： dx/dt=r*xx (0)=x0根据假设可得函数如下.yl=xo*exp(r*x)编辑程序如下：clearclcxdate=1981:1:2016;tl=1981:l:2020;ydate=

34、10.007 10.165 10.301 10.436 10.58510.751 10.930 11.103 11.27 11.433.11.582 11.717 11.852 11.985 12.112.12.239 12.363 12.476 12.579 12.674.12.763 12.845 12.923 12.999 13.076.13.145 13.213 13.280 13.345 13.409.13.474 13.540 13.607 13.678 13.746.13.827;p=polyfit(xdate,log(ydate),1)运行结果：p =0.0090-15.481

35、3将此函数转化为线性关系为in (yl) =ln (x0) +r*x与y=at+b对应，则利用线性拟合即可求解, 过程如下.1.将x与y的数据先进行线性拟合，由结果可知a=0.0147, b=-26.7783,则微分方程的解析 式为：y=exp(-154 813)*exp(00090*x);当 x=2020 时，y=14.8604 ;当 x=2025 时，y=15.5444 ;当 x=2030 时，y=16.2599 ;下面是函数拟合前后的过程编辑程序：clearclcsyms x y x2 y2 x3 y3 x=1981:l:2016;y=10.00710.16510.30110.43610

36、.585.10.75110.93011.10311.2711.43311.58211.71711.85211.98512.112.12.23912.36312.47612.57912.674.12.76312.84512.92312.99913.07613.14513.21313.28013.34513.40913.47413.54013.60713.67813.74613.827 ;p, s=polyfit(x,log (y),1) x2=2020:5:2030;y2=exp(-15.4813)*exp(0.0090*x2); x3=1981:l:2030;y3=exp(-15.4813)*exp(0.0090*x3); plot(x,y) hold onplot(x2,y2z '*) hold onplot (x3,y3, 1) hold on 拟合曲线如下图17.,模型二：clearclct=1981:l:2016;tl=1981:l:2030;x=10.00710.16510.30110.43