Mathematica在实际问题中的若干应用王晓波.docx

学年论文（设计）题 目 Mathematica在实际问题中的若干应用 学生姓名 王晓波 学号 1109014090 所在学院 数学与计算机科学学院 专业班级 数学与应用数学1102班 指导教师 彭严 _ _ 完成地点 陕西理工学院 Mathematica在实际问题中的若干应用王晓波陕西理工学院 数学与计算机科学学院 数学与应用数学 数应1102班指导老师 彭严摘要：在简要介绍Mathematica的基础上，着重讨论mathematica在实际问题中的应用。利用mathematica的强大的数据处理能力在实际学习和生活中帮助人们进行更快的计算，进而帮助人们更快的解决问题。关键词：Mathematica，振动波问题，量子力学，测绘，物理极值。一 引言Mathematica是美国 Wolfrmn 公司研制开发的著名数学计算软件系统，很好地结合了数值和符号计算引擎、图形系统、编程语言、文本系统、和与其他应用程序的高级连接。自 1987 年发布系统的1.0 版本开始便迅速广为流传，后经不断改进和完善，1991 年与 1997 年又先后推出 2.0 版和 3.0 版，1999 年推出 4.0 版本，后来推出了 5.0( 2003) ，6.0( 2007) ，7.0( 2008) ，直到 2011 年推出了 8.0 中文汉化版本。该版本增加了 500 多个新函数，功能涵盖更多应用领域，并拥有更友好更高质量的中文用户界面、中文参考资料中心及数以万计的中文互动实例，使中国用户学习和使用 Mathematica 更加方便快捷。自从 20 世纪 60 年代以来，在数值、代数、图形、和其它方面应用广泛。Mathematica 是世界上通用计算系统中最强大的系统，现在，它已经被应用于科学的各个领域物理、生物、社会学、和其它。二 Mathemtica在实际问题中的应用1 振动波问题振动波问题中,涉及很多关于函数的极限、导数、微分、积分的运算,并且涉及到许多复杂的技巧变换和运算,学生往往要花费大量的时间来计算和练习,以致有的学生产生了畏惧和逃避心理. 学生通过老师的讲解,能够对概念的本质、来龙去脉和其中的思想方法达到一定程度的理解和领会后,而涉及表达式的运算完全可以借助Mathematica软件的命令来完成.1.1 求导数求显函数的低阶导数、一元显函数的高阶导数和二元显函数的偏导数分别有格式命令如下:D表达式,求导变量,m用于求m阶导数;Df,x,求f对x的偏导数dfdx;Df,x1,x2,.,求f对x1,x2,.的高阶偏导数dnfdx1dx2dxn;Df,x,m,求f对x的m阶偏导数mfxm.如下例DC1 Cosbt*x+C2 Sinbt*x+C3 Coshbt*x+C4 Sinhbt*x,x输出为-bt*C1*Cosbt*x+bt*C2 Sinbt*x-bt*C3 Coshbt*x+bt*C4 Sinhbt*x1.2 求不定积分和定积分求积分一直是学生反映比较花费时间的问题,Mathematica的积分命令会使学生感觉积分易如反掌.1)求不定积分的格式为Integrate被积表达式,积分变量.2)求定积分的格式为Integrate被积表达式,积分变量,下限,上限.3)求二重积分的格式为IntegrateIntegratef(x,y),y,h(x),g(x),x,a,b可类推到多重积分.也可用格式:Integrate积分表达式,积分变量1,下限,上限,积分变量2,下限,上限.如计算Dxy2d,其中D是由0x1和y= x2围成.可输入IntegrateIntegratex* y2, y, x2,2x),x,0,1)直接运行即得结果59120.1.3 求系统固有频率和主振型Mathematica的Eigenvalues和Eigenvectors两个命令可用来求系统固有频率和主振型.如图1所示系统,已知m1=m2=m3=m,k1=2k,k2=k3=k,试求出系统的固有频率和正则振型矩阵.图1 多自由度质量弹簧系统多自由度系统的广义特征值问题中,K=2M ,其中特征矩阵为M-1K,2是特征矩阵的特征值,是特征矩阵的特征向量,M为质量矩阵,K为刚度矩阵.图1系统中M矩阵和K矩阵如下M=m1000m2000m3=m000m000mK=k1+k2-k20-k2k2+k3-k30-k3k3=3k-k0-k2k-k0-kk运行如下命令 M=DiagonalMatrixm,m,m;K=3k,-k,0,-k,2k,-k,0,-k, k;A=InverseM.K;EigenvaluesA输出固有频率的结果为:2km,2k-3km,2k+3km再输入EigenvectorsA输出主振型的结果为: -1,-1,1,2-3,-1+3,1,2+3,-1-3,11.4惯性式测振仪的基本原理时,幅频特性曲线如下Ya=112-12+22 相对阻尼系数N取0,0.3, 0.4, 0.5, 0.6, 0.7,1,2, 5和10时的图形,可通过输入如下命令得到.0,0.3,0.4,0.5,0.6,0.7,1,2,5,10,lmda,0.3,10运行后可得到图2.图2 以对数坐标表示的幅频特性曲2 在测绘中的应用2.1 编程计算测绘数据计算中往往会遇到许多测绘计算符号,而利用Mathematica进行编程时可以基本保留原有的测绘符号,这给测绘编程计算带来非常大的便利,同时也增加了程序的易读性。以测绘中经常遇到的空间前方交会计算为例,如图1的前方交会图形,利用Mathematica的编程计算如下:图1 空间前方交会图形In=Xp=XACot+XBCot+YB-YACot+CotYp=XACot+XBCot+XB-XACot+Cot只需输入XA,YA,XB,YB,A,B的值即可进行计算。2.2 平差计算测绘数据处理中经常会遇到数据平差。而最小二乘法是使插值函数在节点上近似地接近测量数据点,使整体误差达到最小。常用的方法是求导数,解线性方程组,其表达式较复杂且计算量相当大,而在Mathematica环境中,可以通过一个函数来实现:FindFitdata, expr,pars, vars。这里data为待拟合的数据系列, expr为表达式,pars为待平差的参数,vars为平差中的自变量。以下举一个多项式形式二次曲线平差的实例。有一组二次曲线数据为: -0.2504,-0.7436,-0.0003,-0.4434,-0.5739,-0.0002 ,-0.6502,-0.4341,0.0006 , -0.8627,-0.3307,0.0003 ,-1.0709,-0.2683,0.0007 , -1.2670,-0.2502,0.0009 , -1.4402,-0.2759,0.0001 ,-1.5859,-0.3451,0.0007 , -1.6953,-0.4543,-0.0006, -1.7636,-0.5992,-0.0005 , -1.7891,-0.7727,0.0006 , -1.7704,-0.9683,-0.0004 ,-1.7093,-1.1769,-0.0008,多项式形式二次曲线表达式为: 1+a1x+a2y+a3x2+a4xy+a5y2。待平差参数为:a1,a2,a3,a4,a5,平差中自变量为x,y。Mathematica的平差计算程序如下:In1: =data = 略FindFitdata,1+a1x+a2y+a3x2+a4xy+a5y2,a1,a2,a3,a4,a5, x,y则得到平差结果:a1=1.17939. a2=1.40695, a3=0.430514, a4=0.405681, a5=0.431796可以看出,利用Mathematica进行平差计算,工作量非常小,且不易出差错。由于常规的测量平差模型v=AX-L也类似这种线性形式,因此同样可以利用FindFit进行平差计算。2.3 曲线拟合Mathematica提供的进行曲线拟合或曲面拟合或超曲面拟合的函数为:Fit data,funs,vars ,其中data表示观测数据,其形式有:data1=y1,y2,yn,此时观测数据为(xi,yi),;data2=x1,y1, x2,y2, xn,yn,即观测数据为(xi,yi);data3= x1,y1,z1, x2,y2,z2, xn,yn,zn ,即观测数据为(xi,yi,zi),此时可以确定z是变量x,y的二元函数。设有一组数据为DD:1,5,10,21,26,50,50,85,91,130,122,210,170,250,260,341,290,455,362,546,通过Fit函数拟合如下:FitDD, 1,x,x2,x,得到拟合方程为5.32018-1.88158x+1.3114x2。拟合结果用Mathematica的绘图函数Plot和ListPlot绘图表示如图2。图2 某曲线拟合图2.4 迭代运算测绘中特别是在大地测量当中经常会遇到迭代运算,虽然利用软件进行迭代编程较为方便,但也较易出错。Mathematica提供了非常丰富的线性和非线性方程解算函数,包括Solve,NSolve,NDSolve和FindRoo,t其中FindRoot函数可以解算非线性方程的数值解,还可以用来进行迭代计算。例如,大地测量中的空间直角坐标与大地坐标系转换似乎是一个业已解决而略显陈旧的问题,但如何方便快捷地进行互换必需进行适当编程和调试。笔者利用Mathematica的FindRoot函数则轻松实现了两组坐标的互换。已知椭球面上某点的大地坐标换算为空间直角坐标(正算)的公式为:X=(N+HAg)cosBcosLY=(N+HAg)cosBsinLZ=N1-e2+HAgsinB式中 N=a1-e2sin2B正算非常简单而又严密。而空间直角坐标换算为大地坐标(反算)就较复杂。利用Mathematica,只需加入编程语句:FinRooX-N+HAgcosBcosL=0,Y-N+HAgcosBsinL=0,Z-N1-e2+HAgsinB=0,B,B0,L,L0,H,H0输入初值B0,L0,H0即得到B,L,H的解。利用教材3中的例子进行计算。已知X=4898979.486, Y=2828427.125, Z=5656854.249。经Mathematica的FindRoot计算得到解为B=45b091 1.2907, L= 30b000 0.0000, H=1632465.5898,这和教材3中的迭代解完全一致。由此看出,大地测量中的类似迭代计算都可以利用Mathematica的非线性方程解算函数来实现。3 在量子力学课程中的应用在量子力学中，力学量算符可以用矩阵表示，如果两个矩阵 ABBA，我们称 A、B 矩阵相互不对易，AB = BA，称 A、B 矩阵相互对易。转置矩阵: 把矩阵 A 的行和列互相调换，所得出的新矩阵称为A 的转置矩阵A。共轭矩阵 A+: A+( A+)mn= (Amn)*= A*nm，m 列 n 行 n 列 m 行转成共轭复数。厄密矩阵: 如果 A+= A，则称 A 矩阵为厄密矩阵( 如果一个矩阵 A 和它的共轭矩阵相等) 。本征值方程在量子力学中有着很重要的作用，我们都知道给定算符如何求本征值与本征函数。大体上分为以下几步: 1) 先求用矩阵表示的本征方程; 2) 代入久期方程求得本征值的解; 3) 本征值代入本征方程求本征函数。在这里，我们举例说明如何利用 Mathematica 来求解以上提到的几个问题:例 1 已知体系的哈密顿算符H与某一力学量算符B在能量表象中的矩阵形式为:H=h-100010001,B=b200001010其中 和 b 为实常数，问1) H 和 B 是否是厄密矩阵;2) H 和 B 是否对易;3) 求算符B的本征值及相应的本征函数 .Mathematica 代码如下:H = h ( _ 1，0，0 ， 0，1，0 ， 0，0，1 _) ;B = b ( _ 2，0，0 ， 0，0，1 ， 0，1，0 _) ;HermitianMatrixQHHermitianMatrixQBH.B B.HEigenvaluesBEigenvectorsB运行结果如下:In15=H=h-100010001;B=b200001010;HermitianMatrix0HHermitianMatrix0BH.B B.HEigenvaluesBEigenvectorsBOut17= FalseOut18= FalseOut19= 0，0，0 ， 0，0，0 ， 0，0，0 Out20= b，b，2bOut21= 0， 1，1 ， 0，1，1 ， 1，0，0 从上面的结果我们可以看出，H 和 B 都不是厄米矩阵，H 和 B 是对易的，直接可以求出 B 的本征值和本征态。例 2 2 个自旋 1/2 的粒子相互作用，它们直接的相互作用的哈密顿量可以写成 H =Jx1x2，假设两个粒子的初态波函数处于 ( 0) = |，由于粒子直接的相互作存在，波函数会随时间演化，如何求解 t 时刻的波函数问题? 下面将用 Mathematica 来求解这个问题。假设 ( 0) = |，时间 t 演化为 ( t) = C1+ C2| + C3| + C4| 在量子力学中，力学量算符可以用矩阵表示，如果两个矩阵 ABBA，我们称 A、B 矩阵相互不对In10=1=0110;2=0110;a=c1c2c3c4;H=JKroneckerProduct1,2.a/ /MatrixFormOut12/ / MatrixFoorm=Jc4Jc3Jc2Jc1In19=DSolvec1t=Jc4t,c2t=Jc3t,c3t=Jc2t,c4t=Jc1t,c10=0,c20=0,c30=0,c40=0c1t,c2t,c3t,c4t,tOut19=c1t0,c4t0,c2tCosJt,c3tSinJt从上面的程序就可以看出态演化的结果。4 在物理极值问题上的应用在计算极值问题时，最简单的方法就是直接用极值函数进行计算或用图形直观表达。（1）Maximizef，cons，x，y，表示求自变量为 x，y，的函数 f 满足条件cons 的最大值；同理，Minimize f，cons，x，y， 表示求函数 f 满足条件cons 时的最小值。（2）NMaximize 与 Maximize 函数具有相同的格式和作用，表示求最大值，只不过输出结果形式不同，如 NMaximizeSinx，x的计算结果是1.， x1.5708。（3）FindMaximumf，x，x0表示从 x0出发求未知量 x 的函数 f 的一个极大值点和极大值；同理，FindMinimumf，x，x0表示求极大值点和极大值。（4）Maxx1，x2，和 Minx1，x2，分别表示求一组数的最大值和最小值。（5）Limitexpr，xx0表示求函数 expr当 xx0时的极限。（6）Plotf(x)，x，a，b表示绘制函数 f(x)在区间a，b上的图形。下 面 就 对 一 个 物 理 实 例 用Mathematica 软件进行分析：如图 1 所示，电源电动势 E=6V，内阻 r=3,研究电源的输出功率随外电阻 R 的变化而变化的情况，并分析在什么情况下电源的输出功率最大。一般解法：输出功率p=I2R=ER+r2 R=E2/(R+r2R+2r)。当 R=r2/R，即 R=r=3时,有最大值 PMax=3W。图1如果此题不能进行理论计算，可以用 Mathematica 进行计算分析，其程序代码和计算结果如下所示（其中括号中内容为笔者说明）：In1=PR=6R+32R(定义函数P=6R+32R)In2=MaximizePR,R,MinimizePR,R0,R(求P的最大值和最小值)Out2=3,R3,0,R0(结果表明，当R=3时,P最大值为3;当R=0时,P最小值为0)In3=LimitPR,R(求R趋向时P的值)Out3=0当R时,P=0(结果表明,当R=3时，P最大值为3;当R=0时,P最小值为0)图2In4=PlotPR,r,0,10,PlotRanegAll,AxesLabelR,P(画出P关于R的图像)Out4=Graphics(从图2中可以看出P随着R的增大而先增大后减小，也可以看出当R=3时，P有最大值3)可见，从理论分析、极值函数计算结果和图形中能得出相同的结论，即输出功率 P 随着 R 的增大而先增大后减小。三 结论Mathem