解偏微分方程

第第第六六六章章章 解解解偏偏偏微微微分分分方方方程程程 1 差差差分分分法法法解解解偏偏偏微微微分分分方方方程程程 2 热热热传传传导导导方方方程程程 ut = a2uxx 的的的差差差分分分公公公式式式 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)+diag(ones(NPTS-1,1),1)+diag(ones(NPTS-1,1),-1); PHI0(1)=0;PHI0(NPTS)=0; for x=2:NPTS-1; PHI0(x)=exp(0.5i*x)*exp(-(x-x0)^2*log10(2)/sigma^2); end PHI(:,1)=PHI0.’; for y=2:time CHI(:,y)=4i*(A\PHI(:,y-1)); PHI(:,y)=CHI(:,y)-PHI(:,y-1); end mo=moviein(time); for j=1:time plot([1:NPTS],abs(PHI(:,j)).^2/norm(PHI(:,j)).^2,[1:NPTS],v/2); mo(:,j)=getframe; end movie(mo) 图 高斯波包贯穿势垒的过程 3 弦弦弦振振振动动动方方方程程程 utt = a2uxx 的的的差差差分分分公公公式式式 3.1 显显显式式式格格格式式式 差分形式是 u(x, t +∆t)− 2u(x, t) + u(x, t−∆t) (∆t)2 = a2 u(x +∆x, t)− 2u(x, t) + u(x−∆x, t) (∆x)2 令x = i∆x, t = j∆t,按中心差分的表达式，可将偏微分变为差分式 表示 ∂2u ∂t2 = ui,j+1 − 2ui,j + ui,j−1 (4t)2 ∂2u ∂x2 = ui+1,j − 2ui,j + ui−1,j (4x)2 得显式差分公式 ui,j+1 = c(ui+1,j + ui−1,j) + 2(1− c)ui,j − ui,j−1 其中 c = a2 (4t) 2 (4x)2 . 当c < 1时,解是稳定的，当c = 1时， 可得到正确的数值解。当c > 1时,解是不稳定的。其截断误差 为O(4x2,4t2)。同样,读者可以自行推导相应的隐式公式。 3.2 初初初始始始条条条件件件 设t = 0时对应j = 1，初始位移φ(x)和初始速度为ψ(x)离散化以后 得φi = φ(xi)和ψi = ψ(xi).则初始位移表示成 ui1 = φi 而初速度用中心差分可表示为 ∂ui,1 ∂t = ui,2 − ui,0 24t = ψi 得ui,0 = ui,2 − 2ψi4t, 把它代入差分方程的最后一项,得 ui,2 = 1 2 [c(ui+1,1 + ui−1,1) + 2(1− c)ui,1 + 2ψi4t] 3.3 例例例题题题 两两两端端端固固固定定定的的的弦弦弦振振振动动动 两端固定的弦, 初速为零，初位移是 u(x, 0) =  h 2/3 x, (0 ≤ x ≤ 2/3) h 1− x 1− 2/3, (2/3 < x ≤ 1) 作图所用程序如下,其中取c = 0.05, l = 1, h = 0.05.这里使用的方程 与初始条件表示方法与上一节相同. N=4000; c=0.05; x=linspace(0,1,420)’; u1(1:420)=0; u2(1:420)=0; u3(1:420)=0; u1(2:280)=0.05/279*(1:279)’; u1(281:419)=0.05/(419-281)*(419-(281:419)’); u2(2:419)=u1(2:419)+c/2*(u1(3:420)-2*u1(2:419)+u1(1:418)); h=plot(x,u1,’linewidth’,3); axis([0,1,-0.05,0.05]); set(h,’EraseMode’,’xor’,’MarkerSize’,18) for k=2:N set(h,’XData’,x,’YData’,u2) ; drawnow; u3(2:419)=2*u2(2:419)-u1(2:419)+c*(u2(3:420)... -2*u2(2:419)+u2(1:418)); u1=u2; u2=u3; end 4 差差差分分分法法法与与与松松松弛弛弛法法法解解解椭椭椭圆圆圆型型型方方方程程程uxx + uyy = 0 4.1 差差差分分分法法法显显显式式式与与与隐隐隐式式式公公公式式式 利用二阶导数的中心差分公式得到二维拉普拉斯方程 uxx+uyy = 0的 差分公式为 u(x +∆x, y)− 2u(x, y) + u(x−∆x, y) (∆x)2 + u(x, y +∆y)− 2u(x, y) + u(x, y −∆y) (∆y)2 = 0 记∆x = ∆y = h, x = ih, y = jh, (i, j = 0, 1, 2, · · · , n− 1),得显式 u(i, j) = 1 4 [u(i + 1, j) + u(i− 1, j) + u(i, j + 1) + u(i, j − 1)] 稳定条件是(∆y) 2 (∆x)2 ≤ 1 . 类似地推导得出隐式格式的公式为 u(i, j + 1)− 2u(i, j) + u(i, j − 1) = −(∆y) 2 (∆x)2 [u(i + 1, j + 1)− 2u(i, j + 1) + u(i− 1, j + 1)] 它的截断误差是O((∆x)2, (∆y)2),该式是恒稳定的. 例题 差分法解平面温度场,设正方形的一边温度为10度,其余各边 的温度为零度,求稳定的温度场。 按照上面的差分公式，任意设定内部各点的温度为零，并假定当 两次迭代计算的值误差小于0.01时，就停止计算。编出如下程序 u=zeros(100,100); u(100,:)=10; uold=u+1; unew=u; for k=1:500 if max(max(abs(u-uold)))>=0.01 unew(2:99,2:99)=0.25*(u(3:100,2:99)+u(1:98,2:99)+... u(2:99,3:100)+u(2:99,1:98)); uold=u; u=unew; end end surf(u) 计算结果得到的图形为 4.2 松松松弛弛弛法法法 如果在差分公式中，随时将上一步算得的格点上的值替代旧值，并 且每次算出的新值也替换成新值与旧值的“组合”，则得到下列松 弛法的计算公式,其中0 < ω < 2。这个公式可以用变分原理去证明。 u(i, j) = (1− ω)u(i, j) + ω 4 [u(i + 1, j) + u(i− 1, j) + u(i, j + 1) + u(i, j − 1)] 例题 松弛法计算平面温度场,已知定解问题如下： uxx + uyy = 0; u(0, y) = 0, u(a, y) = µ sin 3piy b ; u(x, 0) = 0, u(x, b) = µ sin 3pix a cos pix a 为了进行数值计算，取µ = 1, a = 3, b = 2. 编出的程序如下 omega=1.5; x=linspace(0,3,30); y=linspace(0,2,20); phi(:,30)=sin(3*pi/2*y)’; phi(20,:)=(sin(pi*x).*cos(pi/3*x)); for N=1:100 for i=2:19 for j=2:29 ph=(phi(i+1,j)+phi(i-1,j)+phi(i,j+1)+phi(i,j-1)); phi(i,j)=(1-omega)*phi(i,j)+0.25*omega*(ph); end end end colormap([0.5,0.5,0.5]); figure(2), surfc(phi) 所得的图形如下：