基于Mathematica和Matlab的典型数学物理方程解分布图像制作.doc

第 29 卷 第 7 期2012 年 7 月吉林化工学院学报Vol 29 No 7JOURNAL OF JILIN INSTITUTE OF CHEMICAL TECHNOLOGYJul2012文章编号:1007-2853( 2012) 07-0041-04基于和的典型数学物理MathematicaMatlab方程解分布图像制作陈殿伟1 ，杨海英2( 1 吉林化工学院 理学院，吉林 吉林 132022; 2 吉林市第五中学 物理组，吉林 吉林 132011)摘要: 通过计算机辅助软件 Mathematica 和 Matlab，把典型数学物理方程解的空间分布制作成三维图像关 键 词: Mathematica; Matlab; 数学物理方程中图分类号: O 411 1文献标志码: A典型数学物理方程包括: 波动方程; 热传导方程; 泊松方程，u = 0 为拉普拉斯方程，二者都是 稳定分布方程 这三个方程构成了“数理方程”的 主要内容 它们都是二阶线性方程，刻画了很多物 理现象的规律1-3波动方程分布图像1波动方程或称波方程是一种重要的偏微分方程，主要描述自然界中的各种的波动现象 其表达 式为:图 1 位移在时间和空间中的分布图像2( 1)utt = a u + f 热传导方程2如弦的 x = 0 端固定，x = l 端受迫作谐振动 F= Asint，弦的初始速度为零，求弦的振动 这个 定解问题是:当一个物体内部各点的温度不一样时，则热量就会从温度高的地方向温度低的地方流动，这2utt a uxx = 0 ，u | x = 0 = 0，u | x = l = Asint ，u | t = 0 = 0，ut | t = 0 = A x = l 端为非齐次边界条件种现象就是热传导问题 由于热传导过程总是表( 2)现为温度随时间和位置的变化，所以，解决热传导问题归结为求物体内部温度的分布问题 如果研究物质的扩散问题，物质的扩散是由浓度高的地这是一个不受外力作用的振动，但它 在 弦x = l的一端，有一个谐振源 F = Asint，而在初始方扩散到低的地方，且服从傅立叶热传导方程相类似的能斯特扩散定理，即物质的扩散与浓度的时刻，速 度 为 A 用行4，得到如下图形:Mathematica直 接 编 程 运变化成正比 与热传导方程类似，可以得到浓度 u满足热传导方程，所以热传导方程也称扩散方程其表达式为图 1 中，a = 2，A = 5，l = 1， = 1，x 0，1 ，t0，10 从图中看出，在弦 x = 0 的一端，由于没 有谐振源，因此，位移为 0; 而在 x = l 的一端，谐振 源 F = Asint，位移 u( x，t) 随着 F 的作用而变化= a2 u u( 3)t对于细杆导热问题 初始时刻杆的一端温度为零度，另端温度为 u ，杆上温度梯度均匀，零度0收稿日期:2012-04-14作者简介:陈殿伟( 1968-) ，男，吉林省吉林市人，吉林化工学院副教授，主要从事凝聚态物理方面的研究的一端保持温度不变，另一端跟外界绝热杆上温度 u ( x，t) 满足下列泛定方程和定解 条件化，因为杆上的每个领域温度已经达到平衡稳定分布方程分布图像32ut a uxx = 0 ，u | x = 0 = 0，ux | x = l = 0 ，在热传导问题中，如果温度分布稳定，热源强度 f( x，y，z) 不随时间变化，热传导持续下去，最终 将达到稳定状态，空间中各点的温度不再随时间( 4)2u | t = 0 = u0 x / l ( a= k / c)泛定方程和边界条件都是齐次的，用分离变数法 求 解 在 这 里 可 以 用 Mathematica 直 接 编 程4，就可以画出杆上温度 u ( x，t) 的分布图像， 如图 2 3 所示3变化 ，即 ut = 0，得到方程u + f( x，y，z)= 0 ，( 5)此方程称为泊松方程，如果没有热源，即 f = 0则得到拉普拉斯方程 在扩散问题中，浓度处于稳 定的状态考虑振动的平衡现象，同样得到稳定分 布方程3 1 拉普拉斯方程如带电的云跟大地之间的静电场近似是匀强 静电场，其电场强度 E0 是竖直的 水平架设的输 电线处在这个静电场之中，如图 4 所示 输电线是 导体圆柱，柱面由于静电感应出现感应电荷，圆柱 邻近的静电场也就不再是匀强的了 不过，离圆柱 “无限远”处的静电场仍保持为匀强的 现在研究 导体圆柱怎样改变了匀强静电场图 2 温度在时间和空间中的分布图 3 温度在时间和空间中的分布图 2 3 的视角不同，x 0，1 ，t 0，5 ，从两图形中可以看出杆上 x = 0 的一端，由于受其 边界条件的限制，u( x，t) 不随时间的变化而变化， 它仍然是零度; 在 x = l 的一端，t = 0 时 u( x，t) 在最高点，随着时间的增大，u( x，t) 迅速减小，变化 的幅度很大，然后趋于平缓，最后几乎没有什么变 化，但 t = 0 时，u ( x，t ) 是从杆上 x = l 的一端向 x = 0的一端传导的，只是变化幅度不是很大 从整 个图形来看，u( x，t) 除了在边界 x = 0 处，它都是 随着时间的增大而迅速减小，然后趋于平缓，到后 来根本没有什么变化 这很明显是一个热传导方 程的图形，温度从高点向低点飞快的传导，随着时 间的增长，杆上的温度越接近均匀，所以我们后来图 4 电场强度在空间的分布首先需要把这个物理问题表明为定解问题取圆柱的轴为 z 袖 如果圆柱“无限长”，那么，这 个静电场的电场强度、电势显然跟 z 无关，我们只 需在 oxy 平面上加以研究就够了如图 4 所示，正是 oxy 平面上的静电场，圆柱 面在 oxy 平面的剖口是圆 x2 + y2 = a2 ，其中 a 是 圆柱的半径柱外的空间中没有电荷，所以电势 u 满足二 维的拉普拉斯方程:uxx + uyy = 0( 柱外)( 6)导体中的电荷既然不再移动，这说明导体中第 7 期陈殿伟，等: 基于 Mathematica 和 Matlab 的典型数学物理方程解分布图像制作43全可以把导体的电势当作零，从而边界条件为面均匀地一层层切割 u 所得到的，所以，由等高线u | 2 + y2 = a2 = 0 ( 7)组成的图形就像是 u 在 oxy 平面上的倒影x按照分离变量法 u( x，y) = X( x) Y( y) 代入拉普拉斯方程固然不难把它分解为两个常微分方 程，但代入上述边界条件却只能得到:X( x) Y( 槡a2 x2 ) ，( 8)不能分解为 X( x) 或 Y( y) 的边界条件 事实上，既然边界是圆，直角坐标系显然是不适当的，必须采用平面极坐标系拉普拉斯方程在极坐标系中的表达式为2 u / 2 + ( 1 / ) u / + ( 1 / 2 ) 2 u / 2 = 0，图 5 极坐标中的电势分布图像( a)，( 9)由( 9) 式推得其电场强度的表达式为式中: 是极径， 是极角 “导体电势为零”就表明为齐次的边界条件:E = ( u cos usin)+ ( u sin + ( 10)u | = a = 0 u cos，在“无限远”处的静电场仍然保持为匀强的E0 由于选取了 x 铀平行于 E0 ，所以在无限远处，Ey = 0，Ex = E0 ，即 u / x = E0 ，亦即 = ( D2cosx + Ea 22 E2 ( sin x cos x) ) +000yysinx E0 cos 因而还有一个非齐次u = E0 x =的边界条件:( D) ( 13)0y由( 13) 可求出 E，利用 Matlab 编程，并绘制( 11)u | = E0 cos 求得柱外的静电势为:图形，如图 6 所示2u( ，) = D0 ln ( / a) E0 cos E0 ( a / )cos ，( 12)其中: E0 cos 项为原匀强静电场中的电2势分布 当 较大时， E0 ( a / ) cos 可以忽略，所以它是对圆柱附近匀强电场的修正项，原因是 受柱面感应电荷的影响 此外，还有 D0 ln ( / a ) 项，它的系数 D0 是任意常数，这说明包含着某种 不确定的因素，这个不确定因素在于问题提出时 根本没有说明导体柱原来所带的电 量，可 见 D0 ln( / a) 正是圆柱原来所带电量的影响( 由静 电学可知，D0 ln( / a) 正是均匀带电圆柱体周围的 静电场中的电势) 根据其最终解( 12) ，用 Matlab 直接编程并运 行5-6，如图 5 所示图 5 中的 D0 = 3，a = 1，E0 = 5，0，2， =1，10，这个图就像是一个圆柱上面顶着一 个圆面，柱和面接触的地方不动时，圆面整个弯曲 了一些，就形成上面的这个图 从图中看到其边界 是个圆，正好符合了题中给出的边界条件是个圆 的要求 从上往下看，电势分布的方向由高向低 的，但是，当电势降低到圆中心处时，圆的中心是空 的，因为图中圆的中心处正好是导体圆柱所在的地 方，导体内部是个等势体，故变化是零 图下面的线 条是电势的等高线，这些等高线是由平行于 oxy 平图 6 极坐标中的电场分布图 6 是关于电场强度方向的图形，其中 D0 = 3a = 1，E0 = 5，0，2，1，10 ，我们看到的这个图形是在极坐标下表示的图形，在这个图中，左边电力线从无穷远处向近处平行的行进着可是到达中心时，没有了电力线 这是因为导体圆柱体正在这个中心上，由于导体内部为等势体，所以电场强度为零，也就没有了电力线，图上的电力线从远到近方向行进，结束于导体的表面，后又从导体的表面出发，向远处行进泊松方程在圆域 0 上求解泊松方程的边值问题3 2u = a + b( x2 y2 ) ，u |根据其最终解:= c ( 14) = 0u = + = c + a ( 2 2 ) + b 2 ( 2 2 )也有两个相等的最小值，而电势是沿着电荷密度分布的，电荷密度越大的地方，其电势也越大，所 以，电势分布图中，有两个相等的最大值，和两个 相等的最小值 沿着 y 轴，有两个最大值，沿着 x 轴有两个最小值 由于边界条件的限制，其电势分 布边界为一个圆，而在这个圆的平面上，电势为 零 下面那些是等高线，它们是平行 oxy 平面切割 得到的00412( 15)cos2 利用 Matlab 直接编程5-6，就可以把函数 u绘成一个三维图形，如图 7 8 所示结论4主要针对数学物理方程中典型的数学物理方程，使用 Mathematica 和 Matlab 软件进行编程绘 制出位移的分布图像; 对齐次和非齐次的热传导 方程，利用了中的 PDE 工具箱，通过相应的数值 计算步骤，画出了温度的分布图像，对于稳定分布 方程，制作出了物理量( 标量) 和相应向量场的分 布图像图 7 电势分布图像参考文献:1梁昆淼 数学物理方法M 北京: 高等教育出版社，2010姚端正，梁家宝 数学物理方法M 北京: 人民教 育出版社，1997黄大奎，舒慕曾 数学物理方法M 北京: 高等教 育出版社，2001吴剑，胡波主 掌握 和 精 通 Mathematica4 0M北京: 北京邮电出版社，2002陆 君 安，等 偏 微 分 方 程 的 MATLAB 解 法M武汉: 武汉大学出版社，2001刘宏友，李莉，彭锋 MATLAB 6 基础及应用M重庆: 重庆大学出版社，20012图 8 电荷密度分布图像3图 7中，c= 3，a = 1，b= 4，0= 5， 80，2，1，10 ，这两个图像都是由极坐标转换成直角坐标，第一个是电势分布问题，而第 二个图像是为了更好地说明第一个图形而添加的 电荷密度分布问题，电荷密度 w = a + b ( x2 y2 ) 转换成极坐标为 w = a + b( 2 cos2) ，从电荷密度 的公式及其图像，可以看出 w 有两个相等最大值，456Fabrication of Distribution Image for Solution of Typical MathematicalPhysics Equation Based on Mathematica and MatlabCHEN Dian-wei1 ，YANG Hai-ying2( 1 College of Sciences，Jilin Institute of Chemical Technology，Jilin City 132022，China; 2 Physics Group，NO 1 Mid