mathematica软件在物理数值模拟问题中的应用举例.doc

本科生毕业论文设计题目 Mathematica软件在物理数值模拟 问题中的应用 作者姓名 指导教师 所在学院 专业（系） 班级（届） 完成日期 年 月 日 目录Mathematica软件在物理数值模拟问题中的应用1摘要1关键字1英文摘要、关键字1一.前言2二.Mathematica软件在求解物理问题上的应用31.利用Mathematica软件求力矩32.Mathematica软件在量子力学中力学量算符上的应用43.用Mathematica软件进行曲线拟合5三.Mathematica软件在物理问题模拟上的应用71.Mathematica建立单摆模型72.Mathematica软件在点电荷电场中的模拟应用123.Mathematica软件在电路分析与运算上的应用154.Mathematica软件在光的波动性实验中的模拟应用19四小结23参考文献24Mathematica软件在物理数值模拟问题中的应用物理科学与信息工程学院物理学专业 指导教师 郭志辉 作 者 杨 帆摘要：本文介绍Mathematica软件在计算力学力矩、量子算符上的几个具体实例，以及在物理实验的数据处理中的应用，尤其是在曲线拟合中的应用。同时我们也简单介绍了该软件在单摆运动、点电荷、电路分析及光谱分析等物理问题数值模拟中的一些应用，展示了Mathematica软件在数值计算和动态可视化方面的卓越能力。关键字：Mathematica软件；计算机模拟；数值处理Examples of the Application of Mathematica Software in Numerical Simulation of Physical ProblemsYang FanAbstract: This atticle shows the outstanding ability of the Mathematica software in numerical calculation and dynamatic visualization through some concret examples of Mathematica software. We show how to calculate the torque, the quantum operator and the data processing of physical experiment especially the curve fitting. In addition, some applications of the software in the numerical simulation of physical problems, such as the motion of a simple pendulum , the distribution of electric field of the system including two charged points , circuit analysis and spectral analysis in optical system. Key words:Mathematica;Computer Simulation;Numerical Processing一.前言随着科学技术的进步，特别是计算机信息处理技术的发展，数值模拟已成为继科学实验和理论分析后人们进行数据分析并认识世界的新手段。近年，人们越来越多的把现代计算机技术如Matlab，Java，Mathematica等软件应用到物理问题的处理中，极大程度上解决了物理问题在模拟和数值处理方面的难题。本文主要通过举例简单介绍Mathematica软件在物理问题数值模拟中的一些应用。Mathematica是由美国Wolfram Research公司于1988年发布的一款科学计算软件,与Maple、Matlab并称三大数学软件。该产品是以符号见长的数学软件，具有强大的数值计算和制图功能8，其操作简单、表现直观易懂。利用它提供的内部函数，可克服并解决物理数值计算中待测量多、数学公式复杂、不确定关系的推导和计算以及繁杂的图线绘制等问题。利用Mathematica软件来进行对物理问题的模拟和仿真，能够形象、直观、动态的展示物理规律和现象。二.Mathematica软件在求解物理问题上的应用1.利用Mathematica软件求力矩问题:已知一个具有单位质量的质点在一为F=(3t2-4t)i+(12t-6)j的力场中运动，其中t表示时间。设该质点在t=0时刻位于原点，并且速度为零。现求在t=2时刻该质点所受的对原点的力矩为多少。题目解析：受力F为矢量，在Mathematica中用列表形式表示力场矢量。由牛顿第二定律可知F=ma， ，其中r为质点的运动轨迹。因此求轨迹r可使用二次积分，先利用v=Fdt求运动速度，再利用r=vdt求运动轨迹。力矩是受力质点相对于原点的位置矢量r与力F矢量的矢积M，由公式M=rF求得1，对于求矢量的叉积，在Mathematica软件中使用的是Cross a , b 这个内置符号来计算。解答本题时，需建立空间直角坐标系。i、j、k分别表示x、y、z方向的单位向量。则受力F在z方向分量为0，表示为F=(3t2-4t)i+(12t-6)j+0k 。根据以上的分析，v、r在z方向的分量也应为0 。而力矩M为r、F的矢积，可以知道M在x、y方向分量均为0 。程序代码编辑： F=3 t2- 4 t , 12 t-6 , 0;/用列表表示矢量，定义受力函数/v=F 1 dt , F 2 dt , F 3 dt;/单位质量，F与a表达式相同， 第一次积分求运动速度/r= , , ;/二次积分求运动轨迹/t=2;/求t=2时刻的力矩/M=Crossr , F/用Cross函数求M=rF/Clear t , F , v , r , M /清除对F、v、r、M的定义/输出结果： 0 , 0 , - 40 /矢量形式，表示结果为- 40k/输出结果表示力矩M=- 40k ，力矩在x、y方向分量均为0 ，并与z轴负方向平行。与分析结果一致。而通过笔算检验后可知计算结果正确。2.Mathematica软件在量子力学中力学量算符上的应用量子力学是反映微观粒子运动规律的理论，随着量子力学的出现，人们对物质微观结构的认识也日益深入2。在量子力学中，由于多维空间和不同表象的存在，力学量算符常用矩阵的形式表现。使用Mathematica软件可以进行矩阵的运算，也能够检验矩阵是否为厄米矩阵。问题：一个已知体系中的哈密顿算符H与一个力学量算符B在能量表象中的矩阵表示形式为 H= , B= 其中和b均为实常数 ，问（1）H和B是否为厄米矩阵；（2）H和B是否対易3。题目解析：量子力学中表示力学量的算符都是厄米算符，表示厄米算符的矩阵就是厄米矩阵。在了解厄米矩阵的定义前，须先了解何为共轭转置矩阵。当矩阵A+满足A+nm=Amn，即矩阵A+的第n行第m列的元为矩阵A的第m行第n列的元的共轭复数时，称A+为A的共轭矩阵。而当A+=A，即一个矩阵的共轭矩阵为自身时称其为厄米矩阵2。当两个矩阵相乘后有AB=BA时，称矩阵A、B対易；当ABBA时，两矩阵不対易。Mathematica软件的内置符号 HermitianMatrixQ 就是专门用来检验矩阵是否为厄米矩阵的指令。虽然本题可根据定义很容易的判断出矩阵H、B并不是厄米矩阵，但是很多问题是不能直接判断的。而两矩阵的対易性可以通过定义来编写程序为HB-BH，根据输出结果是否为零矩阵判断。如果所遇到的矩阵较为复杂，还可用 RowReduce 函数给出行约化形式。程序代码编辑：H= (-1 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 1) ;/定义H矩阵形式/B=b(2 , 0 , 0 , 0 , 0 , 1 , 0 , 1 , 0) ;/定义B矩阵形式/HermitianMatrixQH /检验H是否为厄米矩阵/HermitianMatrixQB /检验B是否为厄米矩阵/RowReduceH B - B H/通过HB-BH是否为零矩阵判断対易性/ClearH , B/清除H、B定义/输出结果为：FalseFalse0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0第一个Flase表示H不是厄米矩阵，第二个Flase表示B不是厄米矩阵。第三行表示HB-BH的结果矩阵，根据输出矩阵可知HB-BH=0，即两矩阵対易。3.用Mathematica软件进行曲线拟合在物理问题中许多物理量都是连续变化的，然而事实上很难捕捉到连续的数值。因此在物理实验中，常常采用的是测量一系列样本数据进行曲线拟合再进一步进行数据分析和现象解释的实验数据处理方法，尤其是二次拟合图像是经常需要的。通过笔算来得到拟合公式并制图较为复杂，利用Mathematica软件可以非常方便的对离散数据进行曲线拟合。利用Mathematica软件绘制图线的方法如下：1.列出需要绘制的点。data =x1 , y1 , x2 , y2 , x3 , y3 , . ;2.绘制点。g= ListPlotdata 还可以用内置符号PlotStyle来规定数据点的颜色、大小等样式。3.曲线绘制。（1）用线段连接需要绘制的点。g0=ListLinePlotdata（2）绘制数据点的直线拟合图像。l1=Fitdata,1,x,x 用来计算数据点一次拟合的直线方程。g1=Plotl1 , x , xmin , xmax 以绘制数据点的直线拟合图像。（3）绘制数据点的二次拟合图像。l2=Fitdata , 1 , x , x2 , x 用来给出数据点进行二次拟合的曲线方程。g2=Plotl2 , x , xmin , xmax 用来绘制数据点的二次拟合曲线。如需要用一张图像来展示多个图像，如需要数据点、二次拟合曲线用一个图像同时展现，可以利用内置符号Show 来实现。图1 曲线拟合程序编辑与输出图像中横、纵坐标分别为x、y ，粗点表示原始数据点，实折线是由初始数据点直接连接所成，点虚线 表示一次拟合图像，虚线 表示二次拟合图像。Out4、Out6分别输出一次、二次拟合方程。上图所展示的是某一组数据的数据点、数据点连线、直线拟合图像、二次拟合曲线图像在同一坐标系内输出的绘制效果。仅用来展示Mathematica软件的曲线拟合功能和绘图功能，其数据与图像并不代表任何一实际物理问题和物理实验数据。一般来说，在具体的物理实验中是不会出现要求绘制如此多的图像的情况的。三.Mathematica软件在物理问题模拟上的应用1.Mathematica建立单摆模型4在所有不同类型的机械系统中，单摆运动是典型而具有代表性的运动形式之一，长久以来人们对单摆运动进行了大量的研究。理想的单摆运动是无阻力的简谐运动，受阻力影响的阻尼摆的研究也有了非常迅速的发展。对于无阻力的理想单摆，最大摆动角度小于90并且摆线不可伸缩，故此摆球是没有沿摆线方向的“径向运动”的，有且仅有与摆线方向垂直的“切向运动”。有运动方程 即 - mg sin(t)=mL(t) ，其中，G= -mg sin(t)是重力的切向方向分力， 是摆球的切向速度，并规定中间的竖直方向表示摆角为零，右边摆角为正，左边为负，切向的正方向是摆角的增大方向。消去摆球质量m后，有 ，这便是最简单的无阻单摆运动方程。可以看出，只要解出摆角(t)，不仅可以知道单摆在各个时刻的位置，还可以得到它运动的速度和加速度等物理量。因此，了解单摆运动的关键就是解出该二阶微分的运动方程。但是，这是一个非线性二阶微分方程，笔算时只能求得在小摆角情况下的近似解，就是指最大不超过5的摆角，此时有sin，这样就实现了方程的线性化，为(t)+2(t)=0，其中 表示单摆运动的圆频率。若选定了运动的初始条件即摆球在t=0时刻的位置和速度，相当于定出振幅m和初始相位角，就可以完全确定摆角(t)了。下面给出一个可计算无阻单摆运动在任意小摆角情况下的数值解的程序。这个程序利用函数NDSolve 解出单摆的微分方程，再用作图指令Plot 画出摆角(t)随着时间变化而变化的图形。为了证明在小摆角情况下有sin，在程序中将近似前后所得到的(t)绘制在同一张图上进行对比。初始条件规定为：初始摆角=0，摆球初始速度v0 = 0.2m/s，摆线长L=1.0m，因此可将初始角速度换算为 。而因为v0 不符合Mathematica软件对符号定义的形式，故使用v0代替。本文中其它类似符号也做相同处理。程序编辑为：g=9.8; L=1.0; =g/L; v0=0.2; 0=v0/L;/定义常数量和初始条件/s=NDSolvet= = -2 Sint , 0= =0 , 0= = 0 , , t , 0 , 10/解微分方程/=/.s1 ;/将微分方程的解替换给/s1=NDSolve1t= = -2 Sin1t , 10= =0 , 10= = 0 , , t , 0 , 10 /解近似后的微分方程/1=1/.s11 ;/将近似微分方程的解替换给1/l1=Plott , t , 0 , 10/绘制(t)图像/l2=Plot1t , t , 0 , 10 , PlotStyle Dashing0.04 , Thickness0.006/绘制1(t)图像/Showl1,l2,AxesLabelt,/将两图像绘在同一坐标系内/Clearg , L , , v0 , 0 , s , s1 , , 1 , l1 , l2/清除函数定义和赋值/输出结果为：图2 无阻摆摆角输出图像Out2表示函数(t)的数值解是用插值函数给出的。同理，Out4是使用sin进行近似后所求得的1(t)。插值函数是通过一个内部规定的多项式函数将指令NDSolve 所解出的各个分立点 (ti , i) 连接起来。输出图像Out8中横纵坐标分别表示时间t、摆角角度。实线图像表示使用标准方程所解出的摆角随时间变化图像(t) ，虚线表示进行近似后所得的摆角图像1(t) 。因为在实验中通常需要测量30或50个周期，修改程序后可输出30个周期的图像，从图3可以看出(t)与1(t)在30个周期的时间差几乎为零，所以在笔算时使用sin进行计算。 图3 30个周期的无阻摆摆角输出图像而在下面中只将初速度改变为2m/s 、4m/s，来观察一下在较大初速度时图像的变化。图4 v0=2m/s时无阻摆摆角输出图像图5 v0=4m/s时无阻摆摆角输出图像图4、图5分别表示初速度v0为2m/s（摆角为37.3）、4m/s（摆角为79.4），其它条件不变时无阻摆摆角的输出图像。观察可知30个周期的时间绝对误差分别为t=2s、8s，相对误差为EN=3.4%、13.6% 。对比图2可以知单摆的近似求解只能在小摆角下完成。用Mathematica软件还可以解析地计算出如单摆的周期、振幅等物理量。在无阻摆情况下，根据能量守恒，有 并以此计算出最大角度m 。程序编辑为：v0=0.2; g=9.8; L=1.0;acos=1- v02/(2g L);Printmax= , ArcCosacos , rad/求余弦的反函数并输出/Printmax= , ArcCosacos*180/ , /将余弦的反函数输出为“度”/Clearv0 , g , L , acos输出结果为： max = 0.0638985 rad max = 3.66112 即为最大摆角角度。还可以通过求插值函数(t)的极大值的方法求角振幅。程序编辑为：g=9.8; L=1.0;=g/L; v0=0.2;0=v0/L;s=NDSolvet= = -2*Sint , 0= =0 , 0= =0 , , t , 0 , 10 ;=/.s1;FindMaximumt , t , 0.5/求局部极大值/Clearg , L , , v0 , 0 , s , 输出结果为：0.0638985 , t0.5019 在这个程序中指令FindMaximum 是用来求函数局部的最大值的，它是求解出角振幅的关键。从图2可以看出摆角的第一个极大值大约出现在0.5 s左右，因此取t值为0.5 。以上两种方法的输出结果分别为0.0638985 rad 和0.0638985 rad ，计算结果高度一致，可见计算的准确性。如果可以求出相邻的两个摆角极大值出现的时间，那它们之间的时间差就是单摆运动的周期。编辑程序为：g =9.8 ; L=1.0; =g/L ; v0=0.2 ; 0=v0/L;s=NDSolvet= = -2*Sint , 0=0 , 0=0 , , t,0,10 ;=/.s1 ;m1 = FindMaximumt , t , 0.5/求第一个摆角极大值/m2 = FindMaximumt , t , 2.5/求相邻的第二个极大值/T=(t/.m22) - (t/.m12) ;/求两个相邻极大值时间差/PrintT= ,T , s/输出时间差，即周期/Clearg , L , , v0 , 0 , s , , m1 , m2 , T输出结果为：0.0638985,t0.50190.0638985,t2.5095T=2.0076s从输出结果可以看出FindMaximum 函数两次求得的极大值是一样的，表明振动是等幅度的，也说明所求T确为单摆的运动周期。对于阻尼摆，计算过程与无阻摆相似，不过单摆在切向方向的运动方程变成了 其中为阻力系数，表示与阻力有关的参数。下面给出一个阻尼摆的摆角程序编辑与输出图像，同样将近似解与标准解通过同一图像表现。而其振幅、周期的计算与计算无阻摆的相应物理量的方法相同，在此不再赘述。图6 阻尼摆摆角输出图像 图像的横纵坐标分别表示时间t、摆角 ，从图像可以看出，阻尼摆的最大摆角越来越小，表示其在最低点的速度也在逐渐减小，最终静止不再摆动。2.Mathematica软件在点电荷电场中的模拟应用电力线是一种表示电场分布的数学上的线，是真实存在却无法触摸的一种特殊的场物质。规定电力线上每点的切线方向为该点的电场强度矢量的方向，因此又被称为电场线。现在普遍都是用电力线表示电场的分布情况。应用Mathematica软件的内置函数和指令可形象直观的模拟出两个点电荷的电场线和电势的分布图形。由于两个点电荷所激发出的电场和电势的分布是关于两电荷的连线旋转对称的，因此只需画出两点电荷在电荷连线的一个平面（设其为xOy平面）的电场和电势分布图形即可表示整个空间的电场和电势分布情况。现设两个点电荷的电荷量分别为q1、q2，相距间隔为2r0，以两电荷连线所在直线为x轴，将两个点电荷连线的中点设为坐标原点并建立xOy平面直角坐标系。由点电荷电势的分布情况易得两个点电荷在这个平面坐标系中的电势场的分布函数为 5 ，其中0为真空介电常量，其近似值为0= 8.8510 - 12C2/(Nm) 。电场的分布函数为负的电势梯度函数即Eu。程序设计思路为：按电势的分布式定义为电荷量及坐标的函数，并且将两点电荷的电荷量q1、q2为动态变量，使程序设置为一个模型，利用Mathematica软件的内置符号 Manipulate 来创建一个交互式控件，实现模型的建立。软件内置指令 ContourPlot 用于生成关于x和y的函数f的等高线图，可用此符号来绘制等势线分布图。另一个内置指令StreamPlot 是用于生成以 x 和 y 的函数表示的矢量场的流线图，可用这个符号绘制负的电势梯度函数（即电场分布函数）所对应的电场分布图。并可通过内置符号Showa,b来实现将等势线和电场线分布图共同绘制在一张图上。因为在执行指令Manipulate 创建控件时无法进行函数和符号的调用，所以只能完整的输入函数。通过计算得到电场在x、y方向的分量为： , 。执行程序时取r0=1.010-10m 。程序编辑为： 图7 两个点电荷电势电场模拟程序输出的图像通过取值可以展示两个点电荷的电势电场情况:图9 两个等量同号点电荷电势、电场分布 图8 两个等量异号点电荷电势、电场分布 图8展示的是电荷量分别为q1= 3，q2= - 3时的等势线和电场线，即两个等量异号的点电荷的等势线和电场线分布图。图9展示的是电荷量为q1= q2= - 1时的等势线和电场线，即两个等量同号的点电荷的等势线和电场线分布图。两图中带有箭头的曲线为电场线，圆环表示等势线。以上两图均以点电荷所在位置为圆心，右边电荷带电荷量为q1，左边电荷带电量为q2。可动态改变q1、q2即两点电荷的电荷量的大小。从而可以得到两点电荷在不同带电量的情况下的电势和电场的分布图像，从而达到模型建立的目的。3.Mathematica软件在电路分析与运算上的应用4刚才通过一个小例子简单展示了Mathematica软件在静电学的应用，Mathematica软件也可应用在“动电”领域，即电路的分析与运算。图10 直流电桥在大学物理电学实验中，电桥实验电路是比较常见的基础电路。图10是一个最简单的直流电桥的电路图。R1R4四个电阻组成电桥的四个臂，点b与点c之间接入的灵敏电流计G用来检验电桥是否进入“平衡状态”，指标为电流计的指数为零。在给出电源电压、四个臂上的电阻R1R4及灵敏电流计的内部电阻Rg的情况下，便可计算出流过各臂的电流i1i4、电流计G的电流ig等其它未知物理量。为了方便解题，使用“电位”来代替“电压”作为基本物理量。定义d点电位为零，Va、Vb、Vc 分别表示点a、b、c的电位，并令Va=V。则R1两端电压为Va-Vb=V - Vb，根据欧姆定律及电荷守恒定律，还可依次得到以下方程组： Vb - 0 = i2*R2 Vc - 0 = i4*R4 V - Vb = i1*R1 V - Vc = i3*R3 Vb - Vc = ig*Rg i1 = i2 + ig i4 = i3 + ig该方程组内含2个未知电位、5个未知电流，共计7个未知量，而组内有7个方程，因此该方程组可以求解。可用Mathematica软件的函数Solve 求解代数多项式来解方程组，并使用FullSimplify对计算结果进行化简。编辑程序为：equ=vb= =i2*R2 , vc= =i4*R4 , vb - vc= =ig*Rg , v - vb= =i1*R1 , v - vc= =i3*R3, i1= =i2+ig ,I 4= =i3+ig;/定义方程组/vars=i1 , i2 , i3 , i4 , ig , vb , vc ;/设定待求量/Solveequ,vars / FullSimplify/求解方程组/Clearequ , vars输出结果为：图11 直流电路输出输出结果是7个待解量的解析表达式，虽然看似较为繁琐，但它们都经函数FullSimplify 化简过，若笔算求解将会发现十分麻烦。下面就图12所展示的交流电路，继续Mathematica软件在电路分析上的应用。在交流电路运行后，不改变输入电压的频率和大小的话，电路将会达到一个稳定的状态，各个物理量达到平衡，即为稳态。此时各个量的解称为稳态解，求该电路的稳态解，实际上就是进行交流电路的稳态分析。图12 交流电路设交流电源电压是余弦变化的，即V(t)=V0*cost 。式中，V0表示交流电的振幅， 表示交流电的圆频率。将各个电路元件上的阻抗和电源电压写成复数形式。5个元器件的阻抗分别为：z1=R1 , z2=R2 , z3=R3 , z4=jL , z5=1/(jC) ，电源电压为V(t)=V0ejt ，其中，j为复数因子。同样用电位替代电压对电路进行分析，取d点电位为0 。可列出方程组： Vb - 0 = i2*R2 Vc - 0 = i4*(jL) V - Vb = i1*R1 V - Vc = i3*R3 Vb - Vc = ic/(jC) i1 = i2 + ic i4 = i3 + ic根据以上方程组，即可编写程序求解7个未知变量。还可定义电容上的电压uc，并求出它随时间变化的关系式。程序编辑与输出结果为：图13 交流电路稳态解 输出结果中Out7是7个未知变量的解析表达式，Out8是电容两端的电压uc的复数表示式，Out10则是uc随时间变化的关系式uc(t)。依旧采用图12所示的电路，利用Mathematica软件对其进行交流电路的瞬态分析。为简化计算，令电路处在“零状态”情况下，即所有物理量的初始值均为0 。那此时电源电压应为正弦变化，设为V(t)=V0*sint。求解电容的电压随时间变化的规律，电容和电感上的电压分别满足以下微分方程： , 。在下面的求解程序中，给出各元件的具体参数，均采用国际单位，电源频率为220Hz，求解0.1s内的电容电压变化规律。因为规定电路从“零状态”出发，所以电源电压和所有的7个函数的初态都为0，并通过制表函数Table 定义。编辑程序为：R1= 10 ; R2= 20 ; R3= 30 ; c = 10-4 ; L= 0.1 ;/给定元件参数/ = 2* *220 ; u0=1 ; tm=0.1 ;/方程条件设定/i10 , i20 , i30 , i40 , ub0 , uc0 , ic0=Table0 , 7;/定义初态为“零状态”/ut_: = u0*Sin*t ;/定义输出公式/equ=c*Dubt-uct, t = = ict, uct= =L*Di4t, t , i1t= =ict+i2t ,i4t= =i3t+ict , ut-ubt= =i2t*R2 , ubt= =i1t*R1 , ut-uct= =i3t*R3 , i10= =i10 , i20= =i20 , i30= =i30 , i40= =i40, ic0= =ic0 , uc0= =uc0 , ub0= =ub0;/定义求解方程与条件/s=NDSolveequ , i1 , i2 , i3 , i4 , ic , uc , ub,t , 0 , tm , MaxSteps;/求解方程/Plot(ubt-uct)/.s1 ,t , 0 , tm , PlotRange0.5,AxesLabelt,uc/绘图/Clear i1 , i2 , i3 , i4 , ic , uc , ub , u , R1 , R2 , R3 , c , L , , u0 , tm 图14 交流电路瞬态分析输出图像： 从输出的图像来看，电容电压在开始时内并不稳定，振幅有大有小，这是一段过渡状态，大约有0.025s。后来才达到稳定振荡的稳态。这样的输出结果符合实际情况。4.Mathematica软件在光的波动性实验中的模拟应用 光，在我们的生活中无所不在。在科学上的定义，光指的是特定波段的电磁波谱。我们都知道光具有粒子性与波动性，这被称为波粒二象性。 光的波动性从十七世纪被荷兰物理学家惠更斯提出后，直到十九世纪由于光的干涉和衍射现象的发现才得到人们的广泛认可。而光的干涉和衍射实验不仅是证明光具有波动性的强有利证据，也因此成为了经典的光学实验。在以往的物理教学中，光的干涉和衍射都是通过实验现象来观察学习的。现在，利用Mathematica软件可以根据理论公式绘制出物理现象的图像，形象、直观地来解释实验。对于光的干涉现象，需要了解有同一相位差的多光束叠加的理论知识。设有频率相同的光N束，它们都沿同一直线振动，而且各相邻的光束之间相位差相同。则各束光用复数表达分别为：A1eit , A2ei(t+) , A3ei(t+2) , A4ei(t+3) , . , ANeit+(N-1) 。实验中的多光束一般都是将一束光通过分光得到的，为简化同时也使计算更符合实验实际，我们完全可以假设每束光都有相同的振幅为A1 = A2 = A3 = A4 = = AN = A0 。这样多束光叠加的结果可以用复数Aei(t+0)来表示，振幅表示为A，初相位表示为，则有 。因为光强正比于振幅（实数）的平方，而上述式子是用复数表示的，所以光强正比于它与它的复共轭的乘积，即 。再根据欧拉公式eiN+e-iN=2cosN，可将上式改写为 6 。Mathematica软件的内置符号DensityPlot 是用来绘制关于x和y的某函数f的密度图形的指令7。根据光强公式，可以用它来绘制光的干涉、衍射现象的图像。为使图像更加美观，利用指令PlotPoint500取500个绘图点，指令AspectRatio1/5设定图像纵横比为1/5。对于光的双缝干涉，实际上是取N=2时的多光束叠加。为方便编程，令振幅A=1，则可编辑程序为：a=Sin2/Sin/22 ;/定义光强公式/Plota , , -8 , 8 , AspectRatio1/5 /输出光强分布图/DensityPlota, , -8 , 8 , y , -1 , 1 , PlotPoint500, AspectRatio1/5/输出干涉图像/Cleara输出图像：图15 双缝干涉光强和图像上图的Out2表示的是杨氏双缝干涉的光强分布图，Out3是双缝干涉现象的图像。可以看出图像与实验结果的一致性，说明了程序编辑即输入公式的正确性，用熟悉的现象反过来证明理论解释的正确。对于光的单缝衍射，衍射角为的所有次波在一个观察点叠加起来的合振幅为 6 ， 其中b表示狭缝宽度，为波长。为方便计算，令 。则光强为I = = sinc2u 。同样取光束的振幅A0=1。编辑程序为：b=Sincu2;Plotb ,u , -8 , 8 , Plot