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 讲授内容
简谐振动是一种无阻尼振动，而实际上，任何振动物体都要受到阻力的作用，这种振动叫阻尼振动。阻尼振动有摩擦阻尼和辐射阻尼两种。实验指出，当物体以不太大的速率在粘性介质中运动时，介质对物体的阻力与物体的运动速率成正比，方向与运动方向相反，即：
式中叫做阻尼系数，它与物体的形状、大小及介质有关。对弹簧振子，在弹性力及阻力的作用下，物体的运动方程为
令后可得
，这里的
为无阻尼时振子的固有角频率，
称为阻尼因子，代入上式
在阻尼作用较小时，即，此微分方程的解为
其中动称为欠阻尼；当
此时阻尼振动的周期为时，称为过阻尼；当
当
时，称为临界阻尼。
时，这种阻尼振
高数在物理上的一项重要应用就是微分方程。在以前的学习中，我们大多解的是一
阶微分方程，而对二阶的总是避而远之。也许在学高数时，解二阶微分方程给你留下了痛苦的回忆，此时如果让你编写一个二阶微分方程的程序（如阻尼振动），你是否有种相逃的感觉？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(p
[t,x]=ode23('zuni',[0:0.001:5*pi],[6,8],[],beita,wo);%调用ode23,[6,8]为初始条件， %即x=6,x'=8,[]说明以后是输入参数，beita,wo 是参数名
subplot(3,1,1),plot(t,x(:,1),'r'),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,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、分别画出对应和俩时刻弦上的波形图。 解：由已知的波函数求波动的特征量，我们一般采用比较系数法，由
上式说明此简谐波向正方向传播，将它与波函数的标准形式相比较得：
,
,
且有：
接下来我们就运用循环和选择语句来编画波形图的这个程序。 打开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 的学习必会在编程这方面长进不少！