《华中科技大学计算方法课件-第七章常微分方程边值问题数值解法.pdf》由会员分享，可在线阅读，更多相关《华中科技大学计算方法课件-第七章常微分方程边值问题数值解法.pdf（63页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、第 七 章 常 微 分 方 程 边 值 问 题 数 值 解 法 一 个 由 常 微 分 方 程 及 其 边 界 条 件 构 成 的 定 解 问 题 称 为 常 微 分 方 程 边 值 问 题，其 常 常 出 现 在 流 体 力 学、材 料 力 学、波 动 力 学、电 磁 学 及 现 代 控 制 理 论 等 学 科 中。一 般 而 言，常 微 分 方 程 边 值 问 题 的 解 析 求 解 是 非 常 困 难 的。为 了 解 决 这 类 实 际 问 题，人 们 通 常 采 用 数 值 解 法。本 节 将 主 要 介 绍 71维 常 微 分 方 程 边 值 问 题 f，（力）=/（力。（力），i
2、C（Q 2）,的 打 靶 法、有 限 差 分 法 及 有 限 元 法。7.1打 靶 法 从 形 式 上 看 边 值 问 题 与 初 值 问 题 之 间 的 差 别 仅 在 于 定 解 条 件 的 不 同，毫 无 疑 问 这 两 类 问 题 有 着 紧 密 的 关 联。事 实 上，基 于 同 一 微 分 方 程 的 边 值 问 题 与 初 值 问 题 有 着 共 同 的 通 解 形 式。因 此，借 助 于 初 值 问 题 数 值 算 法 来 构 造 边 值 问 题 的 数 值 算 法 是 一 个 很 直 接 的 思 路，而 前 者 有 很 多 有 效 的 数 值 算 法 可 供 选 择。打 靶
3、法 的 基 本 思 想 是 从 不 同 的 初 值 出 发 用 已 知 的 初 值 问 题 数 值 算 法 解 相 同 的 微 分 方 程 直 到 解 的 轨 迹“打”到 给 定 的 右 边 界 点。为 简 单 起 见，本 节 我 们 仅 考 虑 如 下 几 维 线 性 边 值 问 题 f yr=A(t)y+q(t),t G(Q,6),B a M+Bby(b)=7 7,()其 中 4 以 和 场 为 X T Z的 矩 阵 q和 为 几 维 向 量，而 所 用 方 法 也 可 拓 展 到 非 线 性 边 值 问 题。7.1.1简 单 打 靶 法 对 于 一 个 常 微 分 方 程 而 言，不 同
4、 的 初 值 条 件 对 应 于 不 同 的 解 的 轨 线，即 解 是 依 赖 于 初 始 条 件 的。令 沙 依。)为 方 程/=力,)满 足 初 值 条 件 g(Q)=。的 解，若 g(力；。)同 时 也 满 足 边 值 问 题(2)的 边 界 条 件，那 么 g(1;c)正 好 就 是 边 值 问 题(2)的 解。简 单 打 靶 法 要 做 的 工 作 正 是 如 何 寻 求 这 样 的 初 值 c,使 得 B Q C+Bby(b c)=q(3)根 据 常 微 分 方 程 理 论，我 们 知 道 线 性 边 值 问 题 的 解 可 表 示 为 y(t)=Y(t)c+v(t),t e(a
5、,6),(4)其 中 Y)，。分 别 为 基 解 矩 阵 及 特 解，为 一 参 数 向 量，且 满 足Y，=4 K 力 G(Q,b),Y(Q)=/;(5)v=A(t)v+q，O(Q)=a,(6)这 里/为 级 单 位 矩 阵，a G 史 通 常 可 取 为 零 向 量。通 过 计 算 初 值 问 题(5)-(6),我 们 可 以 得 到 边 值 问 题 的 解 的 表 达 式(4)。为 了 确 定 参 数 向 量 c,将(4)代 入 问 题 的 边 界 条 件 可 得 Qc=力，(7)其 中 Q：=Ba+B6y(5),rj：=rj-Bav(a)-Bhv(b).若 边 值 问 题(2)存 在
6、唯 一 解，则 Q 是 非 奇 异 的，因 此 c是 存 在 唯 一 的。当 a为 零 向 量 时，。正 好 为 边 值 问 题 的 初 值 条 件，在 程 序 中 算 出 c后，我 们 再 次 计 算 一 次 初 值 问 题 得 到 原 边 值 问 题 的 数 值 解。_ 简 单 打 靶 法 的 计 算 程 序 如 下，其 中 bas说 上 八 和 par税 分 别 表 示 函 数 4 1 汝 和 4“g+q(,Ba,6b与 况 a表 示 边 鼻 参 数，tspan为 求 解 区 间。算 法 7.1 简 单 打 靶 法 _t,Y=sh o o tin g(b a s is fu n,p a
7、rtifu n,B a,Bb,tspan,eta)n=Ien gth(B a);I=e y e(n);T=zero s(n);t,y=od e45(p a rtifu n,tspan,z e r o s(n,1)；eta 1=eta Bb*y(le n g th(t),:)for j=1:n t,y=ode45(b a sisfu n,tspan,1(:,j);T(:,j)=y(length(t),:)5;endQ=Ba+Bb*T;c=Q eta 1;t,Y=o d e 4 5(p a rtifu n,tspan,c);例 7.1试 分 别 取 X=1,7,15,30,应 用 简 单 打 靶
8、法 求 解 0,1 上 的 边 值 问 题/0 1 0/00 0 1 什 02A3 A2 2A/(7r3+A27r)sin(7rZ)+(2A3+2A7r2)COS(T T)/I 0 0/0 0 0/(e-2A+2e-A+3)/(e-A+2)0 0 0 7/(0)+1 0 0 7/(1)=00 0 0/0 1 0/(A(3-e-A)/(e-A+2)要 求 画 出 数 值 解 和 精 确 解 的 对 比 图，方 程 精 确 解 为 g=（，”,（力）丁，其 中=gA(i-l)+e2A(i1)+ext2+e-x+COS(7T力).解 分 别 取 1,7,15,30,直 接 应 用 简 单 打 靶 法
9、 7.1可 计 算 出 相 应 的 数 值 解。Z521.510.50-0.50 0.5 1t图 7.1a A=1.21.51n 0.50-0.5-10 0.2 0.4 0.6 0.8 1t图 7.1b A=7.20151050-50 0.2 0.4 0.6 0.8 1t图 7.1c A=15.43210.八 20 x 105n0 0.2 0.4 0.6 0.8 1t图 7.1d A=30.从 上 图 可 见，随 着 M直 的 增 大，打 靶 法 的 误 差 也 增 大，出 现 该 现 象 的 原 因 在 于：一 方 面，相 应 的 矩 阵 Q=&+石 万 越 来 越 病 态；另 一 方 面，
10、直 愈 大，相 应 初 值 问 题 的 刚 性 也 愈 强，从 而 导 致 数 值 解 在 人 取 相 对 大 的 值 时 会 出 现 较 大 误 差。7.1.2多 重 打 靶 法 如 果 初 值 问 题 的 积 分 区 间 较 长，简 单 打 靶 法 的 计 算 效 果 将 会 降 低。为 克 服 简 单 打 靶 法 的 这 一 缺 点，下 面 我 们 介 绍 多 重 打 靶 法，其 基 本 思 想 是 首 先 将 求 解 区 间 作 网 格 分 划：a=力 N_i=A然 后 在 每 个 小 区 间 比 T 句 上 构 造 微 分 方 程 的 逼 近 解，最 后 由 解 的 连 续 性 及
11、边 界 条 件 可 将 这 些 数 值 解 修 补 在 一 起 得 到 一 个 全 局 解。考 虑 线 性 问 题(2),在 区 间 期 _1,切 上，边 值 问 题 的 解 可 以 表 示 为 y(t;Ci-i)=匕 Q 1+4，。=1,2,，N,(8)其 中 匕(力)，3 分 别 为 基 解 矩 阵 及 特 解，其 满 足 Y：=4 1)匕，匕(力 I)=I,(9)3 9 i)=0.(10)计 算 初 值 问 题(9)-(10)得 区 间 比 _ i,句 上 边 值 问 题 的 解 的 表 达 式(8)。为 了 确 定 参 数 向 量 c：=(%c f,，C(1)T，将(8)代 入 问 题
12、 的 边 界 条 件 得 g 匕(幻 4 一 1=%(切，1 Z?/-1,(11)Baco+B I Y N CN-X=T)Bi)VN(b).(12)(11)-(12)可 写 成 如 下 线 性 系 统 Lc=r,(13)其 中/-H(力 i)I B(力 2)/L=.,.,YN-ID N-I)I B Q B I Y N Jr=(。1(右)2(力 2),，n-BwN(b).解 线 性 系 统(13)即 得 网 格 点 力 上 的 数 值 解 纳=Q。多 重 打 靶 法 求 解 线 性 边 值 问 题 的 计 算 程 序 如 下，其 中 输 入 量 表 示 意 义 与 简 单 打 靶 法 基 本 相
13、 同，差 别 仅 在 力 帆 es/i为 求 解 区 间 上 的 全 体 网 格 点。算 法 7.2 多 重 打 靶 法 _Y=m u ltish o o tin g(b a sisfu n,parti fu n,B a,Bb,tmesh,eta)n=length(B a);I=eye(n);T=zeros(n);N=length(tmesh)1;L=zeros(N*n);r=zero s(N*n,1);Y=zeros(N+1,n);for k=1:N t,y J=ode45(p a rtifu n,tm esh(k),tm esh(k+l),zeros(n,1);r(n*(k l)+l:n*
14、i)=y(le n g th(t)for j=1:n t,y J=ode45(b a sisfu n,tmesh(k),tmesh(k+1),1(:,j);T(:,j)=y(length(tendL(n*(k l)+l:n*k,n*(k 1)+1:n*k)=T;if k小 L(n*(k l)+l:n*k,n*k+l:n*(k+l)=I;endendL(n*(N-l)+l:n*N,1:n)=Ba;L(n*(N-1)+1:n*N,n*(N 1)+1:n*N)=Bb*T;v=r(n*(N 1)+1:n*N);r(n*(N-1)+1:n*N)=eta Bb*v;c=Lr;for i=1:NY(i,:)
15、=c(n*(i l)+l:n*i)endY(N+1,:)=(T*c(n*(N l)+l:n*N)+v)将 多 重 打 靶 法 应 用 于 例 7.1中 1=3 0的 情 形，可 计 算 出 相 应 的 数 值 解（见 图 7.2）。由 图 7.2可 见，在 计 算 精 度 方 面，多 重 打 靶 法 比 简 单 打 靶 法 具 有 明 显 的 优 势。21n 0.50-0.51.5-10 0.2 0.4 0.6 0.8 1t图 7.2例 7.1。=30）的 精 确 解（虚 线）及 其 多 重 打 靶 解（点 戈 IJ线）.7.2有 限 差 分 法 有 限 差 分 法 也 是 一 种 基 于 初
16、 值 问 题 数 值 方 法 的 边 值 问 题 方 法，但 其 不 同 于 打 靶 法，有 限 差 分 法 在 利 用 初 值 问 题 数 值 方 法 离 散 边 值 问 题 后，将 求 解 区 间 上 所 有 网 格 点 处 的 数 值 解 作 为 离 散 格 式 的 未 知 量，从 而 通 过 求 解 表 征 该 格 式 的 线 性 或 非 线 性 代 数 方 程 组 而 获 得 边 值 问 题 的 数 值 解。在 本 节，我 们 仍 然 考 虑 一 阶 微 分 方 程 系 统 的 边 值 问 题(1)-(2)。为 了 获 得 边 值 问 题 的 系 列 离 散 格 式，我 们 对 求
17、解 区 间 作 如 下 剖 分：n：Q=力 1 12 力 N_1 且 记 步 长 hi=tj 纳 为 g(切 的 逼 近 值。7.2.1中 点 法 与 梯 形 法 在 本 节，我 们 先 考 虑 简 单 的 有 限 差 分 格 式，即 中 点 法 和 梯 形 法。分 别 应 用 中 点 法 和 梯 形 法 到 非 线 性 边 值 问 题（1）-（2）,我 们 依 次 可 得 其 问 题 的 离 散 格 式 yt yi-i(友 一 1/2,1(%+Vi-1，Q=1,2,.,N),g(y0,yN)=0(14)产 t=如+/佐 1,统 一 1）（分=1,2,，N）,。（如,yN）=0.（15）上 述
18、 计 算 格 式 可 视 为 以 向 量 J=（端，端，弧）T为 未 知 量 的（N+1）维 非 线 性 方 程 组，其 可 用 Newton迭 代 法 求 解。为 叙 述 简 单 见，这 里 我 们 仅 给 出 梯 形 法 的 Newton迭 代 格 式，中 点 法 的 Newton迭 代 格 式 类 似 可 得。记 NhVi：=纳 统）+/侑 1,%1）,i=l,2,，NF：=（Mg 即 式,Nh曲 g g 加 丁）丁,则 中 点 法 的 Newton迭 代 格 式 为 卜 旷）E=-F 铲）,r+1=r+,m=0,l,2,.-,M（16）其 中 尸（。为 歹。的 Jacobi矩 阵，加
19、为 迭 代 次 数。若 引 入 记 号：e=（品，册,4 切=力 2 n&=&（鳄,第）,Bb=次 端 第）则 迭 代 格 式（16）可 进 一 步 写 成（幻 金+41一 1后 一 1=N h g.N&+及 菰=一。（婚，第）;、川+1=次+&.=0,，N.（17）在 格 式（17）的 每 步 迭 代 中，通 过 前 二 组 方 程 求 出&然 后 将 其 代 入 第 三 组 方 程 求 出 当 前 逼 近 值 非+1=（附 F（必+乎,，（或 甲）斗。其 计 算 程 序 如 下：算 法 7.3 中 点 方 法 yj=m id p o in t(rfun,prfun,gfun,pgfun 1
20、,pgfun2,tmesh,yO)N,n=size(yO);N=N 1;I=eye(n);L=zeros(n*(N+1);r=zeros(n*(N+1),1);tol=1.d 12;xi=o n e s(n*(N+l),1)；y=yO;while(norm(x i)tol)for j=1:Nhj=tmesh(j+1)tmesh(j);A=feval(prfun,tmesh(j)+hj/2,(y(j,:)+y(j+l,:)/2);L(n*(j l)+l:n*j,n*(j 1)+1:n*j)=1/hj*1 A/2;L(n*(j l)+l:n*j,n*j+l:n*(j+1)=1/hj*1 A/2;r
21、(n*(j-l)+l:n*j)=(y(j+1 y(j,:)/hj+feval(rfun,tmesh(j)+hj/2,(y(j,:)+y(j+1,:)/2);endL(n*N+l:n*(N+l),l:n)=f e v a l(p g fu n l,y(l,:)5,y(N+l,:)5);L(n*N+l:n*(N+l),n*N+l:n*(N+l)=fe v a l(p g fu n 2,y(1,:)5,y(N+lr(n*N+l:n*(N+1)=feval(gfun,y(1,:)5,y(N+1,:)*);x i=L r;for k=1:N+1y(k,:)=y(k,:)+x i(n*(k l)+l:n*
22、k);endend end上 述 程 序 中，输 入 量 r f u n 和 p r f u n 分 别 表 示 函 数 gfun，pgfunl，pgfun2 分 别 表 示 函 数 q(%o),呼 2迦 弗 tmesh 表 示 求 解 区 间 上 所 有 网 格 点 同 工 1,.工 N，yO表 示 网 格 点 上 的 初 始 迭 代 值。例 7.2试 取 步 长 九=1/(10*2Z)(i=0,1,2,3),分 别 应 用 中 点 法 和 梯 形 法 求 解 非 线 性 两 点 边 值 问 题 力 oQ,_一 O,),康 e?-=+)/-/O(1X要 求 给 出 数 值 解 的 全 局 误
23、 差 及 收 敛 阶，其 方 程 的 精 确 解 为 解 令 u(t)2 Incosh(1 1/2)6*/2cosh(8/4),0=1.5171646.则 其 边 值 问 题 可 转 化 为 一 阶 形 式 yr=O 力 IB a y 6+牖=0.取 步 长 九=1/（10*2，）（，=0,1,2,3）,初 始 迭 代 值 为 零，分 别 应 用 中 点 法 和 梯 形 法 求 解 该 边 值 问 题 可 计 算 出 其 数 值 解 的 全 局 误 差 及 收 敛 阶（见 表 7.1）,其 中 全 局 误 差 及 收 敛 阶 的 表 达 式 依 次 为/八、I 1 e r r(h)err(ri
24、):=max(七 一%,p logo-J oi；v1 k 小 82 167 Td/2)表 7 中 点 法 和 梯 形 法 解 例 7.2的 数 值 结 果.h中 点 法 梯 形 法 err(h)Perr(h)P1/10 6.3341e-4-5.8442e-4-1/20 1.5938e-4 1.9907 1.4632e-4 1.99791/40 3.9867e-5 1.9992 3.6594e-5 1.99941/80 9.9698e-6 1.9996 9.1516e-6 1.9995由 表 7.1可 见，中 点 法 和 梯 形 法 均 只 有 2阶 精 度，并 且 二 者 的 计 算 效 果
25、相 近。7.2.2 Runge-Kutta方 法 由 上 可 知，中 点 法 和 梯 形 法 是 二 种 低 精 度 方 法。为 获 得 高 精 度 的 边 值 问 题 方 法，我 们 考 虑 修 改 求 解 初 值 问 题 的 Runge-Kutta方 法，使 其 适 再 于 边 值 问 题-（2）。将 Runge-Kutta方 法 应 用 于 边 值 问 题 可 得 Nhfi3：=/o-f=0,l i N,1 j s,SNhVi：=%Vi-i-hi bjfij=0,1 i N g(y&yN)=0,)=1(18)这 里 玲=+吗。该 方 程 组 中 的 未 知 向 量 y：=（诵 J f,对
26、，届，不，我）fi：=（总 虑，忌）,可 通 过 Newton迭 代 法 获 得，其 计 算 格 式 为 F(=F句 ym+1=r+,m=(19)其 中 尸 为 9 3 的 Jacobi矩 阵，机 为 迭 代 次 数，且 尸：=（N/J。NhV，Nhf Nh 状、,Nh代,NhV g0,V N Y、人 r r-i r r r r r r iNdi=（Nhfl N 麻,，明 博）丁.若 进 一 步 引 入 记 号：=（品 源 卓，苏 8）7=（4，脸 T,Q=-此 人 sBa=g y：y那,为=心 以,或），A g）=弘 熄+儿 口,1=1（4（力）（。114/讥），%=：，用=f：T 4（友
27、s）J Q S1A（力 谢）则 迭 代 格 式（19）可 写 成 Q 1 S4(力 讥)、:,。=(瓦/,C LssA.(tis I(%=Q?:,1 i N,&-1 h,Da,=0,1 S，S N)&a+&切=-g(鳄,。口 一+1=y+&-0 2 TV,，承+1=印+%I i TV.(20)既 然(20)中 的 第 一 个 方 程 可 写 为 口=叱 一 十%1+叱 一 崂 力 1 2 N,因 此 将 其 代 入(20)中 的 第 二 个 方 程 可 得&-1=陌，i i N,其 中 八 口=1+hiDWVi,n=hiDWQi.从 而 迭 代 格 式(20)可 写 成&一-1=片，1 z A
28、,rBa&+1+Bb&N=-g(y/,y用,=r+6：,o z 7 v,)举+1=优+叱-1%1+用 一 七 力 l i N.在 格 式(21)的 每 次 迭 代 中，先 通 过 前 两 个 方 程 计 算 出&然 后 代 入 后 两 个 方 程 完 成 一 次 迭 代，其 计 算 程 序 见 算 法 7.4。在 该 算 法 中，输 入 量 基 本 与 算 法 7.3的 同 名 输 入 量 相 同，唯 一 不 同 的 是 这 里 多 了 一 个 输 入 量/0,它 表 示 右 端 函 数 在 整 个 区 间 上 级 点 处 的 初 始 迭 代 值(足，井)丁。算 法 7.4 Runge-Kut
29、ta 方 法 fu n ctio n y=R u n g e_ K u tta(rfu n,prfun,gfun,pgfun 1,pgfun2,tmesh,yO,fO)N,nJ=siz e(yO);N=N 1;I=eye(n);L=zeros(n*(N+1);r=z e ro s(n*(N+l),1);to l=l.d 12;xi=o n e s(n*(N+l),1);y=yO;f=fO;alpha=l/4,l/4-sq rt(3)/6;1/4+sq rt(3)/6,1/4;beta=1/2,1/2;rho=l/2 sq rt(3)/6,1/2+sq rt(3)/6 5;s=len gth(b
30、eta);w hile(norm(x i)to l)for i=1:Nh _i=tm esh(i+1)tmesh(i);for j=1:st _ i j=tm esh(j)+h _ i*r h o(j);y _ ij=y(i+h _i*kron(alpha(j,:),I)*f(s*n*(i l)+l:s*n*i);A _ ij=fe v a l(p rfu n,t_ ij,y _i j);V_i(n*(j-1)+1:n*j,:)=A_ij;W _i(n*(j-1)+1:n*j,:)=kron(alpha(j,:),A_ij);Q_ij=fe v a l(rfun,t_ ij,y _i j)f(
31、s*n*(i l)+n*(j 1)+.1:s*n*(i l)+n*j);Q _ i(n*(j-l)+l:n*j,l)=Q _ ij;endw u e y e(s*n)h 1i*w i)u u w(1)*V；p T W i(l)*Q i;U(二 i)HU i;p(：i;r-(n*(i l)+l:n*i。n*(i l)+=n*ill(I+hL.*kron(befa。I)*w i(l)*v i)L(n*(i L)+l:n*i.n.i+l I n A i+c v I-r o A i1)+1-n*i h l i*kron(befa,I)*w i(l)*Q i;endL(n*N+l:n*(N+l)，一:n
32、)ufeval(pgfunl“y(一：)，y(N+l“：)”)；L(n*N+l-n*z+l)，n*N+l:n*z+l)=feval(pgfun5,y(l，y(N+l:.)；r(n*N+1n*z+l)n feval(gfun“y(l,:)r y z+lJ)x i u L-r”for kul:Ny Q l T y.:)+xi(n*(kll)+=n*k)二 f(s*n*(kll)+l-s*n*knf(s*n*(kll)+l:s*n*k)+:.u(r:,k)*xi(n*(kll)+=n*k)+p(：.k)；endy(N+l.)u y(N+l.)+xi(n*N+ljn*(N+l);end end分 别
33、取 步 长 九=1/(10*2，)(，=0,1,2,3),应 用 二 级 RadauIA方 法 和 RadauIIA方 法 于 例 7.2,可 获 得 其 数 值 解。在 表 7.2中，我 们 给 出 了 数 值 解 的 整 体 误 差 e rr=maxQiNI(切-Ui及 收 敛 阶 p=log2err(/z)/err(/z/2).由 表 7.2可 见，随 着 网 格 步 长 的 减 小，其 精 度 逐 步 提 高，且 二 种 方 法 均 约 有 3阶 精 度。表 7.2二 级 Radau IA、IIA方 法 解 例 7.2的 数 值 结 果.h2-stage Radau IA 2-stag
34、e Radau IIAerr(h)Perr(h)P1/10 1.8488e-5-1.6173e-5-1/20 2.3025e-6 3.0053 2.0038e-6 3.01281/40 2.8776e-7 3.0002 2.4989e-7 3.00341/80 3.6426e-8 2.9819 3.1646e-8 2.98127.2.3线 性 多 步 法 本 节，我 们 考 虑 修 改 求 解 初 值 问 题 的 线 性 多 步 法，使 其 可 应 用 于 边 值 问 题(1)-(2)。边 值 问 题 线 性 多 步 法 的 基 本 思 想 是：首 先 在 求 解 区 间 二 端 及 其 附
35、近 分 别 给 出 岛 个 和 个 逼 近 值，然 后 在 求 解 区 间 中 部 的 每 个 网 格 点 友 处 用 线 性 多 步 法 离 散，将 关 于 逼 近 解 纳 的 方 程 组 与 边 界 方 程 g(go,UN)=0联 立，通 过 适 当 方 法 求 解 该 联 立 方 程 组 而 获 得 原 边 值 问 题(1)-(2)的 数 值 解。边 值 问 题 线 性 多 步 法 的 一 般 形 式 为 上 2 卜 2:j+kiVi+j=h:f(i+j,Vi+j),i=k i,.,N k2,(22)j=-ki j=ki其 中 A；1+k2=k,to=a,ti=t i-1+h,h=(b-
36、a)/N,y,为 y g的 逼 近。为 了 求 解 方 程 组(2 2),我 们 仍 需 要 k 个 逼 近 值 班 加，1，以-%+1，以，其 可 由 k-1 个 与(22)同 阶 的 线 性 多 步 公 式 k i k if y i+j=h f 吃 J g j,%+/),，=1,岛 一 1,(23)j=i j=ik i k i i+jl/N-k+i+j=h f(t N-k+i+j l/N-k+i+j，=N 一 卷+1,.,Nj=T j=-i(24)及 边 界 条 件 go,UN)=O(25)给 出。至 此，我 们 即 得 到 求 解 边 值 问 题(1)-(2)的 线 性 多 步 法(22
37、)-(25)。为 实 现 格 式(22)-(25)的 计 算，我 们 将(22)-(24)写 成 矩 阵 形 式 其 中 y=(湍,对，,成 咒 F=(f也,y o F f g 班 产,VN)T)T,Qg Q IQ。a。(s-fc2+l)Q l-Qka尸 2+1).a f f+i)矩 阵 B 与 4 有 相 同 形 式，仅 需 将 A 中 的 内,a?分 别 代 之 以 如 邛 即 可 获 得。联 立(26)与(25),并 用 Newton迭 代 法 解 之，即 可 获 得 彷 侑 问 题(1)-(2)的 数 值 解?工下 面，我 们 介 绍 几 个 常 用 的 边 值 问 题 线 性 多 步
38、 法，其 推 导 及 相 关 理 论 可 参 见 Brugnano和 Trigiante的 专 著。三 阶 广 义 向 后 微 分 公 式(GBDF)(12%+1+3%6%_1+yn-2)=hfi,z=2,3,.,A-1,?/3+6y2 3Vl 2yo)=hfi,|(11W-18W-i+9VN-2 2J N _3)=hfN.五 阶 广 义 Adams方 法(GAM).一 统 一 1=4(一 19-2+346力 _ L+456力-74，+i+H+2),=2,.N-2yi-yo=含(19/4+106/3 264/2+646/1+251/0),yN-i U N-2=悬(19加+346/T V-I+4
39、56力 v2 74册 3+H/T V-4),、yN VN I=含(251/N+646/-i 264%2+106/N-3 19/N一 4).四 阶 扩 展 的 第 二 类 梯 形 方 法(ETR2)(盍(%+1+9%2)=亨(力+力 i=2,，N-1,击(一 期+9y2+9 阴 17 次)=亨(3一+/o),点(N3%N 2-9UNI+17g可)=亨(3加 i+于 6 六 阶 方 法(TOM)4(ll%+i+27yL 27t/i-i 11%2)=(力+i+9力+9/i-i+力-2),i=2,.,7V 1,以 如=合(27乃 173力+482/3-798为+1427力+475/0),VN VN-I
40、=y40(27/-5 _ 173加 _4+482加-3 798加-2、+1427加-1+475人).线 性 多 步 法 求 解 非 线 性 边 值 问 题(1)-(2)的 计 算 程 序 见 算 法 7.5,其 中 前 七 个 输 入 量 的 意 义 与 算 法 7.3相 同。算 法 7.5 线 性 多 步 法 _fu n ctio n y=bvm(rfun,prfun,gfun,pgfun 1,pgfun2,tmesh,yO,A,B)N,n=siz e(yO);N=N 1;I=eye(n);h=tm esh(2)tmesh(1);L=zer o s(n*(N+l);J=L;r=z e r o
41、 s(n*(N+l),1);to l=1.d-16;xi=ones(n*(N+1),1);y=yO;w hile(norm(x i)to l)for j=1:N+1Y(n*(j-l)+l:n*j,l)=y(j;F(n*(j l)+l:n*j,:)=fe v a l(rfun,tmesh(j),y(j,:)5);J(n*(j-1)+1:n*j,n*(j l)+l:n*j)=fe v a l(prfun,tmesh(j),y(jendL(1:N*n,:)=kron(A,I)h*kron(B,I)*J;r(1:N*n)=kron(A,I)*Y+h*kron(B,I)*F;L(n*N+l:n*(N+l
42、),l:n)=f e v a l(p g f u n l,y(l,:)?,y(N+lL(n*N+l:n*(N+1),n*N+l:n*(N+l)=fe v a l(pgfun2,y(1,:)5,y(N+1,:)5);r(n*N+1:n*(N+1)=fe val(gfun,y(l,y(N+l x i=L r;for k=1:N+1y(k,:)=y(k,:)+x i(n*(k l)+l:n*k)endend end为 测 试 其 算 法 的 有 效 性，我 们 分 别 取 步 长 九=1/(5*2，)(，=0,1,2),依 次 应 用 上 述 四 种 计 算 格 式 于 例 7.2中 的 边 值 问
43、 题，其 数 值 解 的 全 局 误 差 err(K):=0111a(幻 Uj和 收 敛 阶 1 r err(h)一 P-=log2 777Lerr(/z/2)J列 于 表 7.3中。由 该 计 算 结 果 可 见，其 算 法 关 于 非 线 性 边 值 问 题 是 非 常 有 效 的。hGBDF-3 ETR2-4 GAM-5 TOM-6errhperr Qi)perr(h)perr(h)p1/5 6.3674e-4-3.1464e-5-1.3510e-5-7.5878e-6-1/10 9.7062e-5 2.7137 2.0233e-6 3.9589 4.9753e-7 4.7631 4.4
44、113e-8 7.42631/20 1.4396e-5 2.7532 1.2918e-7 3.9693 1.7827e-8 4.8026 9.7162e-10 5.5047表 7.3线 性 多 步 法 解 例 7.2的 数 值 结 果.7.2 Ritz-Galerkin方 法 源 于 力 学 中 的 变 分 原 理 有 着 重 要 的 理 论 和 实 践 价 值，变 分 原 理 主 要 有 极 小 位 能 原 理 和 虚 功 原 理。在 特 定 的 函 数 空 间 设 置 下，微 分 方 程 边 值 问 题 可 转 化 为 与 之 等 价 的 变 分 问 题，通 过 求 解 相 应 的 变 分
45、 问 题，我 们 可 获 得 原 微 分 方 程 边 值 问 题 的 解。其 中，基 于 极 小 位 能 原 理 产 生 的 方 法 称 为 R itz方 法，而 基 于 虚 功 原 理 产 生 的 方 法 则 称 为 Galerkin方 法，有 时 人 们 也 笼 统 地 称 这 二 种 方 法 为 Ritz-Galerkin 方 法。7.3.1抽 象 空 间 为 了 方 便 本 节 及 后 面 第 九 章 的 应 用，这 里 我 们 首 先 介 绍 一 些 抽 象 空 间 概 念 和 常 用 的 积 分 不 等 式。设 有 区 域 Q u 及 巴 9Q为 其 光 滑 边 界，c=为 其 闭
46、 包，则 我 们 可 给 出 如 下 抽 象 空 间 概 念 及 其 常 用 结 论：0 软。)空 间：在。上 无 穷 次 可 微，在 3 Q 的 某 邻 域 上 恒 为 零 的 全 体 函 数。C7(0)空 间：在。上 加 次 连 续 可 微，在 3 Q 的 某 邻 域 上 恒 为 零 的 全 体 函 数。U(0)空 间(1 p 00)：在。上 夕 次 绝 对 可 积 的 全 体 函 数。易 知,空 间 L0)=/:|/|LP(Q)00是 一 个 Banach空 间。当 p=2 时，L2(0)空 间 中 的 内 积 和 范 数 可 分 别 定 义 为此 外，在 LP(Q)空 间 中，下 列
47、不 等 式 成 立：Minkowski 本 等 式:设/,g L/?(Q),则 11/+g|iL?(o)-II/I|LP(Q)+Holder 不 等 式：设 f e 叫 0),g e V(Q),+1=1,贝 l j 方 GV(Q),且 ll/g|u(o)W-(出 IlglLg).当 夕=q=2 时，Holder不 等 式 即 为 Cauchy-Schwarz不 等 式 ll/g|iM(c)叫 0)。一 阶 广 义 导 数：设/C L 2(Q),若 存 在 g e L 2(Q),使 得 g(c)砂(为 加=一 j/(劣)袈 加，V 矽 e C 软 C),JQ J.c，xi其 中 冗=(的,叼,，
48、喝)丁，则 称 g是/关 于 的 一 阶 广 义 导 数，且 记 为 小/(力)。侬 阶 广 义 导 数:设/6 L2(Q),若 存 在 g 6 L2(Q),使 得 g(x)x)d(j=(一 1)f(x)dad(y,V 协 e C 1(Q),Jo Jo其 中=(1 1,比 1,，n)T,且 Q=(1,Q2,an)e r 同=%,a冲=.端 J=I 1则 称 q为/的|a|阶 广 义 导 数，且 记 为 显 然 广 义 导 数 有 别 于 经 典 导 数，前 者 描 述 的 是 全 局 性 质，而 后 者 描 述 的 是 局 部 性 质。二 者 之 间 有 如 下 关 系：若/e C 3(0),
49、则 f 的 广 义 导 数 与 经 典 导 数 相 同。例 7.3设。=(-1,1),求 函 数 fx=团 3 G Q)的 一 阶 广 义 导 数。解 V。e Cg(Q),由 分 部 积 分 可 得,I PO P1/(x)dx=/(x)dx+/(x)dxJ-l J-l Jo/0 pl pl=/叭 x)dx/叭 x)dx=g(x)3(x)d*,J-i Jo J-i其 中 0-J l，。/11-1,-1 rr 0.因 此，根 据 广 义 导 数 的 定 义 可 知：d1f=g,但 其 经 典 导 数 不 存 在。在 广 义 导 数 的 基 础 上，我 们 可 定 义 加 阶 Sobolev空 间
50、间 叫。)如 下:Hm(Q)=/:VQ:|Q|m,df e L2(Q),其 内 积 与 范 数 可 定 义 为：(/,g E(加/即 匕)近 加 n/iim=.J。am特 别，H(0)即 为 L 2 g)空 间。在 下 面，我 们 也 将 用 到 Sobolev空 间 印 九(0)的 如 下 子 空 间 邮(0)：=/e 那(。)：/)|一=0.7.3.2 Ritz方 法 考 虑 一 维 椭 圆 边 值 问 题(dLy：=-p xyx+q(x)y(x)=f(x),7 6。：=(Q,b),M=0,y(b)=0,(27)其 中 p,q e C 中),P(X)po：=mi。M1),q(i)N o(%