第六章解偏微分方程.PDF

第六章 解偏微分方程 1 差分法解偏微分方程 2 热传导方程 ut = a 2uxx 的差分公式 2.1 显式格式 热传导方程可以写成差分形式(右边取t时刻的值计算) u(x, t + ∆t)− u(x, t) ∆t ≈ a2u(x + ∆x, t)− 2u(x, t) + u(x−∆x, t) (∆x)2 即 u(x, t+∆t) ≈ u(x, t)+ ∆t (∆x)2 a2 [u(x + ∆x, t)− 2u(x, t) + u(x−∆x, t)] 令x = i4x, t = j4t, i, j = 0, 1, 2, · · ·n − 1, r = ∆t (∆x)2 a2, 上式可写 成显式差分公式 u(i, j + 1) = (1− 2r)u(i, j) + r[u(i + 1, j)− 2u(i, j) + u(i− 1, j)] 稳定条件为∆t ≤ (∆x) 2 2a2 ,截断误差为O((∆x)2, ∆t). 例 细杆传热问题 定解问题是 ut = a 2uxx u(0, t) = 0, u(l, t) = 0 u(x, t = 0) = ϕ(x) 其中0 ≤ x ≤ 20, a = 10且 ϕ(x) = { 1, (10 ≤ x ≤ 11) 0, (x 10, x 11) 根据上面的公式,编写如下程序 x=0:20; t=0:0.01:1; a2=10; r=a2*0.01; u=zeros(21,101); u(10:11,1)=1; for j=1:100 u(2:20,j+1)=(1-2*r)*u(2:20,j)+r*( u(1:19,j)+ u(3:21,j)); plot(u(:,j)); axis([0 21 0 1]); pause(0.1) end surf(u) 2.2 隐式格式 热传导方程还可以写成如下差分形式(右边取t +4t时刻的值计算) u(x, t + ∆t)− u(x, t) ∆t ≈ a2u(x + ∆x, t + ∆t)− 2u(x, t + ∆t) + u(x−∆x, t + ∆t) (∆x)2 引入与上面相同的足标,且令 Hui,j+1 = a 2ui+1,j+1 − 2ui,j+1 + ui−1,j+1 (∆x)2 则得到隐式公式为 ui,j+1 − ui,j ∆t = Hui,j+1 令 Hui,j = a 2ui+1,j − 2ui,j + ui−1,j (∆x)2 则显式公式写成 ui,j+1 − ui,j ∆t = Hui,j 将显式与隐式相加,得平均公式 ui,j+1 − ui,j ∆t = 1 2 Hui,j + 1 2 Hui,j+1 即是 ui,j+1 = 1 + 1 2 H∆t 1− 1 2 H∆t ui,j 这个公式对任何步长都是恒稳定的. 含时间的一维薛定谔方程 一维运动的粒子的含时间的薛定格方程可以化成如下形式的抛物 型问题(设方程右边的系数为1), ∂ϕ ∂t = i ∂2ϕ ∂x2 利用前面推导的平均公式上面定义的算符H ,得 ϕi,j+1 = 1 + i12H4t 1− i1 2 H4t ϕi,j 改写成下述形式 ϕi,j+1 =  2 1− i1 2 H4t − 1 ϕi,j ≡ χ− ϕi,j 中间函数χ由下式定义， (1− i1 2 H4t)χ = 2ϕi,j 即是 (2i + H4t)χ = 4 i ϕi,j 在使用MATLAB编程时,先计算出每一时刻的φi,j对应的χ, 程序中是 把χ的系数表示为A, 然后利用MATLAB的矩阵方程求解的功能求χ， 再用它去求下一时刻的φi,j+1. 程序在一个220点的格子上，定义解析 形式的位势(方形势阱或势垒，Gauss势阱或势垒，位势台阶或抛物 线势阱)，建立一个初始的Gauss波包或Lorentz波包，它具有规定的 平均位置、动量和空间宽度，并在保持格子的端点上ϕ为零的边界 条件下演化．用动画在每一时间步长上都显示概率密度|ϕ|2和位势， 以及波包在一指定点的左边和右边两部分每―部分的总概率和平均 位置． 程序使用的数据是，选位势为一个方形位垒，高度是0.18, 位于格 点编号为105和115 的位置。初始的波包是一高斯波包，其平均波数 和宽度为k0 = 0.6, σ = 10, 中心位置x0在编号为40的格点上。时间 步长为1，演化时间为130。动画演示了波包向势垒行进，在势垒上 发生透射和反射，最后分为透射波和反射波二部分向不同的方向运 动。下图给出了这个过程中的几个画面。 NPTS=220; sigma=10;k0=0.6; x0=40; time=130; v(NPTS)=0;v(105:115)=+0.18; A=diag(-2+2i+v)+di