精通MATLAB科学计算（第3版）(王正林)14-3r.pdf

《精通MATLAB科学计算（第3版）(王正林)14-3r.pdf》由会员分享，可在线阅读，更多相关《精通MATLAB科学计算（第3版）(王正林)14-3r.pdf（47页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、第 章 偏 微 分 方 程 求 解 科 学 与 工 程 中 的 大 量 现 象 与 规 律 都 可 以 用 偏 微 分 方 程 来 描 述，偏 微 分 方 程 的 定 解 问 题 成 了 表 述 科 学 与 工 程 中 的 现 象 和 规 律 的 重 要 数 学 工 具，应 用 极 其 广 泛。MATLAB提 供 了 简 单 易 用、功 能 强 大 的 工 具 箱 来 求 解 偏 微 分 方 程，同 时，还 可 以 通 过 MATLAB编 程 实 现 一 些 偏 微 分 方 程 的 求 解。通 过 本 章 学 习，读 者 不 仅 能 掌 握 MATLAB的 偏 微 分 方 程 工 具 箱(PD

2、E Toolbox)的 使 用，而 且 还 能 掌 握 利 用 MATLAB对 椭 圆 型 偏 微 分 方 程、抛 物 线 型 偏 微 分 方 程 和 双 曲 线 型 偏 微 分 方 程 这 三 大 类 偏 微 分 方 程 进 行 数 值 求 解。14.1 偏 微 分 方 程 概 述 所 谓 偏 微 分 方 程 是 指 含 两 个 或 两 个 以 上 自 变 量 的 微 分 方 程。本 章 主 要 讨 论 含 自 变 量 X和)，的 二 阶 偏 微 分 方 程，它 可 表 示 成 下 式 的 形 式：.O2M N.、du.d2u du duA(x,y)+B(x,y)+C(x,y)=f(x,y,

3、u,)d x cxdy d y ox dy其 中，自 变 量 取 值 范 围 是：xo x Xf,yn y yf边 界 条 件 为：“(x,yo)=by0(x),M(x,j/)=by,(x),u(xQ,y)=bXa(y),w(x/,y)=%,(y),根 据 偏 微 分 方 程 中 的 8、/、。系 数 之 间 的 关 系，可 以 把 偏 微 分 方 程 分 为 3 类：(1)椭 圆 型 偏 微 分 方 程，它 满 足：S2-4JC0(2)抛 物 线 型 偏 微 分 方 程，它 满 足：B2-4/C=0(3)双 曲 线 型 偏 微 分 方 程，它 满 足-第 1 4 章 偏 微 分 方 程 求

4、解 S2-4 T 4C0这 三 种 方 程 分 别 对 应 稳 定 状 态、扩 散 状 态 和 振 荡 状 态。由 于 偏 微 分 方 程 的 解 析 解 很 难 获 得，本 章 主 要 讲 述 它 们 采 用 MATLAB进 行 的 数 值 求 解 法。14.2 椭 圆 型 偏 微 分 方 程 以 下 仅 以 一 类 特 殊 的 椭 圆 型 偏 微 分 方 程 Heknhotz方 程 为 例，介 绍 椭 圆 型 偏 微 分 方 程 的 求 解。1 4.2.1 常 规 Helmholtz方 程 的 数 值 解 Helmholtz方 程 是 椭 圆 型 偏 微 分 方 程 中 的 一 种，它 满

5、 足：2(x,y)+g(x,y)u(x,y)=鳖)+，号)+g(x,y)u(x,y)=f(x,y)dx dy其 自 变 量 取 值 范 围 D 为：X o x Xf,yo y y f边 界 条 件 为：u(x,y0)=by,(x)r“(x,力)=by f(x),u(x0,y)=bx(y)-u(xf,y)=by f(y),如 果 Helmholtz方 程 中 g(x j)=0,则 它 被 称 为 Possion方 程。如 果 同 时 满 足 g(x,y)=0 和/(x,y)=0,则 它 被 称 为 Laplace方 程。1.求 解 原 理 描 述 为 了 采 用 不 同 的 方 法 求 解，对

6、区 域 D 进 行 分 割：将 沏”x 町 范 围 内 的 工 轴 等 分 为 M x段，每 段 长 为：Ax=(x/-XQ)I Mx；同 理 将 歹 轴 分 割 成 必，段,每 段 长 为 丁=(力-N o)/%。在 Helmholtz方 程 中，利 用 三 点 中 心 差 分 估 计 二 阶 导 数：d2u(x,y),2w；;y+Uij-dx2 1 A?其 中：Xj=X o+JAx,%=Jo+Z A yd2u(x,y)UM J 际 乂 R 7 1dy Ay 355精 通 MATLAB科 学 计 算(第 2 忡-tii.j=uxj,yi)因 此，对 于 自 变 量 区 间 上 的 每 一 个

7、 内 点(勺,),可 以 得 到 如 下 等 式:ui,j+-+jj-l,uM,j 27+Ui-ij F Z p 8i-jUi-j=加 其 中：Ui,j=u(X j,yi),fiJ=g(X j,yi)在 区 域 D 上 共 有(朋；-1)(%-1)个 内 点，因 此 可 以 建 立(-1)(a-1)个 方 程。但 当 八 较 大 时，联 合 求 解 方 程 组 显 然 不 现 实，迭 代 法 可 以 解 决 这 个 问 题。迭 代 法 首 先 把 边 界 条 件 变 形 为 如 下 形 式：ui.j=a，(ij+l+M/J-l)+rr(M/+1.7+)+xy-fi.j)U i.o=bxo(y。

8、,；*=%/(g),u0 J=byo(Xi),uM y J=by f(xj)其 中：Ay2.Ax2 _ Ar2 Ay2.2(Ax2+Ay2)2(Ax2+A j2)2(Ar2+Ay2)迭 代 算 法 需 要 迭 代 初 值 的 经 验 知 识，采 用 边 界 值 的 平 均 数 作 初 值 具 有 较 好 的 效 果。2.求 解 程 序 MATLAB的 实 现 在 MATLAB中 编 程 实 现 的 Helmholtz方 程 求 解 程 序 为：Helmholtzo功 能：迭 代 法 求 解 Helmholtz方 程。调 用 格 式:u,x,y=Helmholtz g,bxO,bxf,byO,h

9、yf,D,Mx,My,MinErr,Maxlter)o其 中，/,g 为 函 数 名；fcd)为 在 边 界 x=xO上 的 函 数 值；应 为 在 边 界 广 上 的 函 数 值；勿 0 为 在 边 界 了=加 上 的 函 数 值；6M为 在 边 界 上 的 函 数 值；D 为 边 界 点 值，长 度 为 四 的 列 向 量,其 元 素 分 别 为 xO,0，川，犷；M r为 沿 x 轴 的 区 间 数；屿 为 沿 j 轴 的 区 间 数；M inErr为 误 差 控 制 因 子；M axlter为 最 大 迭 代 次 数；u,x,y:方 程“伉 初 在(x,y)点 的 函 数 值。356-

10、第 1 4 章 偏 微 分 方 程 求 解 迭 代 法 求 解 Helmholtz方 程 的 MATLAB程 序 代 码 如 下：function u,x,y=Helmholtz(f,g,bxO,bxfzbyO,byfz D,Mx,My,MinErr,Maxlter)通 方 程:u_xx+u_yy+g(x,y)u=f(xz y)%f,g:函 数 名；%bxO:在 边 界 x=xO上 的 函 数 值；%bxf:在 边 界 x=xf上 的 函 数 值；%byO:在 边 界 y=yO上 的 函 数 值；%byf:在 边 界 y=yf上 的 函 数 值；%D:边 界 点 值，长 度 为 四 的 列 向

11、 量,其 元 素 分 别 为 xO,xf,yO,yf；%Mx：沿 x 轴 的 区 间 数；%My:沿 y 轴 的 区 间 数；%MinErr：误 差 控 制 因 子；%MaxIter:最 大 迭 代 次 数；%u,x,y:方 程 u(x,y)在(x,y)点 的 函 数 值。xO=D(1);xf=D(2);yO=D(3);yf=D(4);dx=(xf-xO)/Mx;x=xO+0:Mx*dx;3构 造 内 点 数 组 dy=(yf-yO)/My;y=yO+0:My*dy;Mxl=Mx+1;Myl=My 4-1;%边 界 条 件 for m=1:Mylu(mz 1 Mxl)=bxO(y(m)bxf(

12、y(m);e左 右 边 界 endfor n=1:Mxlu(1 Myl,n)=byO(x(n);byf(x(n);当 上 下 边 界 end右 边 界 平 均 值 作 迭 代 初 值 sum_of_bv=sum(sum(u(2:My,1 Mxl)u(1 Myl,2:Mx)；u(2:My,2:Mx)=sum_of_bv/(2*(Mx+My-2);for i=1:Myfor j=1:MxF(i,j)=f(x(j)zy(i);G(i,j)=g(x(j),y(i);endenddx2=dx*dx;dy2=dy*dy;dxy2=2*(dx2+dy2);rx=dx2/dxy2;ry=dy2/dxy2;rx

13、y=rx*dy2;for itr=1:MaxIterfor j=2:Mxfor i=2:Myu(iz j)=ry*(u(i,j+l)+u(i,j-1)+rx*(u(i+lz j)+u(i-1,j)+rxy*(G(i,j)*u(i,j)-F(i,j);%迭 代 公 式 endendif itr 1&max(max(abs(u-uO)MinErr 茗 循 环 结 束 条 件 357精 通 MATLAB科 学 计 算(第 2 版 i-break;enduO=u;end3.实 例 分 析 下 面 给 出 求 解 Helmholtz方 程 以 及 其 特 例 Possion方 程 和 Laplace方

14、程 的 实 例。【例 14-1】迭 代 法 求 解 Helmholtz方 程 应 用 实 例。求 以 下 Helmholtz方 程 的 数 值 解：2 2dx2 dy2自 变 量 取 值：0 X,4 r O 八 4边 界：”(0,_y)=y 2 r(4,y)=16cos(y)r(x,0)=x 2 r(x,4)=16cos(x)解：可 知/(x,y)=x 2+/,g(x,y)=fx o在 MATLAB中 编 写 函 数 ex140l来 实 现 求 解，代 码 如 下。function exl401()带 调 用 函 数 Helmholtz,求 满 足 一 定 初 始 条 件 和 边 界 条 件

15、的 Helmholtz方 程 的 数 值 解 f=inline(,xA2+yA2 x y1);g=inline(*sqrt(x)x y1);xO=0;xf=4;yO=0;yf=4;告 自 变 量 取 值 范 围 Mx=30;My=30;年 等 分 段 数 bxO=inline(1yA2,1y);%边 界 条 件 bxf=inline(4A2*cos(y),T y*);byO=inline(xA2,1x);byf=inline(4A2*cos(x),1x);D=xO xf yO yf;Maxlter=400;MinErr=le-4;U,x,y=Helmholtz(f,g,bxO,bxf,byO,

16、byf,D,Mx,My,MinErr,Maxlter);elf,mesh(x,yzU);xlabel(1x)ylabel(1y)zlabel(1U)运 行 结 果 如 图 14-1所 示。358 第 14章 偏 微 分 方 程 求 解 图 1 4-1【例 14-1 的 Helmholtz方 程 的 数 值 解 若 g(x j)=0,演 变 为 Possion方 程。此 时，只 需 将 程 序 exl401.m中 的 代 码“g=inlineCsqrt(x)?x?y)”改 为“=旬 加(。,乂 厅);”即 可。可 得 其 数 值 解 的 图 形 如 图 14-2所 示。若 二 者 均 为 0,则

17、 演 变 为 L aplace方 程，即 将“g=inline(sqrt(x),x,y);改 为 g=inline(O,x,y)；将 f=i n l i n e(xA2+yA2 I x J y，)；改 为 f=i n l i n e(O,x,y);,可 得 其 数 值 解 如 图 14-3所 示。图 14-2 1例 14-1】的 Possion方 程 的 数 值 解 359精 通 M ATLAB科 学 计 算(第 2 蝌 图 1 4-3 1例 14-1 的 Laplace方 程 的 数 值 解 14.2.2 满 足 牛 顿 边 值 条 件 的 Helmholtz方 程 接 下 来 将 讨 论

18、具 有 牛 顿 边 值 条 件 的 Helmholtz方 程，所 谓 牛 顿 边 值 条 件 是 指：Su(x,y)_-x=xo-(y)dx1.求 解 原 理 描 述 对 于 满 足 牛 顿 边 值 条 件 的 Helmholtz方 程，将 左 边 界 微 分 方 程 用 差 分 方 程 估 计，可 以 得 到：*%。()，1.-1*M/.i-2&,(j,)Ax,z=1,2-A/V-lj2 Ar结 合 边 界 点，可 以 得 到：Mi.o=fy(u,j+U,-I)+rv(w,+i,o+Mf-I.o)+fxy(gi,0W),0-fi.O)=fy(W),l+-2砥)(%)Ax)+rx(M,+l,0

19、+H,-I.o)+rXy(g(,OW,.O fi,0)=2rrw,j+rx(Ui+i,o+T0)+%(gi,0%,0 fi,o 一 2砥(%)/Ax)j类 似 地，如 果 在 y=外 处 也 满 足 牛 顿 条 件，可 以 得 到：uoj-fy(woj+i+uoj-i)+2rxuj+rxy(go.jUoj foj 2hy(x,)/Ay),i=特 别 地，在 边 界 点(x0,%)处，可 以 得 到：oe=2(,oj+cM,o)+&j，(go.oo.o-%,0-23；(yo)/Ar+2S；“(xo)/Ay)i2.求 解 程 序 MATLAB实 现 在 MATLAB中 编 程 实 现 满 足 牛

20、顿 边 值 条 件 的 Helmholtz方 程 求 解 程 序 为：Helmholtz。功 能：求 解 满 足 牛 顿 边 值 条 件 的 Helmholtz方 程。360-第 1 4 章 偏 微 分 方 程 求 解 调 用 格 式：w,x,y=Helmholtz_Newton(4 g,dbx,bxO,bxf,byQ,byf,D,Mx,My,MinErr,M ax I ter)o其 中，7,g 为 函 数 名；流)x 为 牛 顿 边 界 函 数；bxO为 在 边 界 x=xO上 的 函 数 值；6切 为 在 边 界 x=切 上 的 函 数 值；byO为 在 边 界 了 可 0 上 的 函 数

21、 值；勿/为 在 边 界 产 切,上 的 函 数 值；。为 边 界 点 值，长 度 为 4 的 列 向 量，其 元 素 分 别 为 xO,犷/0,守；牍 为 沿 x 轴 的 区 间 数；明，为 沿 y 轴 的 区 间 数；M inErr为 误 差 控 制 因 子；M axlter为 最 大 迭 代 次 数；u,x,),为 方 程“应 仞 在(x,y)点 的 函 数 值。迭 代 法 求 解 满 足 牛 顿 边 值 条 件 的 Helmholtz方 程 的 MATLAB程 序 代 码 如 下：function uz xz y=Helmholtz_Newton(f,gz dbx,bxO,bxfzby

22、O,byf,D,Mx,My,MinErr,Maxlter)传 解 方 程：u_xx+u_yy+g(x,y)u=f(x,y)%dbx：牛 顿 边 界 函 数 学 其 他 条 件 同 Helmholtz函 数 xO=D(1);xf=D(2);yO=D(3);yf=D(4);dx=(xf-xO)/Mx;x=xO+0:Mx*dx;与 构 造 内 点 数 组 dy=(yf-yO)/My;y=yO+0:My1*dy;Mxl=Mx+1;Myl=My 4-1;for i=1:Myfor j=1:MxF(i,j)=f(x(j),y(i);G(iz j)=g(x(j),y(i);endendfor i=l:MyD

23、BX(1,i)=dbx(x(1),y(i);enddx2=dx*dx;dy2=dy*dy;dxy2=2*(dx2+dy2);rx=dx2/dxy2;ry=dy2/dxy2;rxy=rx*dy2;带 边 界 条 件 for m=1:Mylu(mz 1 Mxl)=bxO(y(m)bxf(y(m);宅 左 右 边 界 endfor n=1:Mxlu(1 Myl,n)=byO(x(n);byf(x(n);当 上 下 边 界 361精 通 MATLAB科 学 计 算(第 2 版 i-end,边 界 平 均 值 作 迭 代 初 值 sum_of_bv=sum(sum(u(2:Myr 1 Mxl)u(1 M

24、ylf 2:Mx),)；u(2:My,2:Mx)=sum_of_bv/(2*(Mx+My-2)for itr=1:MaxIterfor m=2:(Myl-1)轴 牛 顿 边 值 条 件 迭 代 u(mz1)=2*ry*u(m,2)+rx*(u(m+1,1)+u(m-1z1)+rxy*(G(mz1)*u(mz1)-2*dbx(x(l),y(m)/dx);endfor m=2:(Mxl-1)轴 牛 顿 边 值 条 件 迭 代 u(l,m)=ry*(u(l,m+l)+u(l,m-l)+2*rx*u(2,m)+rxy*(G(lrm)*u(lzm)-F(lzm)-2*dbx(x(m),y(1)/dy);

25、endfor j=2:Mxfor i=2:Myu(iz j)=ry*(u(iz j+1)+u(i,j-1)+rx*(u(i+lz j)+u(i-l,j)+rxy*(G(iz j)*u(i,j)-F(i,j);与 迭 代 公 式 endendif itr 1&max(max(abs(u-uO)MinErr 名 循 环 结 束 条 件 break;enduO=u;end3.实 例 分 析【例 14-2】迭 代 法 求 解 满 足 牛 顿 边 值 条 件 的 Helmholtz方 程 应 用 实 例。求【例 1 4-1 中 满 足 如 下 牛 顿 边 值 条 件 的 数 值 解：Su 2 Su 2

26、x=x,=y ly=x o=xox oy解：在 MATLAB中 编 写 函 数 ex 1402来 实 现 求 解：function exl402()超 调 用 Helmholtz_Newton解 满 足 一 定 边 值 条 件 的 方 程：f=inline(”人 2+y八 2，y);g=inline(1sqrt(x)x1f y1);xO=0;xf=4;yO=0;yf=4;Mx=30;My=30;dbx=inline(*xA2+yA2,1x,*y);bxO=inline(,yA2,y*);bxf=inline(14A2*cos(y),T y1);byO=inline(1xA21zx1);byf=

27、inline(14A2*cos(x),1x);u_xx+u_yy+g(x,y)u=f(x,y)当 构 造 函 数 f(x,y)=(xA2+yA2)多 构 造 函 数 g(.r,y)=1x超 自 变 量 取 值 范 围 等 分 段 数 多 牛 顿 边 值 务 边 界 条 件 362 第 14章 偏 微 分 方 程 求 解 D=xO xf yO yf;Maxlter=100;MinErr=le-4;U,x,y=Helmholtz_Newton(f,g,dbx,bxO,bxfzbyO,byfz D,Mx,My,MinErr,Maxlter);elf,mesh(x,yz U)xlabel(fx)yla

28、bel()zlabel(*U)运 行 结 果 如 图 14-4所 示。图 14-4 牛 顿 边 值 的 Helmholtz方 程 数 值 解 通 过 对 比 图 14-4与 图 14-1可 以 发 现，两 者 的 区 别 也 主 要 体 现 在 边 界 上，这 是 由 于 后 者 具 有 牛 顿 边 值 条 件 造 成 的。1 4.3 抛 物 线 偏 微 分 方 程 一 维 热 传 导 方 程 就 是 一 个 典 型 的 抛 物 线 偏 微 分 方 程，它 用 下 式 的 形 式 表 示：d2u(x,t)_ du(x,t)dx2=dt其 中：0 x”xf,0 t T此 外，还 需 给 出 边

29、界 条 件 和 初 值:=h0(t)u(xj,t)=hx f(t),H(X,O)=M O(X)本 节 讨 论 的 一 维 抛 物 线 型 方 程 的 数 值 解 法 有 显 式 前 向 欧 拉 法、隐 式 后 向 欧 拉 法、Grank-Nicholson法，在 本 节 最 后 将 介 绍 二 维 抛 物 线 的 求 解。1 4.3.1 显 式 前 向 欧 拉 法 363精 通 MATLAB科 学 计 算(第 2 忡-下 面 讨 论 利 用 显 式 前 向 欧 拉 法 解 抛 物 线 型 偏 微 分 方 程。1.原 理 描 述 对 于 抛 物 线 偏 微 分 方 程，仍 然 可 采 用 有 限

30、 项 的 差 分 方 法，因 此 对 时 间 轴 和 x 轴 进 行 划 分:将 0,盯 等 分 成 段，每 段 长 度 为 x=xf/M；将 时 间 轴 等 分 为 N 段，每 段 长 度 为 t=T/N o将 抛 物 线 偏 微 分 方 程,S2w(x,Z)du(x,t)A-=-dx2 dt的 左 边 用 中 心 差 分 估 计，右 边 用 前 向 差 分 估 计，就 可 以 得 到 下 面 的 差 分 方 程：一 2 夕 一 力 A-；-=-Ax2 A z采 用 显 式 前 向 欧 拉 算 法 对 上 式 进 行 迭 代 求 解，即 得 到 形 如 下 式 的 变 形：夕”=r(u+u-

31、)+(1 2F)W/其 中：r=A*,k=l,2 M-1Ax22.求 解 程 序 MATLAB实 现 在 MATLAB中 编 程 实 现 的 显 式 前 向 欧 拉 法 求 解 抛 物 线 方 程 的 程 序 为：EF_Eulero功 能：显 式 前 向 欧 拉 法 求 解 一 维 抛 物 线 方 程。调 用 格 式:w,x,t-EF_Euler(/,xf,T,itO,bxO,hxf,M,N)o其 中，A 为 抛 物 线 偏 微 分 方 程 系 数；必 为 x 取 值 下 限；7 为 f取 值 下 限；讯)为 在 边 界 仁 0 上 的 函 数 值；6x0为 在 边 界 x=0上 的 函 数

32、值；6切 为 在 边 界 x=./上 的 函 数 值；M 为 沿 x 轴 的 等 分 段 数；N 为 沿/轴 的 等 分 段 数；u/,4：方 程“(X,。在(X,/)点 的 函 数 值。显 示 前 向 欧 拉 法 求 解 抛 物 线 方 程 的 MATLAB程 序 代 码 如 下：function u,x,t=EF_Euler(A,xf,T,itO,bxO,bxf,M,N)364 第 14章 偏 微 分 方 程 求 解 集 解 方 程 A u_xx=u_t,0=x=xf,0=t 0.5)disp(1r0.5,unstability);endfor k=1:Nfor i=2:Mu(iz k+1

33、)=r*(u(i+1,k)+u(i-1,k)+rl*u(iz k);endend3.实 例 分 析【例 14-3】显 式 前 向 欧 拉 法 求 解 一 维 抛 物 线 方 程 应 用 实 例。求 满 足 下 列 条 件 的 热 传 导 方 程 的 数 值 解：A=0.5,0 x 2,0 t 0.1解：热 传 导 方 程 为 典 型 的 一 维 抛 物 线 偏 微 分 方 程，如 下 所 示:d2u(x,t)du(x,t).A-=-,(x,0)=sm 兀 x x(-00,00)dx dt代 入 未 知 参 量，并 考 虑 所 给 变 量 取 值 范 围。在 MATLAB中 编 写 函 数 ex

34、 1403来 实 现 求 解:function exl403()他 调 用 函 数 EF_Euler解 热 传 导 方 程 A=0.5;3方 程 系 数 itO=inline(sin(pi*x)r,x1);2初 始 条 件 bxO=inline(0 1);bxf=inline(1 0 1);3边 界 条 件 xf=2;M=50;T=0.1;N=100;ul,x,t=EF_Euler(A,xf,Tz itO,bxO,bxf,M,N);figure(1),elf,mesh(t,xz ul)365xlabel(fx)ylabel(fy)zlabel(U)title(0.5)M=50;ul,xz t=

35、EF_Euler(Az xf,Tz itO,bxO,bxfz M,N);figure(2),elf,mesh(tz x,ul)xlabel(1x*)ylabel(yf)zlabel(Uf)title(r0.5二 二 H:二 K:一 二 0 5 0.02y 0 0 x图 1 4-5 利 用 显 式 前 向 欧 拉 法 得 到 的 抛 物 线 方 程 数 值 解：r0.51：i：1 1!j1广，、r：i.-4-11 i-.：i.1r:.i i L i Xn,:小 二 y 0 0 x图 1 4-6 利 用 显 式 前 向 欧 拉 法 得 到 的 抛 物 线 方 程 数 值 解 7i0.08不 稳 定

36、 区 域:0.08稳 定 区 域 366 第 14章 偏 微 分 方 程 求 解 1 4.3.2隐 式 后 向 欧 拉 法 下 面 讨 论 利 用 隐 式 后 向 欧 拉 法 解 抛 物 线 型 偏 微 分 方 程。1.求 解 原 理 描 述 采 用 后 向 差 分 估 计 式 52W(X,Z)du(x,t)的 右 边，左 边 仍 然 由 中 心 差 分 估 计，可 以 得 到：，哈-2*+血”尸 A-=-At写 成 迭 代 式：-ruf-i+(I+2r)“(一=w/其 中：r=A，,i=-1Ax2若 遮，焉 的 值 满 足 Dirichlet型 边 界 条 件，则-2+(1+2力/-m 3=

37、尸 式 可 写 成 矩 阵 形 式：l+2r-r 0 0 o-L ku+ruQ1+2厂-r 0 0后 芯 T0 T 1+2厂-0 0 茜=“铝 0 0 0 l+2r-rkUM-2k-UM-20 0 0 0 1+2二 kWA/-1r _ i kUM-+下 面 讨 论 牛 顿 边 值 情 形，即 满 足：du.,.|.v=o=尻 ox用 差 分 估 计 微 分，可 以 得 到：等 式 一 忆 3+(1+2 rM-%=中 1中 取 4 0,得 至 I J：367精 通 MATLAB科 学 计 算(第 2 版 i-fU-+(1+2r)W Q-t*ti=UQ综 上 可 得：(1+2r)o 2ru=uo-

38、2rbo(k)Ax此 时，上 述 矩 阵 式 可 以 改 写 为：l+2r r 0 0 0 w o-1-2rbo(k)x-r l+2r-r 0 0 U 尸 0-r 1+2厂 0 02=0 0 0-1+2r-rk“M-2k-UM-20 0 0 0 l+2r_ kUM-W.w-l+m y上 述 两 个 矩 阵 式 方 程 都 是 三 对 角 矩 阵，且 对 角 占 优 保 证 了 算 法 的 收 敛 性，因 此 用 后 向 隐 式 欧 拉 法 解 这 类 问 题 具 有 较 高 效 率。2.求 解 程 序 的 MATLAB实 现 在 MATLAB中 编 程 实 现 的 隐 式 后 向 欧 拉 法

39、求 解 抛 物 线 方 程 的 程 序 为：IB_Euler。功 能：隐 式 后 向 欧 拉 法 求 解 一 维 抛 物 线 方 程。调 用 格 式:u,x,=IB_Euler(Z,xf,T,itO,bxO,bxf,M,N)a其 中，A 为 抛 物 线 偏 微 分 方 程 系 数；切 为 x 取 值 下 限；讯)为 在 边 界 U 0 上 的 函 数 值；为 在 边 界 x=0 上 的 函 数 值；W 为 在 边 界=犷 上 的 函 数 值；用 为 沿 x 轴 的 等 分 段 数；N 为 沿/轴 的 等 分 段 数；u,X,4为 方 程 u(x,/)在(XJ)点 的 函 数 值。隐 式 后 向

40、 欧 拉 法 求 解 抛 物 线 方 程 的 MATLAB程 序 代 码 如 下：function u,x,t=IB_Euler(A,xf,T,itO,bxO,bxf,M,N)义 解 方 程 Al u_xx=u_t,0=x=xf,0=t=T软 初 值：u(x,0)=itO(x)%边 界 条 件：u(0,t)=bxO(t),u(xf,t)=bxf(t)%M：x 轴 的 等 分 段 数%N：t 轴 的 等 分 段 数 dx=xf/M;x=0:M*dx;368 第 14章 偏 微 分 方 程 求 解 dt=T/N;t=0:N*dt;for i=1:M+1,u(i,1)=itO(x(i);endfor

41、 n=1:N+1,u(1 M+1,n)=bxO(t(n);bxf(t(n);endr=A*dt/dx/dx;r2=1+2*r;for i=1:M-1P(i,i)=r2;%构 造 三 对 角 矩 阵 if i 1P(i-l,i)=-r;P(i,i-1)=-r;endendfor k=2:N+1b=r*u(1,k);zeros(M-3,1);r*u(M+1,k)+u(2:M,k-1);u(2:M,k)=linsolve(P,b);end3.实 例 分 析【例 14-4】隐 式 后 向 欧 拉 法 求 解 一 维 抛 物 线 方 程 应 用 实 例。利 用 隐 式 后 向 欧 拉 法 求 满 足 下

42、 列 条 件 的 热 传 导 方 程 的 数 值 解：A=0.50 x 2_,0 t 0.1即 利 用 隐 式 后 向 欧 拉 法 求 解【例 l4-31o解：热 传 导 方 程 为 典 型 的 一 维 抛 物 线 偏 微 分 方 程，如 下 所 示：d2u(x,t)8u(x,t)A-=-,(x,0)=sin/rx X G(-C O,O O)dx dt代 入 未 知 参 量，并 考 虑 所 给 变 量 取 值 范 围。在 MATLAB中 编 写 函 数 ex 1404来 实 现 求 解：function exl404()%调 用 IB_Euler解 热 传 导 方 程 A=0.5;超 方 程

43、系 数 itO=inline(sin(pi*x)1,*x1);%初 始 条 件 bxO=inline(1 0 1);bxf=inline(0 1);%边 界 条 件 xf=2;M=80;T=0.1;N=100;ul,xz t=IB_Euler(A,xf,T,itO,bx0f bxf,M,N);mesh(t,x,ul)xlabel(1x)ylabel(1y)zlabel(P)运 行 结 果 如 图 14-7所 示。369精 通 MATLAB科 学 计 算(第 2 蝌 图 1 4-7 利 用 隐 式 后 向 欧 拉 法 得 到 的 抛 物 线 方 程 数 值 解 从 图 14-7可 以 看 出，无

44、 不 稳 定 的 r 区 间，因 此 与【例 14-31相 比，采 用 隐 式 后 向 欧 拉 法 具 有 更 好 的 鲁 棒 性。14.3.3 Grank-Nicholson 法 下 面 讨 论 利 用 Grank-Nicholson法 解 抛 物 线 型 偏 微 分 方 程。1.原 理 描 述 本 节 讨 论 对 隐 式 后 向 欧 拉 法 的 一 种 改 进。对 于 等 式：-A-_-Ax2 Z其 左 边 的 差 分 估 计 在 时 间 点 木 进 行；如 果 把 右 边 的 差 分 估 计 看 做 是 时 间 步 长 为 A/2的 中 心 差 分，则 估 计 在 点 A和 h l 的

45、中 点 进 行。等 式 左 右 的 估 计 点 不 在 同 一 个 时 刻。为 了 消 除 这 个 不 一 致 性，需 要 将 等 式 左 右 的 估 计 放 在 同 一 个 时 间 点 进 行。在 时 间 点 发 和 好 1对 等 式 两 边 进 行 差 分 估 计，可 以 得 到：A u 1-+纤 2心+*?(Z?+后 人 K-对 上 式 进 行 变 形 就 得 到 Grank-Nicholson迭 代 方 法：ru+2(1+*ru-=ru+2(1+ru 3其 中，r=A-Ax2如 果 方 程 在 刈/XM处 具 有 牛 顿 边 值 条 件，则 可 把 迭 代 式 写 成 三 对 角 矩

46、阵 形 式：370 第 1 4 章 偏 微 分 方 程 求 解 该 算 法 具 有 高 效 率 和 恒 稳 定 性，即 对 于 任 意/,算 法 均 收 敛。-2+2r-r0-r2+2厂 r0-r2+2厂 0000000Ui002-2r00r000 2+2 00-r2+2r_o-kU r(o+I W o)r 2-2 r r 0 0 2 00 r 2-2 r 0 0 00 0 0 2-2 r rkW.w-l+00 0 0 0 2-2 r_ k_2”必(攵+1)+”必(行).2.求 解 程 序 的 MATLAB实 现 在 M ATLAB中 编 程 实 现 的 Grank-Nicholson法 求

47、解 抛 物 线 方 程 的 程 序 为：Grank_Nicholson0功 能：Grank-Nicholson法 求 解 一 维 抛 物 线 方 程。调 用 格 式：,x,t=Grank_Nicholson(y4,xf,T,itO,bxO,bxf,M,N)。其 中，A 为 抛 物 线 偏 微 分 方 程 系 数；为 x 取 值 下 限；7 为 取 值 下 限；加)为 在 边 界 仁 0 上 的 函 数 值；6x0为 在 边 界 x=0上 的 函 数 值；次 为 在 边 界 x R 上 的 函 数 值；用 为 沿 x 轴 的 等 分 段 数；N 为 沿/轴 的 等 分 段 数；u,t：方 程 w

48、(x,f)在(X,。点 的 函 数 值。Grank-Nicholson法 求 解 抛 物 线 方 程 的 MATLAB程 序 代 码 如 下：function u,x,t=Grank_Nicholson(A,xf,Tz itO,bxOz bxfzM,N)售 解 方 程 A u_xx=u_t,0=x=xf,0=t=T%初 值：u(x,0)=itO(x)%边 界 条 件：u(0,t)=bxO(t),u(xf z t)=bxf(t)%M：x 轴 的 等 分 段 数%N：t 轴 的 等 分 段 数 371精 通 MATLAB科 学 计 算(第 2 版 i-dx=xf/M;x=0:M*dx;dt=T/N

49、;t=0:N*dt;for i=1:M+1,u(i,1)=itO(x(i);endfor n=1:N+1,u(1 M+1,n)=bxO(t(n);bxf(t(n);endr=A*dt/dx/dxrl=2*(1-r);r2=2*(1+r);for i=1:M-1P(i,i)=rl;Q(i,i)=r2;if i 1P(i-1,i)=-r;P(i,i-1)=-r;Q(i-l,i)=r;Q(i,i-1)=r;endendfor k=2:N+1b=Q*u(2:Mr k-1)+r*(u(1,k)+u(lr k-1);zeros(M-2,1);u(2:MZ k)=linsolve(P,b);end372 第

50、 14章 偏 微 分 方 程 求 解 3.实 例 分 析【例 14-5 Grank-Nicholson法 求 解 一 维 抛 物 线 方 程 应 用 实 例。利 用 Grank-Nicholson法 求 满 足 下 列 条 件 的 热 传 导 方 程 的 数 值 解：4=0.5,0”x”2,0”f”0.1即 利 用 Grank-Nicholson法 求 解【例 14-3】。解：热 传 导 方 程 为 典 型 的 一 维 抛 物 线 偏 微 分 方 程，如 下 所 示：d2u(x,t)du(x,t)A-=-,w(x,0)=sin 7vx x G(-O O,CO)。dxz 8t代 入 未 知 参