数学物理问题的 Mathematica 求解.pdf

72 第七章 数学物理问题的 Mathematica 求解 7 1 偏微分方程的解析解 1 一阶偏微分方程的通解 Mathematica 中的 DSolve 命令也可以用来求偏微分方程的通解 使用格式为 DSolve 偏微 分方程 未知函数 自变量 例如 求一阶偏微分方程0 xy uu 的通解 命令语句为 DSolve D u x y x D u x y y 0 u x y x y 输出结果为 u x y C 1 x y 注意其中的 C 1 是一个任意函数 而不是常微分方程通解中的任意常数 上面的结果可以写成 我们更习惯的方式 1 u x yC yx 上例中的一阶偏微分方程是线性齐次的 现在我们再看一个线性非齐次的偏微分方程 1 xy uuxy 求解的命令语句为 DSolve D u x y x D u x y y 1 x y u x y x y 输出结果为 u x yD Log xD Log yD xC 1D x yD yC 1D x yD x y 上面的结果可以写成 11 1 lnln lnln xyx C yxy C yxxy u x yC yx xyxy 结果中的第一项为对应齐次方程的通解 第二项为非齐次方程的特解 下面再看一个非线性的一阶偏微分方程 2 xy uuu 求解的命令语句为 DSolve D u x y x D u x y y u x y 2 u x y x y 输出结果为 u x yD 1 x C 1D x yD 上面的结果可以写成 1 1 u x y xC yx 73 2 二阶偏微分方程的通解 现在 我们来看二阶偏微分方程的求解 例如求波动微分方程 2 ttxx ua u 的通解 命令语 句为 DSolve D u t x t t a 2 D u t x x x u t x t x 输出结果为 u t xD C 1DB a2t xF C 2DB a2t xF 注意现在通解中有 C 1 和 C 2 二个任意函数 由于 Mathematica 不知道参数 a 的性质 因而未 对 2 a进行化简 设 a 0 写成习惯的形式为 12 u x tC xatCxat 按照方程的要求 这两个任意函数都应该存在二阶偏导数 我们也可以用 DSolve 来求拉普拉斯方程0 xxyy uu 的通解 命令语句为 DSolve D u x y x x D u x y y y 0 u x y x y 输出结果为 u x y C 1 x y C 2 x y 注意现在通解中的二个任意函数 C 1 和 C 2 都应该是解析函数 不是所有的偏微分方程都可以用 DSolve 来求出通解的 例如 我们所熟悉的热传导方程 2 txx ua u 虽然存在形式上的通解 2 4 1 2 x s at u x ts eds at 但是用 DSolve 命令 DSolve D u t x t a 2 D u t x x x u t x t x 得到的结果为 DSolve uH1 0L t xD a2uH0 2L t xD u t xD 8t x D Mathematica 把输入的内容原样输出 表明它无法按要求执行 DSolve 通常也不能求出非齐次二阶线性偏微分方程的通解 例如 非齐次波动方程 2 ttxx ua ux 的求解命令为 DSolve D u x t x x D u x t t t x u x t x t 输出结果为 DSolve uH0 2L x tD uH2 0L x tD x u x tD 8x t D 3 特解的计算 74 一般来说 DSolve 命令也不能求出偏微分方程的特解 例如 对于波动方程的柯西问题 22 00 0 2 ttxx xx ttt uux ueuxe u t xD 8t x F Mathematica 也是把输入的内容原样输出 又如 对于拉普拉斯方程问题 2 0 0 2 xxyy y uu ux 用 DSolve 命令 DSolve D u x y x D u x y y 0 u x 0 2x 2 u x y x y 得到的结果为 DSolve 8uH0 1L x yD uH1 0L x yD 0 u x 0D 2x2 u x yD 8x y D 7 2 数学物理问题的数值求解 1 典型案例 虽然 Mathematica 的 DSolve 命令在求偏微分方程的通解或者特解中表现不佳 但是数值求 解微分方程的命令 NDSolve 却可以较好地处理偏微分方程的定解问题 例如 在时间范围 0 0 3 t 内数值求解热传导方程的定解问题 01 0 01 0 0 1 txx xx t uux uu uxx 命令语句为 NDSolve D u x t t D u x t x x u x 0 x 1 x u 0 t 0 u 1 t 0 u x 0 1 t 0 0 3 输出的结果为 75 u InterpolatingFunction 0 1 0 0 3 这表明 Mathematica 得到了一个在 0 1 0 0 3 xt 范围内的插值函数 u x t 在此基础上 我们可以用三维作图命令 Plot3D 来画出插值函数的图像 命令语句为 Plot3D Evaluate u x t First x 0 1 t 0 0 3 输出结果为 0 0 25 0 5 0 75 1 0 0 1 0 2 0 0 1 0 2 0 0 25 0 5 0 75 我们也可以用一系列的时间抽样图像来显示温度 u 随着时间的演化情况 如果抽样的间隔 为 0 1 命令语句如下 Table Plot Evaluate u x t First x 0 1 t 0 0 3 0 1 输出结果为 0 20 40 60 81 0 05 0 1 0 15 0 2 0 25 t 0 20 40 60 81 0 02 0 04 0 06 0 08 t 76 0 20 40 60 81 0 005 0 01 0 015 0 02 0 025 0 03 0 035 t 0 20 40 60 81 0 002 0 004 0 006 0 008 0 01 0 012 t 虽然所显示的温度分布的图形相似 但是数值的标度范围在迅速变小 这说明物体中的温度随 时间迅速下降 我们可以调整 Plot 中的参数来保持标度范围不变 如果要取作图的纵坐标范围固定为 0 0 25 可以设置参数为 PlotRange 0 0 25 命令语句为 Table Plot Evaluate u x t First x 0 1 PlotRange 0 0 25 t 0 0 3 0 1 输出的结果为 0 20 40 60 81 0 05 0 1 0 15 0 2 0 25 t 0 20 40 60 81 0 05 0 1 0 15 0 2 0 25 t 77 0 20 40 60 81 0 05 0 1 0 15 0 2 0 25 t 0 20 40 60 81 0 05 0 1 0 15 0 2 0 25 t 这样就可以明显地看出不同时刻的温度分布情况的变化 即温度分布随着时间的演化情况 我们也可以固定某个位置 观察温度随着时间变化的规律 如果固定位置的坐标为 0 1 0 5 0 9 对应的命令语句为 Table Plot Evaluate u x t First t 0 0 3 PlotRange 0 0 25 x 0 1 0 9 0 4 输出的结果为 0 050 10 150 20 250 3 0 05 0 1 0 15 0 2 0 25 x 0 050 10 150 20 250 3 0 05 0 1 0 15 0 2 0 25 x 78 0 050 10 150 20 250 3 0 05 0 1 0 15 0 2 0 25 x 上面的 表示引用刚刚结束的操作结果 因此对应的命令必须紧邻 如果隔开了就会出错 为了提高用 Mathematica 解题的灵活性 我们可以给操作结果起一个名 以后可以在任何时候 按名引用 例如 上面的命令组可以改写为 Solution NDSolve D u x t t D u x t x x u x 0 x 1 x u 0 t 0 u 1 t 0 u x 0 1 t 0 0 3 Table Plot Evaluate u x t First Solution t 0 0 3 PlotRange 0 0 25 x 0 1 0 9 0 4 2 其它热传导定解问题 下面 我们再考虑几个热传导方程的定解问题 例 1 在时间范围 0 1 t 内数值求解第二类边界条件齐次泛定方程问题 0 2 0 0 0 0 sin txx xxxx t uux uu ux 命令语句为 solution NDSolveA9 tu x tD x xu x tD u x 0D Sin xD2 uH1 0L 0 tD 0 uH1 0L tD 0 u 8x 0 8t 0 1 E Plot3D Evaluate u x tD First solutionDD 8x 0 8t 0 1 D 输出的结果为 0 1 2 3 0 0 2 0 4 0 6 0 8 1 0 0 25 0 5 0 75 1 0 1 2 79 例 2 在时间范围 0 1 t 内数值求解第一类边界条件非齐次泛定方程问题 1 2 0 2 0 sin2 0 0 0 sin txx xx t uutx uu ux 命令语句为 solution NDSolveA9 tu x tD x xu x tD Sin 2 tD 2 u x 0D Sin xD2 u 0 tD 0 u tD 0 u 8x 0 8t 0 1 E Plot3D Evaluate u x tD First solutionDD 8x 0 8t 0 1 D 输出的结果为 0 1 2 3 0 0 2 0 4 0 6 0 8 1 0 0 25 0 5 0 75 1 0 1 2 3 波动问题中的应用 同样地 我们可以 Mathematica 的数值求解微分方程的 DSolve 命令来求解波动方程的定解 问题 例如 在时间范围 0 1 t 内数值求解第一类边界条件波动问题 0 00 0 0 0 sin2 0 ttxx xx ttt uux uu ux u 的命令语句组为 solution NDSolveA9 t tu x tD x xu x tD u x 0D Sin 2 xD uH0 1L x 0D 0 u 0 tD 0 u tD 0 u 8x 0 8t 0 1 E Plot3D Evaluate u x tD First solutionDD 8x 0 8t 0 1 D 输出的结果为 80 0 1 2 3 0 0 2 0 4 0 6 0 8 1 1 0 5 0 0 5 1 0 1 2 又如 在时间范围 0 1 t 内数值求解第二类边界条件下的波动问题 0 00 0 0 0 cos2 0 ttxx xxxx ttt uux uu ux u 的命令语句组为 solution NDSolveA9 t tu x tD x xu x tD u x 0D Cos 2 xD uH0 1L x 0D 0 uH1 0L 0 tD 0 uH1 0L tD 0 u 8x 0 8t 0 2 E Plot3D Evaluate u x tD First solutionDD 8x 0 8t 0 2 D 输出的结果为 0 1 2 3 0 0 5 1 1 5 2 1 0 5 0 0 5 1 0 1 2 再如 在时间范围 0 6 t 内数值求解周期性边界条件下的波动问题 81 2 66 00 66 0 ttxx xx x ttt uux uu ueu 的命令语句组为 solution NDSolve D u t x t t D u t x x x u 0 x Exp x 2 Derivative 1 0 u 0 x 0 u t 6 u t 6 u t 0 6 x 6 6 Plot3D Evaluate u t x First solution t 0 6 x 6 6 结果得到 0 2 4 6 5 2 5 0 2 5 5 0 0 25 0 5 0 75 1 0 2 4 如果对图像的精度不满意 可以通过改变默认的插值点数来提过精度 例如 上面问题中 可以插入参数 PlotPoints 50 NDSolve D u t x t t D u t x x x u 0 x Exp x 2 Derivative 1 0 u 0 x 0 u t 6 u t 6 u t 0 6 x 6 6 Plot3D Evaluate u t x First t 0 6 x 6 6 PlotPoints 50 得到的结果成为 0 2 4 6 5 2 5 0 2 5 5 0 0 25 0 5 0 75 1 0 2 4 82 Mathematica 的 DSolve 命令还可以用来求解非齐次的波动问题 例如 在时间范围 0 1 t 内数值求解定解问题 0 22 00 0 1 0 1 0 ttxx xx ttt uux uu uxu 的命令语句组为 solution NDSolveA9 t tu x tD x xu x tD u x 0D 1 x2 2 uH0 1L x 0D 0 u 0 tD 1 u tD 0 u 8x 0 8t 0 1 E Plot3D Evaluate u x tD First solutionDD 8x 0 8t 0 1 D 输出的结果为 0 1 2 3 0 0 25 0 5 0 75 1 0 0 25 0 5 0 75 1 0 1 2 4 其它问题 NDSolve 命令也可以用来求解稳定问题 例如 在空间范围 0 1 y 内数值求解求定解问 题 01 00 0 01 0 1 1 cos2 0 xxyy xx yyy uuxy uu ux u 的命令语句组为 solution u FirstANDSolveA9 x xu x yD y yu x yD 0 u x 0D Cos 2 xD uH0 1L x 0D 0 u 0 yD 1 u 1 yD 1 u 8x 0 1 8y 0 1 EE Plot3D solution x yD 8x 0 1 8y 0 180 结果得到 0 0 25 0 5 0 75 1 0 0 25 0 5 0 75 1 2 109 0 2 109 0 0 25 0 5 0 75 NDSolve 命令还可以用来求解非线性问题 例如 求定解问题 01 0 01 sin2 tx xx t uuux uu ux 的命令语句组为 solution u First NDSolve 8 tu x tD u x tD xu x tD u x 0D Sin 2 xD u 0 tD u 1 tD 8u 8x 0 1 8t 0 5 DD Plot3D solution x tD 8x 0 1 8t 0 5 D 得到的结果为 84 0 0 2 0 4 0 6 0 8 1 0 0 1 0 2 0 3 0 4 0 5 1 0 1 0 0 2 0 4 0 6 0 8 又如 求非线性定解问题 2 2 1010 00 1 12 1010 0 ttxx xx x ttt uuuux uu ueu 的命令语句组为 solution NDSolve D u t x t t D u t x x x 1 u t x 2 1 2u t x u 0 x Exp x 2 Derivative 1 0 u 0 x 0 u t 10 u t 10 u t 0 10 x 10 10 输出的结果为 u InterpolatingFunction 0 10 10 10 追加一个绘图命令 Plot3D Evaluate u t x First solution t 0 10 x 10 10 PlotPoints 80 得到输出结果为 0 2 4 6 8 10 10 5 0 5 10 1 0 1 0 2 4 6 8 我们也可以用下面的命令来绘制一张平面的灰度图 DensityPlot Evaluate u 10 t x First solution x 10 10 t 0 10 PlotPoints 200 Mesh False 85 结果得到 10 50510 0 2 4 6 8 10 7 3 拉普拉斯变换方法 1 典型案例 利用拉普拉斯变换也可以求解偏微分方程 例如 下面的热传导问题 04 0 04 0 0 6sin 2 3sin txx xx t uux uu uxx 对时间变量 t 作拉普拉斯变换 得到 6sin 2 3sin 04 0 0 4 0 sU xxxUxx UU 其中 U x为 u x t的拉普拉斯变换 由上式解出 U x后 再作拉普拉斯变换的逆变换 得到 u x t 这个过程的 Mathematica 语句如下 eq tu x tD x xu x tD e1 LaplaceTransform eq t sD e2 e1 9u x 0D 6 Sin x 2D 3 Sin xD LaplaceTransform u x tD t sD U xD LaplaceTransformAuH2 0L x tD t sE D U xD 8x 2 D ans DSol