MATLAB在物理中的应用

MATLAB在物理中的应用 0 序 科学昌盛的今天，计算机应用无所不在。利用计算机技术来研究和分析物理学，对有大量的数学分析和计算的学科来说，尤为重要。这里推出的MATLAB程序，是专用于科学和工程计算的程序语言，是目前世界上应用最为广泛的计算语言之一。通过学习，你会觉得编程原来也是轻而易举；物理计算不再会觉得深奥复杂；物理图像更是信手沾来。 在国外有许多理工科类的籍和教材已将MATLAB作为专用的科学计算语言融入到专业内容之中，并得到广泛认可。MATLAB程序具有以下几方面特点： 1.MATLAB语言与人们的思维方式和书写习惯相适应，操作简单、方便。 2.可以方便地将运算结果用图形、图象、声音、动画等形象地表达出来。 3.无需编译，键入命令即可解释运行。 4.能自动选择合适的坐标范围。 5.有功能强大的工具箱。 6.有算法先进的数值计算和符号计算功能。 7.数据类型是矩阵，用户不必定义变量和数据类型，且矩阵大小可任意改变。 本讲座通过MATLAB在大学物理中的应用实例讲解，将大家引进MATLAB之门。旨在帮助同学在学习物理的同时也能掌握一门计算机技术，并能熟练应用计算机技术来学习物理和分析物理。为从简单到复杂逐渐将MATLAB语言的应用介绍给同学，本讲座的入门篇安排在振动与波的章节中，因此，同学从该章节开始学习。 衷心地希望你在本讲座的引导下，能够喜欢和使用MATLAB语言，更希望你能通过MATLAB来分析物理和研究物理。 1 简谐振动 1.1 讲授内容 振动是自然界中普遍存在的一种运动形式。例如，钟摆的运动、气缸中活塞的往复运动、心脏的跳动等。任一物理量（位移、速度、电场强度、电压等）在一定值附近作周期性的变化都可以称为振动。振动的运动形式很多，而且比较复杂，最简单、最基本的是简谐振动，其它任何复杂振动都可以分解为若干简谐振动的叠加。 简谐振动是指物体相对于平衡点的位置随时间按余弦（或正弦）规律变化的运动，其数学表达式为 作简谐振动的系统，称为谐振子。弹簧振子、单摆、复摆等都是谐振子。物体作简谐振动时,速度为 加速度为 1.2 程序编写 无论你有没有编过程序，都将觉得要编写一个物理学方面的程序，会是一件很费力的事，但选择MATLAB可以让你惊奇，在不到一分钟的时间，极少的语句你就能完成用其它语言数倍时间、数倍语句的程序。够酷吧。 还等什么，赶快在你的电脑上安装MATLAB6.1程序吧，相信你安装不会有问。下面，我们要做一个简谐振动的例题，先看题目要求。 例一、已知振动方程,求t=8秒时，位移、速度及加速度各位多少？绘制振动曲线和速度、加速度曲线。 题目有两个小题，先看第一问，计算位移、速度及加速度。 准备好了吗，现在请你跟我一起做： 第一步：启动matlab程序，进入程序主窗口； 第二步：在命令窗口中输入下列程序代码： >>t=8;%赋值 >>y=2*cos(3*t+pi/5);%确定y、t的涵数关系 >>v=-2*3*sin(3*t+pi/5);%确定v、t的涵数关系 >>a=-2*3*3*cos(3*t+pi/5)%确定a、t的涵数关系 第三步：按Enter,可得运算结果。 感觉怎么样？就是这么简单。 下面解决题目的第二问。在命令窗口中直接输入下列程序代码，绘制振动曲线： >>clc;%清除指令窗 >>t=0:0.001:10;%确定时间变量的变化范围及步长 >>y=2*cos(3*t+pi/5);%确定y、t的涵数关系 >>v=-2*3*sin(3*t+pi/5);%确定v、t的涵数关系，pi表示圆周率 >>a=-2*3*3*cos(3*t+pi/5);%确定a、t的涵数关系 >>plot(t,y,'-',t,v,'-.',t,a,':')%绘制图形 %说明：单引号内的符号表示所画曲线的线型，如-（实线）、-.（点划线）、：（虚线）、——（破折线）、p（五角星线）、h（六角形行）等。 按Enter,可得运行结果。 相信现在的你对MATLAB已有了感觉... 2 振动的合成 2.1 讲授内容 实际生活中，一个系统往往会同时参与两个或更多的振动。例如悬挂在颠簸船舱中的钟摆，两列声波同时传入人耳等。一般的振动合成显然是比较复杂，下面仅讨论几种间单情况的简谐振动合成。 一、同方向同频率简谐振动的合成 若两个同方向的简谐振动，频率都是，它们的运动方程分别为 因振动是同方向的，所以这两个谐振动在任意时刻的和位移应在同一直线上，且等于这两个振动位移的代数和，即 合位移仍为简谐振动 二、两个同方向不同频率简谐振动的合成拍 如果两个简谐振动的振动方向相同而频率不同，那么合成后的振动仍与原振动方向相同但不再是简谐振动。现设两简谐振动的振幅都为A，初相位为零，它们的振动方程分别为 合成振动方程为 若两个分振动的频率都较大且其差很小时，即，合振动可看作为振幅随时间缓慢变化的近似谐振动，振幅随时间变化且具有周期性，表现出振动或强或弱的现象，称拍，变化的频率称拍频，变化的振幅为 变化的频率为 三、相互垂直的简谐振动的合成李萨如图 如果两个简谐振动分别在x轴和y轴上进行，他们的振动方程分别为 合成后，可得质点的轨迹为椭圆方程 若两分振动有不同的频率，且两频率之比为有理数时，则合成后的质点运动具有稳定、封闭的轨迹。称其为李萨如图形。 2.2 程序编写 我们已经在第一讲中体验了matlab的编程，可是你一定会生出这样的问号，辛辛苦苦在命令窗口写的一大堆代码怎么不保留？不用担心，matlab程序和其他编程工具一样，也有专门的文件，称m文件，文件名形式为“文件名.m”。你可以用matlab自带的编辑器来输入你的程序代码，当然你也可以用其它编辑器或最经济的文本编辑器，不过别忘记添加文件名的后缀“.m”。 下面，请跟我一起用m文件编辑器来编写matlab程序。 例题：两个振动方向相同而频率不同的简谐振动方程分别为 ， 合成后的方程是 ，请用matlab程序描述合成波和拍频现象。 编程： 第一步：点击matlab图标，打开程序窗口。 第二步：选file—new—m-file,打开编辑器。 第三步：在编辑器中输入下列程序代码： t=1:0.001:10;%设定时间的范围从1到10 a1=input('振幅1=');w1=input('频率1=');%要求用户输入简谐振动1的状态参量 a2=input('振幅2=');w2=input('频率2=');%要求用户输入简谐振动2的状态参量 y1=a1*cos(w1*t+phi1);%简谐振动1的运动方程 y2=a2*cos(w2*t+phi2);%简谐振动2的运动方程 y=y1+y2;%简谐振动1、2的合成谐振3 subplot(3,1,1),comet(t,y1),ylabel('y1')%绘制谐振1、2、3的振动曲线 subplot(3,1,2),comet(t,y2),ylabel('y2') subplot(3,1,3),comet(t,y),ylabel('y'),xlabel('t') sound(y1);pause(2), sound(y2);pause(2),sound(y) 注意：拍频现象要求两个振动的频率都很大，且相差很小。 第四步：按工具栏上的“保存”按钮将程序保存在当前目录下（可以重命名），再按F5运行此程序进入命令窗口根据窗口上的提示输入数据，每次输完一个数据按"Enter"结束， 据所输入的数据（在命令窗口中） 4 阻尼振动 4.1 讲授内容 简谐振动是一种无阻尼振动，而实际上，任何振动物体都要受到阻力的作用，这种振动叫阻尼振动。阻尼振动有摩擦阻尼和辐射阻尼两种。实验指出，当物体以不太大的速率在粘性介质中运动时，介质对物体的阻力与物体的运动速率成正比，方向与运动方向相反，即： 式中叫做阻尼系数，它与物体的形状、大小及介质有关。对弹簧振子，在弹性力及阻力的作用下，物体的运动方程为 令，这里的为无阻尼时振子的固有角频率，称为阻尼因子，代入上式后可得在阻尼作用较小时，即，此微分方程的解为 其中此时阻尼振动的周期为当时，这种阻尼振动称为欠阻尼；当时，称为过阻尼；当时，称为临界阻尼。 4.2 程序编写 高数在物理上的一项重要应用就是微分方程。在以前的学习中，我们大多解的是一阶微分方程，而对二阶的总是避而远之。也许在学高数时，解二阶微分方程给你留下了痛苦的回忆，此时如果让你编写一个二阶微分方程的程序（如阻尼振动），你是否有种相逃的感觉？Wait amoment!上帝是青睐你的。Matlab只要用一个函数ode23就能解决你所有的问题。那我们现在就开始学吧！ode23函数可以用来解微分方程组，但只能是一阶的。所以对于二阶的微分方程，我们可以将其分解成两个一阶的方程来解。而且使用此函数，必须提前编写一个ode文件（ode文件）。 ode文件是一种函数M文件，即：在m-file编辑器中，文件必须以function开头，后面跟的函数名必须与将来保存的M文件名一致。函数文件第一句的具体形式为：functiony=average(x) 其中为y输出变量名，x为变量，average为函数名，保存时文件名必须是average.m。当函数有一个以上的输出函数时，输出参数包含在方括号里，变参包含在圆括号里，即函数文件的第一行的具体形式为： function[out1,out2,…]=average(in1,in2,…) 例：计算数组元素的平均值 解：functiony=average(x) y=sum(x)/length(x)注：上下两式等号前的符号"y"必须相同 ），再调用，其调用格式为：ode23('FUN',TSPAN,YO)其中，FUN为字符串，表示微分方程的ode文件名，TSPAN=[TO,TFINAL]表示积分区间，YO为初始条件。此函数表示在初始条件下Y0从T0到TFINAL对微分方程y'=F(t,y)进行积分。函数在图形窗口返回具体的数值解，则可以采用以下的格式： [T,Y]=ode23('FUN',TSPAN,Y0) 这样在T向量返回积分点，Y返回对应积分点的积分函数值。 例：编写m-file描述物体作阻尼振动的振动曲线，分欠阻尼、过阻尼和临界阻尼三种情况。 解：首先将阻尼振动的微分方程 分解为一阶微分方程：令 则原方程化为 打开m-file编辑器， 在其中编写ode文件为 functionf=zuni(t,x,flag,beita,wo); f=[x(2); -2*beita.*x(2)-wo.^2*x(1)]; 保存文件名为zuni.m 另外打开一个新m-file，在其中编写程序代码： beita=input('欠阻尼系数p='); wo=input('固有频率wo(pwo)='); [t,x]=ode23('zuni',[0:0.001:5*pi],[6,8],[],beita,wo); subplot(3,1,2),plot(t,x(:,1),'b'),title('过阻尼')%取x(1),绘制t-x(1)位移曲线，兰色 beita=input('临界阻尼系数p='); wo=input('固有频率wo(p=wo)='); [t,x]=ode23('zuni',[0:0.001:5*pi],[6,8],[],beita,wo); subplot(3,1,3),plot(t,x(:,1),'g'),title('临界阻尼')%取x(1),绘制t-x(1)位移曲线，绿色然后保存（如保存其名为zunizhend.m）， 再按F5运行。移至Matlab主窗口，输入数据， 按Enter键即得振动曲线 感觉如何？相信你通过以上的学习，已经揭开了编程的神秘面纱了.Goon! 5 受迫振动共振 5.1 讲授内容 一、受迫振动 振动系统在周期性的外力作用下，其所发生的振动称为受迫振动，这个周期性的外力称为策动力。假设策动力有简单的形式：,为策动力的角频率。由于同时受到弹性力和阻力的作用，物体受迫振动的运动方程为 令，，则上式可写成 式(1) 在阻尼较小的情况下，上述微分方程的解为 此解表示受迫振动是两个振动的合成。解的第一项为减幅振动，它随时间很快衰减，而第二项为一个稳定的等幅振动。一段时间后,第一项衰减到可忽略不计，所以受迫振动的振动表达式为 将此式代入式（1）可求得 式(2) 二、共振 由式（2）可知，稳态受迫振动的位移振幅随策动力的频率而改变。当策动力的频率为某一特定值时，振幅达到极大值。由式（1）利用求极值的方法可得振幅达到极大值时的角频率为 振幅的最大值为 可见，在阻尼很小（）的情况下，若策动力的频率近似等于振动系统的固有频率，位移振幅将达到最大值。这种现象就称为位移共振。用类似的方法可得，当策动力的频率正好等于振动系统的固有频率时，速度夫达到最大值，称之为速度共振。即：在小阻尼的情况下，当策动力的频率正好等于振动系统的固有频率时，振动系统的振幅最大，这种现象就称为共振现象。 5.2 程序编写 在上节我们学会了如何用Matlab编写阻尼振动的程序，这节我们主要编的是受迫振动的。由于受迫振动的方程只是在阻尼振动的方程上加了一个策动力F,所以二者的编程过程极其相似。相信通过这次的巩固，你以后就可以对这方面的编程更加得心应手了。例：弹簧振子在策动力、弹性力和阻尼力的作用下做受迫振动，其振动方程为 令，，则上式可写成 编写m-file描述物体作受迫振动的振动曲线，并讨论共振现象（）。 解：同样先将受迫振动的微分方程分解为一阶微分方程：令 则原方程化为 打开m-file编辑器， 在其中编写ode文件为 functiony=shoupozd(t,x,flag,f,wo,w,beita) y=[x(2); f*cos(w*t)-wo^2*x(1)-2*beita*x(2)]; 保存文件名为shoupozd.m。 另外打开一个新m-file，在其中编写程序代码： f=input('单位质量物体所受的最大策动力fo=');%输入函数所需的参数 w=input('策动力的频率w='); wo=input('振动物体的固有频率wo='); beita=input('阻尼因子beita='); [t,x]=ode23('shoupozd',[0:0.001:10],[6,8],[],f,w,wo,beita);%调用ode23,[0:0.001:10]为t的积分区间， %[6,8]为初始条件（这里只是任取的），即x=6,x'=8,[]说明以后是输入参数，f,w,wo,beita是参数名，其顺序必须与所调用的%ode文件的参数顺序相同 plot(t,x(:,1))%取x(1),绘制t-x(1)位移曲线 xlabel('t') ylabel('x') title('受迫振动') 然后保存其名为shoupozdong.m， 再按F5运行。移至Matlab主窗口，输入数据， 按Enter键即得振动曲线 从图中我们可以看到当时，开始的一段时间内，振幅缓慢减小，直到达到稳定状态，振动曲线变成一余弦曲线，振幅不再变化。而当，振动系统的能量得到了最大的补充，产生了共振。 我们只需在命令窗中输入shoupozdong就可以运行它，此次在输入数据 时，让即得共振的振动曲线. 我们可以很明显的看到，产生共振的系统的振幅渐渐增大。 Well,加上图的形象性，想必你对这节的学习在理论和编程上都没问题了吧。 6 非线性振动、混沌 6.1 讲授内容 现代科学技术研究发现，非线性是真实世界的普遍特性，非线性问题大量出现在自然科学、社会科学和工程科学中，并起着重要的作用。混沌的研究是20世纪物理学的重大事件，自从麻省理工学院的Lorenz教授在60年代进行了开创性的研究以来，已有更多的学者深入探索，逐步揭示了混沌运动的基本特征，即确定性中包含的非周期性和不可预测性，对初值的敏感性等等。 “混沌”译自英语中“chaos”一词，原意是混乱、无序，在现代非线性理论中，混沌则是泛指在确定体系中出现的貌似无则的、类随机的运动。简单地说，混沌理论所研究的是非线性力学混沌，目的是要揭示貌似随机的现象背后所隐藏的规律。 混沌的研究表明，一个完全确定的系统，即使非常简单，由于自身的非线性作用，同样具有内在的随机性。绝多数非线性动力学系统，既有周期性，又有混沌运动，而混沌即不是具有周期性和对称性的有序，又不是绝对的无序，而是可用种复杂的有序。总之，混沌是非周期的有序性。 例如，我们在学简谐振动时常用到的一种理想模型--单摆，在摆角很小的情况下，单摆做周期性的简谐振动，而在大摆角时就会出现混沌现象。 6.2 程序编写 单摆运动的微分方程为：,为摆长，为摆角。 (1)当较小时，运动方程可写为 (2)大角度单摆运动的微分方程为 编写m-file文件描述单摆在大角度和小角度的情况下的振动曲线()及相图() 解：step1.当较小时，先将此时的微分方程分解成一阶微分方程： ，令则原方程化为： 打开一个新的m-file，在其中编写ode文件为 functionT=danbaixfun(t,theta,flag,g,l) T=[theta(2); -g*theta(1)/l]; 保存文件名为danbaixfun.m。 step2.当为大角度时，同样先将微分方程分解成一阶微分方程： 令则原方程化为 另打开一个m-file，在其中编写ode文件为 functionT=danbaidfun(t,theta,flag,g,l) T=[theta(2); -g*sin(theta(1))/l]; 保存文件名为danbaidfun.m。 step3.再打开一个m-file,在其中编写程序代码,为了方便对比，我们将两种情况下的振动放在一起： %单摆在摆角较小的情况下作简谐振动其微分方程为ml*D(theta)^2=-mgtheta,l为摆长,theta为摆角 %编写文件描述单摆的振动曲线和相图 g=input('重力加速度g='); l=input('摆长l='); theta0=input('初位置的摆角theta(theta为小角度)='); [t,theta]=ode23('danbaixfun',[0:0.01:3*pi],[theta0,0],[],g,l);%调用ode23,[0:0.001:3*pi]为t的积分区间， %[theta0,0]为初始条件,[]以后是输入参数 subplot(2,2,1);plot(t,theta(:,1),'r'),xlabel('t'),ylabel('theta'),title('小角度振动曲线') subplot(2,2,2);plot(theta(:,1),theta(:,2),'b'),xlabel('theta'),ylabel('omega'),title('小角度相图') %单摆在摆角较大的情况下不再作简谐振动，其微分方程为ml*D(theta)^2=-mg*sin(theta),l为摆长,theta为摆角 %编写文件描述大角度单摆的振动曲线和相图 g=input('重力加速度g='); l=input('摆长l='); theta_0=input('初位置的摆角theta(theta为大角度)='); [t,theta]=ode23('danbaidfun',[0:0.01:3*pi],[theta_0,0],[],g,l);%调用ode23,[0:0.001:3*pi]为t的积分区间， %[theta_0,0]为初始条件,[]以后是输入参数 subplot(2,2,3);plot(t,theta(:,1),'r'),xlabel('t'),ylabel('theta'),title('大角度振动曲线') subplot(2,2,4);plot(theta(:,1),theta(:,2),'b'),xlabel('theta'),ylabel('omega'),title('大角度相图') 然后保存其名为danbai.m， 再按F5运行。 step4.移至Matlab主窗口，输入数据， 按Enter键即得振动曲线。 为验证混沌运动的无序性，我们重新输入一组数据， 得到另一 振动曲线， 比较第一次和第二次的曲线，可得单摆在摆角很小的情况下，其振动有规律可循，而在大角度时，振动呈无序状。为进一步说明问题，我们在输入一组数据（将小角度的振动初值改变，而大角度的振动初值只将初角位移加上0.001）， 得到第三组振动曲线。 比较第二次和第三次的曲线，你是否会奇妙的发现，对于大角度振动，仅是0.001的细微差别具有迥然不同的结果，这就是混沌的无序性和队初值的敏感性。 7 平面简谐波 7.1 讲授内容 一、平面简谐波 振动状态的传播过程称为波。若波源作简谐振动，这时的波动即为简谐波。当简谐波传播时，介质中的各个质点都作简谐振动，若其波阵面是平面，就称为平面简谐波。因为在理想的无吸收的均匀无限大的介质中，同一波阵面上的各点振动状态相同，所以在研究平面简谐波传播的规律时，只要讨论与波阵面垂直的任意一条波线上的传播规律就可以了。 设平面简谐波沿轴正方向传播，波速为u。取任意一条波线为轴，并取O为即坐标原点。已知原点O处的质点作简谐振动，适当的选取计时零点可令其初相为零，故其振动表达式为 表示处质点离开其平衡位置的位移。考察波线上另一任意点P处的质点运动情况，P点的坐标为。若波是在理想的无吸收的均匀无限大的介质中传播，则P点处的质点将以相同的振幅和频率重复O点处的振动，但时间要晚一点。振动从O点传播到点P需时，因而O点时刻的振动状态在时刻传到点P。即使说，点P处质点时刻离开平衡位置的位移等于O点处的质点在时刻离开平衡位置的位移，即： 上式给出了波线上任一点处的质点位移随时间的变化规律，通常称为波的波动表达式。至于一般情形，即已知轴上某一定点的振动表达式为 则相应的在轴上传播的简谐波的波函数为 上式中负号对应简谐波向轴正方向传播，而正号对应向轴负方向传播的简谐波。 二、平面简谐波和简谐振动 波是振动状态的传播，振动的波源是产生波的必要条件。平面简谐波是弹性介质中大量质点参与的一种集体运动形式。即某一时刻，由大量偏离平横位置的质点组成波形图。此外，波是向前传播的，而波线上的每一振动点并不随波迁移，它们只在各自的平衡位置往复振动。将x=代入波函数即得处点的振动表达式： 其中 7.2 程序编写 在前面的学习中，我们重点讲了如何用Matlab解微分方程,但总体用的都是顺序结构，即按先后顺序逐条执行命令。你可能认为Matlab易学正是因为这一点，如果你真这么想，那你就犯了一个的大错。Matlab易学并不等于它功能弱，因而它理所当然具有循环、选择及其他语句。现在我们就通过距离子，让你更多的了解它吧！ 例：一余弦波在弦上传播，其波函数为 式中和的单位为，t的单位为 1、试求其振幅、波长、频率、周期和波速。 2、分别画出对应和俩时刻弦上的波形图。 解：由已知的波函数求波动的特征量，我们一般采用比较系数法，由 上式说明此简谐波向正方向传播，将它与波函数的

标准

excel标准偏差excel标准偏差函数exl标准差函数国标检验抽样标准表免费下载红头文件格式标准下载

形式相比较得： ,, 且有： 接下来我们就运用循环和选择语句来编画波形图的这个程序。 打开m-file编辑器， 在其中编写程序代码： x=-0.25:0.001:0.25;%设定x的取值范围 fori=0:2%用循环语句令i分别取0,1,2 t=0.0025*i; y=0.002*cos(pi*5*x-200*pi*t);%此为时刻t的波函数 ifi==0;%用选择语句分别用不同的颜色线型画不同时刻的波形图 plot(x,y,'k-')%用黑色实线画t=0时刻的波形图 holdon%保存图形的命令，否则后一幅图会覆盖前一幅 gridon%绘制网格 elseifi==1;plot(x,y,'r--')%用红色虚线画t=0.0025时刻的波形图 elseifi==2;plot(x,y,'b-.')%用蓝色虚点线画t=0.005时刻的波形图 end%结束if语句 end%结束for语句 然后保存（如保存其名为bodong.m）， 再按F5运行,即得波动曲线。 此时，你是否对“万事都有联系”这句话有了更深的理解呢？因为计算机编程的各种语言都有相似的地方，相信你通过Matlab的学习必会在编程这方面长进不少！