数学实验1
- 格式:doc
- 大小:285.50 KB
- 文档页数:22
MATLAB数学实验报告1Matlab数学实验报告⼀、实验⽬的通过以下四组实验,熟悉MATLAB的编程技巧,学会运⽤MATLAB的⼀些主要功能、命令,通过建⽴数学模型解决理论或实际问题。
了解诸如分岔、混沌等概念、学会建⽴Malthu模型和Logistic 模型、懂得最⼩⼆乘法、线性规划等基本思想。
⼆、实验内容2.1实验题⽬⼀2.1.1实验问题Feigenbaum曾对超越函数y=λsin(πx)(λ为⾮负实数)进⾏了分岔与混沌的研究,试进⾏迭代格式x k+1=λsin(πx k),做出相应的Feigenbaum图2.1.2程序设计clear;clf;axis([0,4,0,4]);hold onfor r=0:0.3:3.9x=[0.1];for i=2:150x(i)=r*sin(3.14*x(i-1));endpause(0.5)for i=101:150plot(r,x(i),'k.');endtext(r-0.1,max(x(101:150))+0.05,['\it{r}=',num2str(r)]) end加密迭代后clear;clf;axis([0,4,0,4]);hold onfor r=0:0.005:3.9x=[0.1];for i=2:150x(i)=r*sin(3.14*x(i-1));endpause(0.1)for i=101:150plot(r,x(i),'k.');endend运⾏后得到Feigenbaum图2.2实验题⽬⼆2.2.1实验问题某农夫有⼀个半径10⽶的圆形⽜栏,长满了草。
他要将⼀头⽜拴在⽜栏边界的桩栏上,但只让⽜吃到⼀半草,问拴⽜⿐⼦的绳⼦应为多长?2.2.2问题分析如图所⽰,E为圆ABD的圆⼼,AB为拴⽜的绳⼦,圆ABD为草场,区域ABCD为⽜能到达的区域。
问题要求区域ABCD等于圆ABC的⼀半,可以设BC等于x,只要求出∠a和∠b就能求出所求⾯积。
(1)参数方程:z=2^2^/2^2^sin y x y x ++(-8<=x<=8,-8<=y<=8) (2)程序:[X,Y]=meshgrid(-8::8);r=sqrt(x.^2+y.^2)+eps;Z=sin(r)./r;Mesh(x,y,z)Axis square(3)程序的输出结果:3:球面,椭球面,双叶双曲面,单叶双曲面1球面: (4):参数方程:⎪⎩⎪⎨⎧===ϕθϕθϕcos *sin *sin *cos *sin *R z R y R x 0π<=θ<2* 0<=ϕ<π (5)程序:u=[0:pi/60:2*pi];v=[0:pi/60:pi];[U,V]=meshgrid(u,v);R=3;X=R*sin(v).*cos(u);Y=R*sin(v).*sin(u);Z=R*cos(v);Surf(x,y,z);axis equal;(3)程序输出结果:2椭球面: (1)参数方程:⎪⎩⎪⎨⎧===ϕθϕθϕcos *sin *sin *cos *sin *c z b y a x 0<=θ<2*π 0<=ϕ<=π (2)程序:ezsurf(‘3*sin(u)*cos(v) ,’3*sin(u)*sin(v)’,’1*cos(u)’,[0,pi,0,2*pi]);(3)程序的输出结果:3单叶双曲面: (1)参数方程:⎪⎩⎪⎨⎧===ϕθϕθϕtan sin *sec *cos *sec *z a y a x 0<=θ<2*π -π/2<ϕ<π/2 (2)程序:ezsurf(‘3*sec(u)*cos(v),’3*sec(u)*sin(v)’,’5*tan(u)’,[-pi/2,pi/2,0,2*pi]);axis auto(3)输出程序结果:4双叶双曲面: (1)参数方程:⎪⎩⎪⎨⎧===ϕθϕθϕsec *sin *tan *cos *tan *c z b y a x 0<=θ<2*π -π<ϕ<3*π/2,ϕ≠π/2(2)程序:ezsurf(‘3*tan(u)*cos(v)’,’3*tan(u)*sin(v)’,’5*sec(u)’,[-pi/2,3*pi/2,0,2*pi]);axis auto(4) (3)输出程序结果:抛物螺线: (1)参数方程:⎪⎩⎪⎨⎧===2^*sin **cos **t c z t t b y t t a x 0<T<+∞ (2)程序:ezplot3(‘2*t*cos(t)’,’2*t*sin(t)’,’t.^2/3’,[0,50]);(3)输出程序结果:(5)马鞍面: (1)参数方程:z=x^2/9-y^2/4 (-25<=x<=25,-25<=y<=25)(2)程序:[X,Y]=meshgrid(-25:1:25);Z=X.^2/9-Y.^2/4;Surf(X,Y,Z)Title(‘马鞍面’)grid off(3)输出程序结果:(6)黎曼函数:(1)程序:n=100;x=[];y=[];k=1;for q=2:nfor p=1:q-1if gcd(q,p)==1 %利用函数gcd(m,n)可求m和n的最大公约数x(k)=p/q;y(k)=1/q;k=k+1;endendendplot(x,y,’.b’); axis([0,1,0,1])(2)程序输出结果:。
实验一 常微分方程1. 分别用Euler 法和ode45解下列常微分方程并与解析解比较: (1) ,(0)1,13y x y y x '=+=<<Euler 法:function [t,y]=euler(Fun,tspan,y0,h) t=tspan(1):h:tspan(2); y(1)=y0;for i=1:length(t)-1y(i+1)=y(i)+h.*feval(Fun,t(i),y(i)); end t=t'; y=y';function f=Fun(x,y) % 常微分方程的右端函数 f=x+y;>> [x,y]=euler('Fun',[0,3],1,0.1)>> [x,y] ans =0 1.0000 0.1000 1.1000 0.2000 1.2200 0.3000 1.3620 0.4000 1.5282 0.5000 1.7210 0.6000 1.9431 0.7000 2.1974 0.8000 2.4872 0.9000 2.8159 1.0000 3.1875 1.1000 3.6062 1.2000 4.0769 1.3000 4.6045 1.4000 5.1950 1.5000 5.8545 1.6000 6.5899 1.7000 7.4089 1.8000 8.3198 1.9000 9.3318 2.0000 10.4550 2.1000 11.7005 2.2000 13.0805 2.3000 14.6086 2.4000 16.2995 2.5000 18.1694 2.6000 20.2364 2.7000 22.5200 2.8000 25.0420 2.9000 27.8262 3.0000 30.8988ode45:>> [x,y]=ode45('Fun',[0,3],1) ans =0 1.0000 0.0502 1.0528 0.1005 1.1109 0.1507 1.17460.2010 1.2442 0.2760 1.3596 0.3510 1.4899 0.4260 1.63610.5010 1.7996 0.5760 1.9817 0.6510 2.1838 0.7260 2.4074实验一 常微分方程0.8010 2.6544 0.8760 2.9264 0.9510 3.2254 1.0260 3.55351.1010 3.9131 1.1760 4.3065 1.2510 4.7364 1.3260 5.20561.4010 5.7172 1.4760 6.2744 1.5510 6.8810 1.6260 7.54061.7010 8.2574 1.7760 9.0359 1.8510 9.8808 1.9260 10.79742.0010 11.7912 2.0760 12.8683 2.1510 14.0351 2.2260 15.29862.3010 16.6664 2.3760 18.1466 2.4510 19.7478 2.5260 21.47962.6010 23.3522 2.6760 25.3764 2.7510 27.5641 2.8260 29.92812.9010 32.4820 2.9257 33.3694 2.9505 34.2796 2.9752 35.21343.0000 36.1711解析解:>> y=dsolve('Dy=x+y','y(0)=1','x') y =2*exp(x) - x - 1(2) 20.01()2sin(),(0)0,(0)1,05y y y t y y t ''''-+===<< Euler 法:实验一常微分方程function f=Fun(t,y)% 常微分方程的右端函数f=[y(2);0.01*y(2)^2-2*y(1)+sin(t)];>> [t,y]=euler('Fun',[0,5],[0,1],0.2)ode45:>> [t,y]=ode45('Fun',[0,5],[0,1])t =0 0.0001 0.0001 0.0002 0.0002 0.0005 0.0007 0.0010 0.0012 0.00250.0037 0.0050 0.0062 0.0125 0.0188 0.0251 0.0313 0.0627 0.0941 0.12550.1569 0.2819 0.4069 0.5319 0.6569 0.7819 0.9069 1.0319 1.1569 1.28191.4069 1.5319 1.6569 1.7819 1.90692.0319 2.1569 2.2819 2.4069 2.53192.6569 2.7819 2.90693.0319 3.1569 3.2819 3.4069 3.5319 3.6569 3.78193.90694.0319 4.1569 4.2819 4.4069 4.5319 4.6569 4.7427 4.8285 4.91425.0000y =0 1.0000 0.0001 1.0000 0.0001 1.0000 0.0002 1.0000 0.0002 1.00000.0005 1.0000 0.0007 1.0000 0.0010 1.0000 0.0012 1.0000 0.0025 1.00000.0037 1.0000 0.0050 1.0000 0.0062 1.0000 0.0125 1.0000 0.0188 1.00000.0251 0.9999 0.0313 0.9998 0.0627 0.9987 0.0941 0.9965 0.1253 0.99340.1564 0.9893 0.2786 0.9632 0.3966 0.9220 0.5085 0.8662 0.6126 0.79670.7072 0.7146 0.7908 0.6210 0.8620 0.5176 0.9198 0.4058 0.9632 0.28760.9915 0.1647 1.0043 0.0392 1.0013 -0.0869 0.9826 -0.2117 0.9485 -0.33310.8996 -0.4490 0.8365 -0.5578 0.7605 -0.6577 0.6725 -0.7471 0.5742 -0.8246实验一 常微分方程0.4669 -0.8889 0.3525 -0.9393 0.2327 -0.9748 0.1095 -0.9950 -0.0154 -0.9996-0.1398 -0.9887 -0.2619 -0.9624 -0.3798 -0.9212 -0.4916 -0.8657 -0.5957 -0.7970-0.6904 -0.7161 -0.7742 -0.6242 -0.8460 -0.5228 -0.9046 -0.4134 -0.9491 -0.2978-0.9789 -0.1777 -0.9934 -0.0549 -0.9945 0.0300 -0.9883 0.1146 -0.9748 0.1985-0.9543 0.28092. 求一通过原点的曲线,它在(,)x y 处的切线斜率等于22,0 1.57.x y x +<<若x 上限增为1.58,1.60会发生什么?function f=Fun(x,y) % 常微分方程的右端函数 f=2*x+y.^2;>> [x,y]=ode45('Fun',[0,1.57],0) x =0 0.0393 0.0785 0.1178 0.1570 0.1963 0.2355 0.2748 0.3140 0.3533 0.3925 0.4318 0.4710 0.5103 0.5495 0.5888 0.6280 0.6673 0.7065 0.7458 0.7850 0.8243 0.8635 0.9028 0.9420 0.9813 1.0205 1.0598 1.0990 1.1383 1.1775 1.2168 1.2560 1.2953 1.3345 1.3738 1.4130 1.4248 1.4367 1.4485 1.4604 1.4722 1.4840 1.4959 1.5077 1.5140 1.5203 1.5265 1.5328 1.5376 1.5424 1.5472 1.5519 1.5543 1.5567 1.5591 1.5614 1.5631 1.5647 1.5664 1.5681 1.5685 1.5690 1.5695 1.5700 y =实验一 常微分方程0 0.0015 0.0062 0.0139 0.0247 0.0386 0.0556 0.0758 0.09920.1259 0.1559 0.1895 0.2266 0.2675 0.3124 0.3615 0.4152 0.4738 0.5378 0.6076 0.6841 0.7679 0.8601 0.9620 1.0751 1.2014 1.3434 1.5045 1.6892 1.9037 2.1557 2.4577 2.8282 3.3003 3.9056 4.7317 5.9549 6.4431 7.0116 7.6832 8.4902 9.4821 10.7170 12.3090 14.4551 15.9220 17.7080 19.9390 22.8164 25.6450 29.2282 33.9673 40.5910 44.9434 50.3088 57.1229 66.1087 74.3108 84.7123 98.4901 117.7875 124.9206 132.9699 142.1268 152.641500.20.40.60.81 1.2 1.4 1.6若x 上限增为1.58,1.60,则超出运算的范围,发生溢出。
数学实验1重庆大学学生实验报告实验课程名称数学实验开课实验室DS1422学院计算机学院年级2013专业班物联网2班学生姓名卢锦晔学号20135195开课时间2014 至2015 学年第2学期总成绩教师签名数学与统计学院制开课学院、实验室:DS1422实验时间:2015 年3月30 日课程名称数学实验实验项目名称MATLAB软件入门实验项目类型验证演示综合设计其他指导教师肖剑成绩实验目的[1] 熟悉MATLAB 软件的用户环境;[2] 了解MATLAB 软件的一般目的命令;[3] 掌握MATLAB 数组操作与运算函数;[4] 掌握MATLAB 软件的基本绘图命令;[5] 掌握MATLAB 语言的几种循环、条件和开关选择结构。
通过该实验的学习,使学生能灵活应用MATLAB 软件解决一些简单问题,能借助MATLAB 软件的绘图功能,对函数的特性进行探讨,广泛联想,大胆猜想,发现进而证实其中的规律。
实验内容1.MATLAB 软件的数组操作及运算练习;2.直接使用MATLAB 软件进行作图练习;3.用MATLAB 语言编写命令M-文件和函数M-文件。
基础实验一、问题重述1.设有分块矩阵⎥⎦⎤⎢⎣⎡=⨯⨯⨯⨯22322333S O R E A ,其中E,R,O,S 分别为单位阵、随机阵、零阵和对角阵,试通过数值计算验证⎥⎦⎤⎢⎣⎡+=22S 0RS R E A 。
2.某零售店有9种商品的单件进价(元)、售价(元)及一周的销量如表1.1,问哪种商品的利润最大,哪种商品的利润最小;按收入由小到大,列出所有商品及其收入;求这一周该10种商品的总收入和总利润。
表1.1货号1 2 3 4 5 6 7 8 9 单件进价7.15 8.25 3.20 10.30 6.68 12.03 16.85 17.51 9.30 单件售价11.10 15.00 6.00 16.25 9.90 18.25 20.80 24.15 15.50 销量 568 1205 753 580 395 2104 1538 810 6943.建立一个命令M-文件:求所有的“水仙花数”,所谓“水仙花数”是指一个三位数,其各位数字的立方和等于该数本身。
例如,153是一个水仙花数,因为153=13+53+33。
4.编写函数M-文件sq.m :用迭代法求a =x 的值。
求平方根的迭代公式为)a (211n n n x x x +=+迭代的终止条件为前后两次求出的x 的差的绝对值小于10-5。
5. 近景图 将x 的取值范围局限于较小的区间内可以画出函数的近景图,用于显示函数的局部特性。
局部放大 在绘图时,把x 的范围逐渐缩小,可把函数的细节部分展现的很清楚.特别是观察极限问题时,这种方法比较便利.远景图 函数的远景图,是把x 的范围取得比较大,使我们能够在大范围内观察函数图像.当研究x 趋向于∞时,这种方法给我们带来方便.1)比较函数 33)(,)(,)(x x h x x x g x x f =+== 在x →0时函数的性态。
观察到什么现象?从观察到的现象,反映了什么结论。
2)在日常生活中我们有这样的经验:与幂函数相比,指数函数是急脾气,对数函数是慢性子。
这就是说,当x →∞时,再小的指数函数也比幂函数变化快,再大的对数函数也比幂函数变化慢。
当x →∞时,比较10x y =与 x y 1.1= 的大小.当x →∞时,比较 001.0x y =与 x y lg 1000= 的大小.3)在同一个坐标下作出y 1=e x ,y 2=1+x,y 3=1+x+(1/2)x 2,y 4=1+x+(1/2)x 2+(1/6)x 3这四条曲线的图形,要求在图上加各种标注,观察到什么现象?发现有什么规律?同时用subplot 分别在不同的坐标系下作出这四条曲线,为每幅图形加上标题。
4).作出下列曲面的3维图形,)sin(22y x z +π=;5). 作出函数y=x 4-4x 3+3x+5 (x ∈[0,6])的图形,用小红点标出其在[0,6]之间的最小值点,并在最小值点附近标出该最小值点的坐标值;二、实验过程(一般应包括实验原理或问题分析,变量说明、程序、调试情况记录、图表等,实验结果及分析)1.➢源程序:E=eye(3,3)R=rand(3,2)O=zeros(2,3)a=[3 2]S=diag(a)A=[E R;O S]B=[E R+R*S;O S^2]IfB==A^2disp('right')elsedisp('wrong’)end➢运行结果:E =1 0 00 1 00 0 1R =0.2785 0.9649 0.5469 0.1576 0.9575 0.9706O =0 0 0 0 0 0a =3 2S =3 00 2A =1.0000 0 0 0.27850.96490 1.0000 0 0.54690.15760 0 1.0000 0.95750.97060 0 0 3.0000 00 0 0 02.0000B =1.0000 0 0 1.11402.89470 1.0000 0 2.18750.47280 0 1.0000 3.83002.91180 0 0 9.0000 00 0 0 04.0000Right分析:在该题中应该注意’*’与’.*’的区别,避免不必要的错误。
题中涉及的矩阵函数有eye、Rand、zeros、diag等,还有矩阵合并和结果输出的有关知识。
2.➢源程序:a=1:9b=[7.15 8.25 3.20 10.30 6.68 12.03 16.85 17.51 9.30]c=[11.10 15.00 6.00 16.25 9.90 18.25 20.80 24.15 15.50]d=[568 1205 753 580 395 2104 1538 810 694]x=c.*dprofit=(c-b).*df=sum(x)g=sum(profit)[income,number]=sort(x)➢运行结果:➢分析:此题最重要的部分是sort函数的应用,实现数组元素的排序和输出其在数组中对应的位置。
另外,从中我还学会了利用数组对一系列的数据进行运算处理。
3.➢计算程序for a=1:1:9for b=0:1:9for c=0:1:9A=a*100+b*10+c;if A==a^3+b^3+c^3disp(A)endendendend➢最终结果:153 370 371 407➢分析:此题使我进一步熟悉了循环语句for、while 等的基本格式和应用方法,值得注意的是每一个循环语句或条件语句的出现,必定要对应的出现一个end,这一点很容易忽视。
4.➢源程序:function f=sq(a)if a>=0x=eps;y=1/2*(x+a/x);while (abs(x-y)>=10^(-5))x=yy=1/2*(x+a/x)endf=xelsedisp('their exist errors')end➢运行结果调用sq(16),输出结果是4.0000➢分析:此题要求掌握利用循环语句进行迭代,实现开方方程的求解,另外x的取值有一定的讲究,在该题题设的情况下,很显然x>=0,但从迭代公式的形式可以得到,x不能等于0,这就要求x从一个较为接近0的数开始取,于是就得引用eps(计算机能辨别的最小值)。
5.(1)➢源程序:x=-1:0.0001:1;y1=x;y2=x.^3;y3=y1+y2;plot(x,y1,x,y2,x,y3)➢运行结果分析当x→0时,f(x)与g(x)很接近,而h(x)与前两个函数都不接近。
(2)➢源程序:x=linspace(5000,8000,500);y1=x.^10;y2=1.1.^x;Subplot(1,2,1),plot(x,y1),xlabel('x');ylabe l('y)'); grid;title('y=x^1^0');Subplot(1,2,2),plot(x,y2),xlabel('x');ylabe l('y)'); grid;title('y=1.1^x');➢运行结果分析:从上图可以看出来指数函数变化快➢源程序:x=linspace(5000,8000,500);y1=x.^0.001;y2=1000.*log(x);Subplot(1,2,1),plot(x,y1),xlabel('x');ylabe l('y)'); grid;title('y=x^0.001');Subplot(1,2,2),plot(x,y2),xlabel('x');ylabe l('y)'); grid;title('y=1000.*log(x)');➢运行结果分析:由以上函数图形可知对数函数变化比幂函数慢。
(3)➢源程序:x=linspace(0,2.50);y1=exp(x);y2=1+x;y3=1+x+0.5.*x.^2;y4=1+x+0.5.*x.^2+1./6.*x.^3;plot(x,y1,'b.'),gtext('y1=exp(x)');hold on, plot(x,y2,'y-'),gtext('y2=1+x');plot(x,y3,'g:'),gtext('y3=1+x+0.5.*x.^2');plot(x,y4,'m--'),gtext('y4=1+x+0.5.*x^2+1./ 6.*x.^3');hold off➢作图结果(4)➢作图程序x=-10:0.6:10;y=x;[X,Y]=meshgrid(x,y);Z=sin(pi.*sqrt(X.^2+Y.^2));surf(X,Y,Z)xlabel('x')ylabel('y')zlabel('z')title('z=sin(pi*sqrt(x.^2+y.^2)')shading flat➢作图结果➢分析:在该题中,需要用到三维作图相关函数meshgrid、mesh、surf,注释函数xlabel、ylabel、zlabel、title和开方函数sqrt,最后还涉及到了幂函数等。
通过这道题,我充分运用并巩固了作图方面的基础知识。
(5)➢作图程序x=0:0.1:6;y=x.^4-4.*x.^3+3.*x+5;plot(x,y)[y1,x1]=min(y);hold onplot(x(x1),y1,'r.','MarkerSize',20)hold offa=['x=',num2str(x(x1))];b=['y=',num2str(y1)];min=char(a,b);text(x(x1),y1+50,min)➢作图结果6.自由发挥:自己提出问题,实验探索,广泛联想,发现规律,大胆猜想。