matlab解偏微分方程
合集下载
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
第六章
1 2
2.1
解偏微分方程
差分法解偏微分方程 热传导方程 ut = a2 uxx 的 差 分 公 式
显式格式
热传导方程可以写成差分形式(右边取t时刻的值计算) u(x, t + ∆t) − u(x, t) u(x + ∆x, t) − 2u(x, t) + u(x − ∆x, t) ≈ a2 ∆t (∆x)2 ∆t 2 a [u(x + ∆x, t) − 2u(x, t) + u(x − ∆x, t)] 即 u(x, t+∆t) ≈ u(x, t)+ (∆x)2 ∆t 2 令x = i x, t = j t, i, j = 0, 1, 2, · · · n − 1, r = (∆ a , 上式可写 2 x) 成显式差分公式
ui,j +1 − ui,j = Hui,j +1 ∆t Hui,j = a2 ui+1,j − 2ui,j + ui−1,j (∆x)2
ui,j +1 − ui,j = Hui,j ∆t 将显式与隐式相加,得平均公式 ui,j +1 − ui,j 1 1 = Hui,j + Hui,j +1 ∆t 2 2
u(i, j + 1) = (1 − 2r)u(i, j ) + r[u(i + 1, j ) − 2u(i, j ) + u(i − 1, j )] (∆x)2 2 , 截断误差为 O ((∆ x ) , ∆t). 稳定条件为∆t ≤ 2a2 例 细杆传热问题 定解问题是 2 ut = a uxx u(0, t) = 0, u(l, t) = 0 u(x, t = 0) = ϕ(x) 其中0 ≤ x ≤ 20, a = 10且 1, (10 ≤ x ≤ 11) 0, (x < 10, x > 11) 根据上面的公式,编写如下程序 ϕ (x ) = 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));
得ui,0 = ui,2 − 2ψi
1 ui,2 = [c(ui+1,1 + ui−1,1) + 2(1 − c)ui,1 + 2ψi t] 2
3.3
例题 两端固定的弦振动
两端固定的弦, 初速为零，初位移是 h x, (0 ≤ x ≤ 2/3) 2 / 3 u(x, 0) = 1−x , (2/3 < x ≤ 1) h 1 − 2/3
图 高斯波包贯穿势垒的过程
3
3.1
弦振动方程 utt = a2 uxx 的 差 分 公 式
显式格式
差分形式是 u(x, t + ∆t) − 2u(x, t) + u(x, t − ∆t) (∆t)2 u(x + ∆x, t) − 2u(x, t) + u(x − ∆x, t) = a2 (∆x)2 令x = i∆x, t = j ∆t,按中心差分的表达式，可将偏微分变为差分式 表示 ∂ 2u ui,j +1 − 2ui,j + ui,j −1 = ∂t2 ( t)2 ∂ 2u ui+1,j − 2ui,j + ui−1,j = ∂x2 ( x )2 得显式差分公式 ui,j +1 = c(ui+1,j + ui−1,j ) + 2(1 − c)ui,j − ui,j −1 2 2 ( t) 其中 c = a . 当c < 1时,解 是 稳 定 的 ， 当c = 1时 ， ( x )2 可 得 到 正 确 的 数 值 解 。 当c > 1时,解 是 不 稳 定 的 。 其 截 断 误 差 为O( x2, t2)。同样,读者可以自行推导相应的隐式公式。
3.2 初始条件
设t = 0时对应j = 1，初始位移φ(x)和初始速度为ψ (x)离散化以后 得φi = φ(xi)和ψi = ψ (xi).则初始位移表示成 ui1 = φi 而初速度用中心差分可表示为 ∂ui,1 ui,2 − ui,0 = = ψi ∂t 2 t t, 把它代入差分方程的最后一项,得
即是
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)
4
4.1
差分法与松弛 法 解 椭 圆 型方 程 uxx + uyy = 0
差分法显式与隐式公式
利用二阶导数的中心差分公式得到二维拉普拉斯方程 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 ) + =0 (∆y )2 记∆x = ∆y = h, x = ih, y = jh, (i, j = 0, 1, 2, · · · , n − 1), 得显式 1 [u(i + 1, j ) + u(i − 1, j ) + u(i, j + 1) + u(i, j − 1)] 4 (∆y )2 稳定条件是 ≤ 1 . 类似地推导得出隐式格式的公式为 (∆x)2 u(i, j ) = (∆y )2 u(i, j + 1) − 2u(i, j ) + u(i, j − 1) = − (∆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 ) +
作图所用程序如下,其中取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
即是
1 1 + H ∆t 2 ui,j +1 = ui,j 1 1 − H ∆t 2 这个公式对任何步长都是恒稳定的. 含时间的一维薛定谔方程 一维运动的粒子的含时间的薛定格方程可以化成如下形式的抛物 型问题(设方程右边的系数为1), ∂ϕ ∂ 2ϕ =i 2 ∂t ∂x
利用前面推导的平均公式上面定义的算符H ,得 1 1+i H t 2 ϕi,j +1 = ϕi,j 1 1−i H t 2 改写成下述形式 2 − 1 ϕi,j ≡ χ − ϕi,j 1 1−i H t 2 中间函数χ由下式定义， ϕi,j +1 = 1 (1 − i H t)χ = 2ϕi,j 2 (2i + H t)χ = 4 i ϕi,j 在使用MATLAB编程时,先计算出每一时刻的φi,j 对应的χ, 程序中是 把χ的系数表示为A, 然后利用MATLAB的矩阵方程求解的功能求χ， 再用它去求下一时刻的φi,j +1. 程序在一个220点的格子上，定义解析 形式的位势(方形势阱或势垒，Gauss势阱或势垒，位势台阶或抛物 线势阱)，建立一个初始的Gauss波包或Lorentz波包，它具有规定的 平均位置、动量和空间宽度，并在保持格子的端点上ϕ为零的边界 条件下演化．用动画在每一时间步长上都显示概率密度|ϕ|2和位势， 以及波包在一指定点的左边和右边两部分每―部分的总概率和平均 位置． 程序使用的数据是，选位势为一个方形位垒，高度是0.18, 位于格 点编号为105和115 的位置。初始的波包是一高斯波包，其平均波数 和宽度为k 0 = 0.6, σ = 10, 中心位置x0在编号为40的格点上。时间 步长为1，演化时间为130。动画演示了波包向势垒行进，在势垒上 发生透射和反射，最后分为透射波和反射波二部分向不同的方向运 动。下图给出了这个过程中的几个画面。 NPTS=220; sigma=10;k0=0.6; x0=40; time=130;
plot(u(:,j)); end surf(u)
2.2 隐式格式
Biblioteka Baidu
axis([0 21 0 1]);
pause(0.1)
热传导方程还可以写成如下差分形式(右边取t +
t时刻的值计算)
u(x, t + ∆t) − u(x, t) u(x + ∆x, t + ∆t) − 2u(x, t + ∆t) + u(x − ∆x, t + ∆t) ≈ a2 ∆t (∆x)2 引入与上面相同的足标,且令 Hui,j +1 = a2 则得到隐式公式为 令 则显式公式写成 ui+1,j +1 − 2ui,j +1 + ui−1,j +1 (∆x)2
ω [u(i + 1, j ) + u(i − 1, j ) + u(i, j + 1) + u(i, j − 1)] 4 例题 松弛法计算平面温度场,已知定解问题如下： uxx + uyy = 0; 3πy u(0, y ) = 0, u(a, y ) = µ sin ; b u(x, 0) = 0, u(x, b) = µ sin 3πx cos πx a 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));
1 2
2.1
解偏微分方程
差分法解偏微分方程 热传导方程 ut = a2 uxx 的 差 分 公 式
显式格式
热传导方程可以写成差分形式(右边取t时刻的值计算) u(x, t + ∆t) − u(x, t) u(x + ∆x, t) − 2u(x, t) + u(x − ∆x, t) ≈ a2 ∆t (∆x)2 ∆t 2 a [u(x + ∆x, t) − 2u(x, t) + u(x − ∆x, t)] 即 u(x, t+∆t) ≈ u(x, t)+ (∆x)2 ∆t 2 令x = i x, t = j t, i, j = 0, 1, 2, · · · n − 1, r = (∆ a , 上式可写 2 x) 成显式差分公式
ui,j +1 − ui,j = Hui,j +1 ∆t Hui,j = a2 ui+1,j − 2ui,j + ui−1,j (∆x)2
ui,j +1 − ui,j = Hui,j ∆t 将显式与隐式相加,得平均公式 ui,j +1 − ui,j 1 1 = Hui,j + Hui,j +1 ∆t 2 2
u(i, j + 1) = (1 − 2r)u(i, j ) + r[u(i + 1, j ) − 2u(i, j ) + u(i − 1, j )] (∆x)2 2 , 截断误差为 O ((∆ x ) , ∆t). 稳定条件为∆t ≤ 2a2 例 细杆传热问题 定解问题是 2 ut = a uxx u(0, t) = 0, u(l, t) = 0 u(x, t = 0) = ϕ(x) 其中0 ≤ x ≤ 20, a = 10且 1, (10 ≤ x ≤ 11) 0, (x < 10, x > 11) 根据上面的公式,编写如下程序 ϕ (x ) = 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));
得ui,0 = ui,2 − 2ψi
1 ui,2 = [c(ui+1,1 + ui−1,1) + 2(1 − c)ui,1 + 2ψi t] 2
3.3
例题 两端固定的弦振动
两端固定的弦, 初速为零，初位移是 h x, (0 ≤ x ≤ 2/3) 2 / 3 u(x, 0) = 1−x , (2/3 < x ≤ 1) h 1 − 2/3
图 高斯波包贯穿势垒的过程
3
3.1
弦振动方程 utt = a2 uxx 的 差 分 公 式
显式格式
差分形式是 u(x, t + ∆t) − 2u(x, t) + u(x, t − ∆t) (∆t)2 u(x + ∆x, t) − 2u(x, t) + u(x − ∆x, t) = a2 (∆x)2 令x = i∆x, t = j ∆t,按中心差分的表达式，可将偏微分变为差分式 表示 ∂ 2u ui,j +1 − 2ui,j + ui,j −1 = ∂t2 ( t)2 ∂ 2u ui+1,j − 2ui,j + ui−1,j = ∂x2 ( x )2 得显式差分公式 ui,j +1 = c(ui+1,j + ui−1,j ) + 2(1 − c)ui,j − ui,j −1 2 2 ( t) 其中 c = a . 当c < 1时,解 是 稳 定 的 ， 当c = 1时 ， ( x )2 可 得 到 正 确 的 数 值 解 。 当c > 1时,解 是 不 稳 定 的 。 其 截 断 误 差 为O( x2, t2)。同样,读者可以自行推导相应的隐式公式。
3.2 初始条件
设t = 0时对应j = 1，初始位移φ(x)和初始速度为ψ (x)离散化以后 得φi = φ(xi)和ψi = ψ (xi).则初始位移表示成 ui1 = φi 而初速度用中心差分可表示为 ∂ui,1 ui,2 − ui,0 = = ψi ∂t 2 t t, 把它代入差分方程的最后一项,得
即是
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)
4
4.1
差分法与松弛 法 解 椭 圆 型方 程 uxx + uyy = 0
差分法显式与隐式公式
利用二阶导数的中心差分公式得到二维拉普拉斯方程 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 ) + =0 (∆y )2 记∆x = ∆y = h, x = ih, y = jh, (i, j = 0, 1, 2, · · · , n − 1), 得显式 1 [u(i + 1, j ) + u(i − 1, j ) + u(i, j + 1) + u(i, j − 1)] 4 (∆y )2 稳定条件是 ≤ 1 . 类似地推导得出隐式格式的公式为 (∆x)2 u(i, j ) = (∆y )2 u(i, j + 1) − 2u(i, j ) + u(i, j − 1) = − (∆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 ) +
作图所用程序如下,其中取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
即是
1 1 + H ∆t 2 ui,j +1 = ui,j 1 1 − H ∆t 2 这个公式对任何步长都是恒稳定的. 含时间的一维薛定谔方程 一维运动的粒子的含时间的薛定格方程可以化成如下形式的抛物 型问题(设方程右边的系数为1), ∂ϕ ∂ 2ϕ =i 2 ∂t ∂x
利用前面推导的平均公式上面定义的算符H ,得 1 1+i H t 2 ϕi,j +1 = ϕi,j 1 1−i H t 2 改写成下述形式 2 − 1 ϕi,j ≡ χ − ϕi,j 1 1−i H t 2 中间函数χ由下式定义， ϕi,j +1 = 1 (1 − i H t)χ = 2ϕi,j 2 (2i + H t)χ = 4 i ϕi,j 在使用MATLAB编程时,先计算出每一时刻的φi,j 对应的χ, 程序中是 把χ的系数表示为A, 然后利用MATLAB的矩阵方程求解的功能求χ， 再用它去求下一时刻的φi,j +1. 程序在一个220点的格子上，定义解析 形式的位势(方形势阱或势垒，Gauss势阱或势垒，位势台阶或抛物 线势阱)，建立一个初始的Gauss波包或Lorentz波包，它具有规定的 平均位置、动量和空间宽度，并在保持格子的端点上ϕ为零的边界 条件下演化．用动画在每一时间步长上都显示概率密度|ϕ|2和位势， 以及波包在一指定点的左边和右边两部分每―部分的总概率和平均 位置． 程序使用的数据是，选位势为一个方形位垒，高度是0.18, 位于格 点编号为105和115 的位置。初始的波包是一高斯波包，其平均波数 和宽度为k 0 = 0.6, σ = 10, 中心位置x0在编号为40的格点上。时间 步长为1，演化时间为130。动画演示了波包向势垒行进，在势垒上 发生透射和反射，最后分为透射波和反射波二部分向不同的方向运 动。下图给出了这个过程中的几个画面。 NPTS=220; sigma=10;k0=0.6; x0=40; time=130;
plot(u(:,j)); end surf(u)
2.2 隐式格式
Biblioteka Baidu
axis([0 21 0 1]);
pause(0.1)
热传导方程还可以写成如下差分形式(右边取t +
t时刻的值计算)
u(x, t + ∆t) − u(x, t) u(x + ∆x, t + ∆t) − 2u(x, t + ∆t) + u(x − ∆x, t + ∆t) ≈ a2 ∆t (∆x)2 引入与上面相同的足标,且令 Hui,j +1 = a2 则得到隐式公式为 令 则显式公式写成 ui+1,j +1 − 2ui,j +1 + ui−1,j +1 (∆x)2
ω [u(i + 1, j ) + u(i − 1, j ) + u(i, j + 1) + u(i, j − 1)] 4 例题 松弛法计算平面温度场,已知定解问题如下： uxx + uyy = 0; 3πy u(0, y ) = 0, u(a, y ) = µ sin ; b u(x, 0) = 0, u(x, b) = µ sin 3πx cos πx a 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));