名字起的太长不太好吧
范文一：金属—绝缘层—半导体结构泊松方程的数值解法
, ,在 , , ,, , , ,, ,,, ，, ,, ，,,,， ,,， ,, ,,,,, , ,, ，, , ,, , , , , 绝 金属 一 缘 层 一 半导 体 结 构 泊松 方 程 的数 值 解 法 吴 应 前 四川 大 学) ( 提 要 详细讨论了金属—绝缘层一 的 半导体结构一维泊捂方程 的数值解法( 推导了非均句结点下, 由 牛 顿迭代 公式 (提 出 了边 界条件 的恰 当形式 ( 作为 具 体应 用的 实例 ， 热氧 化 ,,,结 构 的高 额 测 量曲 线进行 了数 值计 算， 获得 了 — 关 系及界 面 陷阱 密度 分布 ( 将这 种 方法 与 常 用 的准 静态 怯 作 了实验 比较 ，两 者 的结 果符 合 良好( 但前者 在 制样 ，测量 方面 都 比后 故 者 简 单得 多， 适台作 为工 艺监 测 的常规 手段 ( 数 美键 词:泊松 方程， 值 计算法 ， 绝 半 金属 — 缘层一 导体 结 构， 表面 势， 界面 陷阱 密度( 一 、 引 言 ,金 绝 在 ,， ( 属一 缘层一半导体) 半 物理中 导体表 面 的基 本电学性质，诸如表面 电场， 表 荷 电荷 密 度 和 表 面 空 间 电 容 等 都 是通
均 匀杂 质 分 布 下 的一 维 泊 松 方 程 而 得 出 舶 ( 但 在 过求 解
实 际 ,， 所 ,器 件 中， 有 的掺 杂 工 艺 无 论 是扩 散 还 是 离 子 注 入 ，都 使 衬 底 杂 质 非 均 匀 分 布 ( 即使 是 单 纯 的热 氧化 工 艺 ，杂 质 的 再 分 布 效 应 也 会 使 杂 质 分 布 不 均 匀( ,器件 的衬底杂质分布实 际上都是 非均 匀的 此时泊松 方程不存在解析解 ， 因此任
只能 用数 值法解(本文详 细推导 了非均匀结点下 的牛顿迭代公何 ,，
式 ，提 出了边界条件 的恰 当形式(这样，在杂 质分 布已知的条件下 就可 以从 泊橙方程算 出任意表面势 下村 底 进 内 的 电势 分 布 ， 而 算 出 有 关 的 表 面 电学 参 数 作 为具 体 应 用 的 例 子 ，对 热 氧 化 ,,,结 , 算 一 关 系以及 界面陷阱分布( 为 了进行 比较 ， 同 构 的高频 ,一 羽量曲 线， 出了 ～ ( 对 一 , 样 品作 了准静态实验 钡量 ，结果证明符台 蘸好 (因为准静态 法制备样 品很麻烦 ，漏 电 面 要求很 苛刻， 高频 ,一 测量对 漏 电不敏感，极易进行(所 以把高频 一 与数值计算 结合起来可 以成为氧 化工 艺的常规 监测手段( 非均 匀结 点 下泊 松 方程 的 解 ,结构来说 ， 对 ,， ,衬底杂质分布 既可用多 种实 验手段测量也可 以按工艺条件通过计 算 机模拟计 算出来(测杂质分布 的实验方 法很多 ，其 中比较常见 的是 准 静 态 ,一 法、 , ,, 脉 冲 ,, 法 和 二 次谐 波 法 ( 如 果 采 用 目前 国 外 最 通 用 、 国 内亦 开始 移 植 的 , ,, ,, ,, 年 ,月 , ,, 年 ,日牧耐 (,, ,月 ,,日收 蓟修 改稿 应 用 科 学 学 报 ,卷 ,, ,, , ,虹, ,,, ,, ,, ,,, ,)(,, ,,,, , , , , ,,,,, , ,,, , ,程 序 ， 那 么 只 需将 所 用 工 艺 参 数 输 入计 算机，也可 以获得 比较 精确的杂 质分布 ( 下 面着重讨论在杂质 分布 ?( 已知 如 的情 况 下 ， 何 由 泊松 方程 解 出任 意表 面 势 下 衬 底 内 的 电势 分布 ( 以绝缘层一 以 () 半导体界 面为 原点建立指 向 衬 底 内部 的一维坐标 ( ? ) ?。 一 代 以 代 则 () 表 处的净施 主杂质 浓度， ) 表 处 的无量纲 电势 ， 泊松 方程为 , , 一 —, : ?( 旦 〔 一 … , 〕 , 百 … () , 、 绝 界 啦其 中 ? 为 半 导 体
( 对 ) 电常 数 ， 为 本 征 载 流 子 浓 度 为 玻尔 兹 曼 常数 ( 先将方程 (), 离散 化(因半导体 中有可 动载 流子， 能把 任何 外加 电压 引起 的衬底 电势 所 电变 化 主 要 集 中在 表 面层 ， 以在 这 个 部 分 应 把 结 点取 得 比较 密 (越 往 衬 底 深 处 ， 势变 化 结 ,离越 慢 ， 点 可 越 取 越 稀 ( 这 样 就能 以足 够 少 的结 点 数 准 确地 将式 () 散 化 (换 言之 ， 与 在 非采 用 均匀 结 点 相 比， 同样 多 的 结 点 总数 下 ， 均 匀 结 点 的离 散 误 差要 小得 多 ( 由 于 这 个 原 因，我们 可 以根据 实 际情 况 的 需 要 ， 按 在 ,,剜 体 内足 够 深 的 范 围 内， 从 密 剜 稀 的顺 序 取, , ,如 ，,个 分离 的 绪 点 :,、 、 ，… …， “ 是 其 中 ‰ 和 ， 两 个 辅 助 结 点 ( 取位 于 一,的界 面 处 ， ， 在衬 底 内 足 够 深 的位 置 上 ( 任 意 两 个 相邻 结 点 ,与 ，之 间 的 距 离 用 代 表 ，郎 , 』, 罩, … ，, ) — 一 一 ( ，, ，, 则 按拉 格朗 日插 值公式可将任意 处 的无量 纲
附近三个 结点 上 的 值 表 示 出 来 一 )与 一 丝 电势 ()用
, , ， 因 此 〔 器 。 丧 一 ， 〕 , ) ( ,在 , 离散，按式 () 各主结点上对方程 () 得差分方程 一 ， ( )， , 机 鲁 ， ( ” ， 旭 ， ( , ) 中其 ， ), ,(如 ; , , 丢 ( “击 , , ， , , ， ， ，, 一 ) 为 我 ? 断 处 的净 施 主
浓 度 ( 为 了简便 ， 们 可 以使 用 下 列矩 阵及 向 量 记号 一 杂 质
, ， , 粤, ,， ， , 一 , , 粤, ,，
, ( , ‰ 一, , 一 ,， 、总 ,善 , ? ～, — ～ 属 屡导 结泊方 的值法 盘 啪缘 半 体 构 襁 程 数 解 ,, , ,吼 广 一, , ， ，( ,, 〕一 , , , , 二 。 , , , …, ,, ,, , , ( ,便 这样式 () 可写
成 (( ;( ，( (( ( ， ( (( , ( 一,譬, ，( ) ( 一 式 ( 为非线性方程，〕 ? ，, , , , , 可用 , , , 法 迭 代求 解 先 用 ( 表 示 余 向量 (然 后在 预 报 值 线性 化 附近 将 , ) ( 即 甩 ， )酱 , ; , () 来代替方程 , ) (其 中, ( 一, 一 , 唧 代表校正值 ( , 知 由式 ( ) 蔷一一 斋 一 , , , 〕 中 其 , , 一 ,, .
范文二：双调和泊松方程的新解法
Ξ
双调和泊松方程的新解法
王 玲
()广东水利电力职业技术学院, 广州, 510635
摘 要 本文用 Fo u r ie r 变换法求解了区域为带形区域和上半平面的双调和 Po isso n 方程的边值问题。
双调和方程 带形区域 上半平面 变换 关键词 Fo u r ie r
The So lut ion of Four ier Tran sf orm to B iharm on ic Equa t ion
W an g L in g
(), , 510635Guangdo n T ech n ica l Co llege o f W a te r R e so u rce s and E lec t r ic E ng inee r ingGuangzho u
A bstrac t T h is p ap e r p re sen t th e bo unda ry va lue p ro b lem so lved a s th e b ih a rm o n ic equa t io n o f upp e r p lane
.and st r ip dom a in w ith Fo u r ie r t ran sfo rm ing m e tho d
.Keyword s b ih a rm o n ic po isso n equa t io n upp e r p lane st r ip dom a in Fo u r ie r t ran sfo rm
1 带形区域的边值问题
考虑下面双调和 Po isso n 方程边值问题:
4 4 4 5 u 5 u 5 u 2 ()()+=f x , y 1 ?u + 2 + - ?
4 4 5 u 5 u ()((2 ) ) ? < x="">< +?0,="" -="" u="" x="" ,="" 0="u" x="" ,="" 1="0," |="" y="|" y="=" 0="1" 445x="" 5x="">
- ?u = w , 则有假设 解方程
() ()- ?w = f x , y ,3 - ? < x="">< +?,="" 0="">< y=""><1>
(() ) ()w x , 0= w x , 1, - ? < x="">< +="" 4="" 作偶延拓:="">
() f x , y , 0 < y="">< 1,="" 3="" ()="" f="" x="" ,="" y="()" f="" x="" ,="" -="" y="" ,="" -="" 1="">< y="">< 0="">
() w x , y , 0 < y="">< 1,="" 3="" ()="" w="" x="" ,="" y=";" ()="" w="" x="" ,="" -="" y="" ,="" -="" 1="">< y="">< 0="" 3="" 3="" ()="" ()="" w="" x="" ,="" 0="w" x="" ,="" 1="0,">
再作周期为 2 的延拓:
Ξ 姚正安教授推荐 收稿日期: 2005 年 12 月 11 日
第 2 期 双调和泊松方程的新解法 107
~3 3 ( ) () f x , y +n = f x , y , n ? z ; 2
~ 3 ~3 () () wx , y +2n = wx , y , n ? z , ~ 3 ~3 () () wx , y + 2n = wx , 2n + 1= 0, n ? z ,0 ~ ~3 3 ( () () wx , 2n = wx , 2n + 1= 0, n ? z () () 从而 3, 4可化为全空间上的 Po isso n 方程
~ 3 ~3 ()- ?w= f (x , y ) , - ? < x="">< +="" -="">< 5="" y="">< +="" ~="" ~3="" 3="" (x="" ,="" 2n="" )="w(x" ,="" 2n="" +="" 1)="0," -="">< x="">< +="" w?,="" n="" z="" ()6="">
() 对方程 5施行 Fo u r ie r 变换, 利用导数性质得:
? ? 1 ~~ 3 3 () () f ()wΚ=Κ7 2| Κ| ? ? ~ ~ ~~3 3 3 3 其中w(Κ) 和f (Κ) 分别是 w(x , y ) 和 f (x , y ) 的 Fo u r ie r 变换
2 2 2 () | Κ| = Κ1 + Κ, 对 7作 Fo u r ie r 逆变换得: 2
? ? ? 1 1 ~3 ~ 3 3 ~() f Κ=() () 3 f x , y wx , y = 22 Κ Κ ||||+ ? 1 + ? 1 1 () [ f Ν, Γln d Γ + d Ν= ? 2 2?- ??02Π n= - ? () ) (x - Ν+ y + Γ 1 1 () ()fd Γ]Ν, Γln 8 22 ?0 ) (() y - Γ x - Ν+
() () 于是 1- 2可化为下列问题
~ 3 () - ?u = wx , y ,?
1 1 3 ~() () u x , y = ln3 wx , y .2 2 2Π y x +
2 上半平面边界值条件问题 我们考虑下面双调和泊松方程 2 ?u = f (, y ) ()x ,9 < +="" y="" ε="" 0,="" -="">< x="">
) ()(?
解方程假设 - ?u = w , 则有
() ()- ?w = f x , y , -11 0 ? < x="">< +="" y="" ε="">
(()) ? < x="">< +="" w="" x="" ,="" 0="0." 12="" -="">
() () 作偶延拓, 可以把 9- 10化为全空间上的 Po isso n 方程:
3 3 ()- ?w = f (x , y ) , 13 y < +="" -="">< x="">< +="" ,="" -="">< 3="" (())="">< x="">< +="" w="" x="" ,="" 0="0," -="" 14="">
3 3 () () () 这里 f 和 w x , y 分别是 f x , y 和 w x , y 的偶延拓, 即
() f x , y , y > 0, 3 () , f x , y = () f x , - y , y < 0="">
()w x , y y > 0, 3 () w x , y = ()w x , - y y <>
108 数学理论与应用 第 26 卷
() 对方程 13施行 Fo u r ie r 变换, 利用导数性质得:
? ? 3 ()f Κ 3 ()() 15 w Κ= 2| Κ| ? ? 3 3 3 3 2 2 2 () () () () 其中w Κ和f Κ分别是 w x , y 和 w x , y 的 Fo u r ie r 变换, Κ = Κ1 + Κ, 作 ||2
f
Fo u r ie r 逆变换得 1 1 3 3 () w , y = lnx 3 ()f x , y 2 2 2Π x + y + ? + ? 1 1 3 () fΝ, Γ ln d Νd Γ = 22 ??2Π- ?- ? (Γ) () y - x - Ν+
() () 于是, 9- 10可化为下面问题:
~3 () ()- ?u = wx , y ,16 - ? < x="">< +0,="" y="" ε="">
(()) 17 u x , 0= 0, y = 0.
() () 同理, 问题 9- 10的解为
+ ? + ? 1 1 1 1 () u , y = lnx 3 () d Ν[ d Γ Ν, Γln f 2 2 2 2?- ? ?02 Π2Π () ) (x + y x - Ν+ y + Γ + ? 1 () d Γ]+ Ν, Γln f 2 2?0 () ) (x - Ν+ y - Γ
3 结语
在科学研究的许多领域, Fo u r ie r 变换对于问题的求解和简化非常有用, 在偏微分方程的 求解与理论研究中, Fo u r ie r 变换亦有着广泛的应用。
2 () 本文对文献 1中 Po isso n 方程: - ?u = f x , , x ? R 的 Fo u r ie r 变换解法进行了推广,
求解了区域为带形区域和上半平面双调和 Po isso n 方程的边值问题。
参考文献
1 严子谦, 王俊禹, 王光烈, 崔志勇 1 数学物理方程 [M ]1 吉林: 吉林大学出版社, 1991176- 781
廖玉麟 1 数学物理方程 []1 武汉: 华中理工大学出版社, 199211190- 2021 M 2
3 吕 涛, 石 济 民, 林 振 宝 1 区 域 分 解 算 法—— 偏 微 分 方 程 数 值 解 新 技 术 [M ] 1 北 京: 科 学 出 版 社, 199911691
4 1 156严镇军 1 数学物理方程 [M ]1 合肥: 中国科学技术大学出版社, 19891154-
范文三：方程根的数值解法
毕 业 论 文
题 目 方程根的数值解法
学生姓名 朱敏敏
专 业 数学与应用数学
2011年 2 月 28 日
1
方程根的数值解法
摘 要
对于线性方程根是很容易求解的，而对于非线性方程，由于方程的多样性，非
线性方程根具有复杂性。利用数值解法解决非线性方程的求根问题。针对非线性
方程根的几种常用数值解法，对方程根的解法进行简单的介绍和探讨。
关键词:非线性方程 增值寻根法 二分法 迭代法 牛顿法 割线法
Abstract
The root of linear equations is very easy to solve, and for nonlinear equations, the diversity of the equation, nonlinear equation with complex roots. Numerical Solution of Roots of nonlinear equations to solve the problem. Roots of nonlinear equations for the numerical solution of several commonly used on the Roots of a simple solution presented and discussed.
Key words: nonlinear equation value dichotomy Root Law Newton secant method iteration method
2
目 录
1增值寻根法..............................................4
1.1确定方程根的存在区间..........................................4
1.2确定方程有根的区间后，利用增值寻根法求方程根的近视值...........4
1.3用增值寻根法求值..............................................4 2二分法..................................................5
2.1确定方程根的存在区间..........................................5
2.2利用二分法求方程根的具体步骤..................................5
2.3用二分法求值...................................................5 3 迭代法..............................................6
3.1简单迭代法
3.2收敛定理.......................................................6
3.3用迭代法求值...................................................6
4 牛顿法..............................................7
4.1牛顿法及其收敛性...............................................7
4.2牛顿法的计算步骤...............................................7
4.3用牛顿法求方程根...............................................7 5 割线法..............................................8
5.1割线法的基本思想...............................................8
5.2用弦割法求值...................................................8
.小结.................................................8 参考文献..............................................9
3
1.增值寻根法
1.1确定方程根的存在区间
* 增值寻根法的基本思想是，非线性方程f(x)=0的根是x，从初值x开始，0按规定的一个初始步长h来增值。令x=x+h(n=0,1,2,.....)，同时计算f(x)。n+1nn+1在增值的计算过程中可能遇到三种情形:
*f(x)=0，此时x即为方程的根x。 (1)n+1n+1
(2)f(x)和f(x)同符号，这说明方程在[x，x]内无根。 n+1nnn+1
(3)f(x)和f(x)异号，即f(x)?f(x)<>
程根的存在区间。X或x均可视为根的近视值。 nn+1
1.2确定方程有根的区间后，利用增值寻根法求方程根的近视值
* 接下来就是要求出方程根x在该区间内更精确的近视值。为此，我们可以再次使用增值寻根法。例如选新步长h=h/10，这样得到的有根区间更小，从而得1*到更精确的x。若精度不够，可重复使用增值寻根法，直到有根区间的长度|xn+1
)或f(x)的值可近似认为是零，而x,x|<ε(ε为所要求的精度)为止。此时f(xnn+1nn+1或x就是满足精度的方程的近似根。如图所示:>
1.3用增值寻根法求值
3例1 用增值寻根法求方程f(x)=2x+3x,1=0的有根区间。
解:取x=-2,h=1,则计算结果如下表: 0
x -2 -1 0 1
f(x) -23 -6 -1 4
所以f(x)的有根区间为(0,1)。再取x=0,h=0.1，计算结果如下表: 0
x 0 0.1 0.2 0.3 0.4
f(x) -1 -0.698 -0.384 -0.045 0.328 有此可知，f(x)=0的有根区间精确为(0.3，0.4)。以此类推，可以得到更为精确的有根区间。
2.二分法
2.1确定方程根的存在区间
4
依据零点定理可知，若函数f(x)在闭区间[a,b]上连续，且f(a)与 f(b)异号(即f(a)× f(b)<>
*点ξ(a<>
二分法的基本思想就是逐步对分区间,a,b,，通过判断两端点函数值乘积的符号，进一步缩小有根区间，将有根区间的长度缩小到充分小，从而求出满足精度的根的近似值，如图所示:
2.2利用二分法求方程根的具体步骤
具体步骤如下:将区间[a,b]的中点记作x,即x=(a+b)/2，并计算f(x)。111*同时记(a，b)=(a,b),如果此时f(x)恰好等于0，则我们已经找到方程的根x=x。1111若f(x)?0，f(a)?f(x)<><>
在(a，b)里面，区间长度是(a，b)的一半。对新的有根区间(a，b)111122实施同样平分区间的方法，判断中点与两端点函数值乘积的符号求出(a，b)，33它包含在(a，b)中，长度是(a，b)的一半，(a，b)的四分之一。重222211
复上述做法，得:(a，b)?(a，b)?…?(a，b)?…，且。1122nn
* x?(a，b)，n=1,2,…，，。？nn
**即序列{x}以等比数列的速度收敛于x。此时我们可以取x作为x的近似值。 nn
2.3用二分法求值
2-3例2 用二分法求方程f(x)=x-x-1=0的正根，精确到10。
解:利用二分法计算结果如下表:
n a b x f(x) nnn
1 1 2 1.5 0.875
2 1 1.5 1.25 -0.296875
3 1.25 1.5 1.375 0.2246094
4 1.25 1.375 1.3125 -0.0515137
5 1.3125 1.375 1.34375 0.0826111
6 1.3125 1.34375 1.328125 0.014576
7 1.3125 1.328125 1.3203125 -0.0187106
5
8 1.3203125 1.328125 1.3242188 -0.0021278
9 1.3242188 1.328125 1.3261719 0.0062089
10 1.3242188 1.3261719 1.3251954 0.0020369 *-3 ?x-x?<(b-a)><10，迭代10次，近似根x=1.3251954即为n101010所求。>
3.迭代法
3.1简单迭代法
迭代法是一种逐步逼近的方法。将f(x)变成另一种等价形式x=φ(x)，设x*为方程的根，选取某一近视值x?[a,b]，则按递推关系x=φ(x)，k=0,1,?，0k+1k产生迭代序列{x}。这种方法称为简单迭代法。 k
3.2收敛定理
压缩不动点定理或压缩印象定理:
若迭代函数g(x)满足(1)g(x)?[a,b],x?[a,b];(2)0<><1,使x',x"?[a,b],，,,有|g(x')-g(x")|?l|x'-x"| 0*则1="" x="g(x)在[a,b]内有唯一解x;">
0 2 对x?[a,b]，由x=g(x)得到{x}[a,b],且; ,0k+1kk,
0 3 有误差估计:;
0* 4 若g'(x)存在，则。
3.3用迭代法求值
例3 用简单迭代法求的近似值。 2
解: 设x=,1，则x(x+2)=1 2
所以求的近似值转化为求转化为求方程x(x+2)=1的正根，列出等价方2
程:
，以为迭代函数，以x=0为初始近似值，得到迭代序列: 0
?
取作为,1的近似值，得=1+ ?1.4142157。. 22
下证序列?收敛于,1，只要证满足上述收2
敛定 理，即证g(x)=在某个区间上满足上述收敛定理的条件。
6
取区间为[0,]，对x?[0,]，g(x)=?[0,]， ,
又u、v ?[0,]， ,
|g(u),g(v)|=。
所以g(x)在[0,]上满足收敛定理。
则迭代法收敛。即1.4142157就是的近似值。
4牛顿法
4.1牛顿法及其收敛性
牛顿法是一种线性化方法，其基本思想是将非线性方程f(x)=0逐步归结为某种线性方程来求解。设已知方程f(x)=0有近似根x(假定f'(x)?0)，将函数f(x)kk在点x展开，有f(x)?f(x)+f'(x)(x,x),于是方程f(x)=0可近似的表示为kkkk
f(x)+f'(x)(x,x)=0。这是个线性方程，记其根为f(x)，则x的计算公式为kkkk+1k+1
(k=0,1,?)，这就是牛顿法。
牛顿迭代收敛定理:设f(x)在区间[a,b]上二阶导数存在，且f(a)f(b)<>
?0,f"(x)不变号，初值x?[a,b]且f(x)f"(x)>0，则牛顿迭代序列{x}收敛于方程000kf(x)=0，在[a,b]的唯一根α。
4.2牛顿法的计算步骤
步骤1 准备 选定初始近似值x,计算f=f(x)，f'=f'(x)。 00000
步骤2 迭代 按公式x=x,f/f'，迭代一次得新的近似值x，计算f=f(x)，1000111f'=f'(x)。 11
步骤3 控制 如果x满足，则终止迭代，以x作为所求的根;否11则转步骤4。这里是允许误差，而其中C是取绝对误差或相对误差的控制常数，一般可取C=1。. 步骤4 修改 如果迭代次数达到预先指定的次数，或者f'=0，则方法失败;否则1以(xf,f')代替(xf,f')转步骤2继续迭代。 1,110,00
4.3用牛顿法求方程根
x例4 用牛顿法解方程xe,1=0
解:这里牛顿公式为，取迭代初值x=0.5，迭代结果列于下0
7
表:
k x k
0 0.5
1 0.57102
2 0.56716
3 0.56714 -x题中所给方程实际上是方程x=e的等价形式，若有不动点迭代到同一精度要迭代17次，可见牛顿法德收敛速度是很快的。
5割线法
5.1割线法的基本思想
在牛顿迭代格式(k=0,1,?)中，将曲线y=f(x)上点(x，kf(x))的切线斜率f'(x)，改为其上两点连线的斜率 kk
则得双点或单点弦割法迭代格式
f(x)(x,x)kkk1,x,x, (k=1,2,?) k1k，f(x),f(x)kk1,f(x)(x,x)kk0x,x,k1k，f(x),f(x)k0
5.2用弦割法求值
x-4例5 用割线法求方程xe-1=0在0.5附近的根(ε=10)。
解:取x=0.5，x=0.6，由迭代公式求得下表 01
k X x-x 1kk-1
0 0.5 0
1 0.6 0
2 0.56754 0.03246
3 0.56715 -0.00036
4 0.56714 -0.00001
*故x=0.56714，满足精度要求。
小结
非线性方程求根的问题一直都是学习的重要一章内容。选择使用某种方法，应保证迭代收敛，收敛速度快，计算量少。实际计算时，最好先判断收敛性;但也不妨先试算几步，检查迭代序列是否收敛。对于收敛慢，甚至发散序列，可用加速技术加速收敛。以上介绍的五种方法，都比较简单，对于我们在学习非线性方程求根方面有很多的帮助。
8
参考文献
[1] 李乃成，邓建中.数值计算方法[M].西安:西安交通大学出版社.2002. [2] 魏毅强，张建国，张洪斌等编.数值计算方法[M].北京:科学出版社.2004. [3] 陈文娟.论方程根的数值解法.现代商贸工业.2009年第十三期.
[4] 李同胜，罗蕴玲.代数方程求根的一个新的数值解法.河北农业大学学报.2001年1月. [5] 郭晓梅.非线性方程数值解法探究.枣庄学院学报.2010年4月.
[6] 薛云峰，王惠书.关于对《数值分析》中用“割线法”求方程根的探讨.山西广播电视
大学学报.2005年2月.
[7] 黄清龙.同时求解非线性代数方程全部根的Newton迭代法。江苏工业学院学报.2004
年12月.
[8] 刘颖芬.一个经典方程的多种数值解法.赣南师范学院学报.2005年.
9
范文四：三维泊松方程数值模拟的多重网格方法
第２４卷第１期
２００９年２月（页码：１５４～１５８）
地球物理
ＩＮ
学进展
ｖ０１．２４，Ｎｏ．１
Ｆｅｂ．２００９
ＰＲｏＧＲＥＳＳＧＥＯＰＨＹＳＩＣＳ
鲁晶津，吴小平，ＫｌａｕｓＳｐｉｔｚｅｒ．三维泊松方程数值模拟的多重网格方法．地球物理学进展，２００９，２４（１）：１５４～１５８
Ｌｕ
ＪＪ，ＷｕＸＰ，ＫｌａｕｓＳｐｉｔｚｅｒ．Ｍｕｈｉｇｒｉｄｍｅｔｈｏｄｆｏｒ３ＤｍｏｄｅｌｉｎｇｏｆＰｏｉｓｓｏｎｅｑｕａｔｉｏｎ．ＰｒｏｇｒｅｓｓｉｎＧｅｏｐｈｙｓ．（ｉｎＣｈｉｎｅｓｅ），
２００９，２４（１）：１５４～１５８
三维泊松方程数值模拟的多重网格方法
鲁晶津１，吴小平¨，Ｋｌａｕｓ
（１．中国科学技术大学地球和空间科学学院，合肥２３００２６；
Ｓｐｉｔｚｅｒ２
２．ＩｎｓｔｉｔｕｔｆｔｉｒＧｅｏｐｈｙｓｉｋ，ＴＵＢｅｒｇａｋａｄｅｍｉｅＦｒｅｉｂｅｒｇ，Ｇｅｒｍａｎｙ）
摘要本文简要介绍了多重网格方法的基本思想和原理，然后应用多重网格（ＭＧ）方法求解三维泊松方程，网格尺度从１７×１７×１７逐次增加至２５７×２５７×２５７，并与不完全Ｃｈｅｌｅｓｋｙ共轭梯度法（ＩＣＣＧ）、Ｇａｕｓｓ直接解法进行比较．结果表明，ＭＧ方法计算速度明显优于ＩＣＣＧ、Ｇａｕｓｓ方法，对于１２９×１２９ｘ１２９网格的三维数值模拟费时４３ｓ，比ＩＣＣＧ法快７倍，而对于２５７×２５７×２５７超大型网格的三维数值模拟也仅需４１关键词
多重网格方法，泊松方程，三维数值模拟
文献标识码
Ａ
２
Ｓ．
中图分类号Ｐ６３１文章编号
１００４—２９０３（２００９）０１—０１５４—０５
Ｍｕｌｔｉｇｒｉｄｍｅｔｈｏｄ
（１．Ｓｃｈｏｏｌ
ｆｏｒ３Ｄ
ＷＵ
ｍｏｄｅｌｉｎｇｏｆＰｏｉｓｓｏｎｅｑｕａｔｉｏｎ
ＫｌａｕｓＳｐｉｔｚｅｒ２
ＬＵＪｉｎｇ—ｊｉｎｌ，
Ｘｉａｏ—ｐｉｎｇ¨，
ｏｆＥａｒｔｈａｎｄＳｐａｃｅＳｃｉｅｎｃｅ，ＵｎｉｖｅｒｓｉｔｙｏｆＳｃｉｅｎｃｅａｎｄＴｅｃｈｎｏｌｏｇｙｏｆＣｈｉｎａ，Ｈｅｆｅｉ２３００２６
２．Ｉｎｓｔｉｔｕｔ
ｆｉｉｒＧｅｏｐｈｙｓｉｋ，ＴＵＢｅｒｇｎ矗ｎｄＰ棚ｉＰＦｒｅｉｂｅｒｇ，Ｇｅｒｍａｎｙ）
Ａｂｓｔｒａｃｔｅｑｕａｔｉｏｎｉｓ
Ｉｎｔｈｉｓｐａｐｅｒ，ｗｅｂｒｉｅｆｌｙｉｎｔｒｏｄｕｃｅｔｈｅｂａｓｉｃｉｄｅａａｎｄｐｒｉｎｃｉｐｌｅｏｆｍｕｈｉｇｒｉｄ（ＭＧ）ｍｅｔｈｏｄ．Ｔｈｅｎ３Ｄｐｏｉｓｓｏｎ
ｔＯ
ｓｏｌｖｅｄｕｓｉｎｇＭｕｌｔｉｇｒｉｄ（ＭＧ）ｍｅｔｈｏｄ。ｗｉｔｈｇｒｉｄｓｃａｌｅｉｎｃｒｅａｓｉｎｇｆｒｏｍ１７×１７×１７
２５７×２５７×２５７．Ｉｎ
ｐａｒｉｓｏｎｗｉｔｈｉｎｐｌｅｔｅＣｈｏｌｅｓｋｙｃｏｎｊｕｇａｔｅｇｒａｄｉｅｎｔ（ＩＣＣＧ）ｍｅｔｈｏｄａｎｄＧａｕｓｓｄｉｒｅｃｔｍｅｔｈｏｄ，ｔｈｅｐｕｔｅｄ
ｓｐｅｅｄｏｆ
ＭＧｉｓｍｕｃｈｆａｓｔｅｒｔｈａｎＩＣＣＧａｎｄＧａｕｓｓｍｅｔｈｏｄ．ＭＧｍｅｔｈｏｄｔａｋｅｓ４３ｓｆｏｒ３Ｄｍｏｄｅｌｉｎｇｗｉｔｈ
ｔｉｍｅｓｆａｓｔｅｒｔｈａｎＩＣＣＧｍｅｔｈｏｄ，ａｎｄｏｎｌｙｔａｋｅｓ４１２ｓｆｏｒ３Ｄｐｒｏｂｌｅｍｗｉｔｈ
ａ
ａ
ｇｒｉｄｏｆ１２９
×１２９×１２９．７２５７．Ｋｅｙｗｏｒｄｓ
ｈｕｇｅｇｒｉｄｏｆ２５７×２５７×
ｍｕｈｉｇｒｉｄｍｅｔｈｏｄ，ｐｏｉｓｓｏｎｅｑｕａｔｉｏｎ，３Ｄｍｏｄｅｌｌｉｎｇ
０
引言
元方法，对方程及边界条件进行离散，最后都归结到解大型线性方程组Ｌｕ—ｆ．直接法（如Ｇａｕｓｓ法）求
地球物理重、磁、电勘探方法中涉及的地球物理场，如重力位、磁位、电位均满足泊松方程［１“ｏ．采用合适的数值方法，解决不同源项和边界条件下的泊
解大型方程组较困难，存在内存需求多、计算量大、
耗时长等问题．而一般的共轭梯度、ＳＯＲ迭代解法收敛速度较慢．吴小平等口￣６３引入了不完全
松问题，可分别实现对重力场、磁场、电场的正演模拟．随着三维勘探的展开，地球物理三维正反演越来
越重要，而首先要解决的就是三维正问题．因此，快
Ｃｈｏｌｅｓｋｙ共轭梯度（ＩＣＣＧ）法求解电法三维泊松问题，加快了收敛速度，明显提高了计算速度，应用于
电阻率三维反演也取得很好的效果．尽管如此，现阶段地球物理三维勘探面临越来越复杂的问题，其要求更精细的三维网格剖分．而ＩＣＣＧ方法随着网格节点的增加收敛速度也相应减慢，其对于应用中更
速、精确地求解三维泊松方程，在计算地球物理中具
有重要意义．
泊松方程三维数值模拟一般用有限差分、有限
收稿日期２００８基金项目
０５
１０；修回日期２００８０８２０．
国家自然基金（４０６７４０３７）和新世纪优秀人才支持计划项目（ＮＣＥＴ一０５—０５５３）资助．
作者简介鲁晶津，女，１９８３生，中德联合培养博士生，从事电磁三维数值模拟研究，现在ＴＵＢｅｒｇａｋａｄｅｍｉｅＦｒｅｉｂｅｒｇ，Ｇｅｒｍａｎｙ＊通讯作者
吴小平，男，１９６７年生，教授，博士生导师，主要从事电磁测深和地球内部物理研究．（Ｅ—ｍａｉｌ：ｗｘｐ＠ｕｓｔｃ．ｅｄｕ．）
万方数据
１期鲁晶津，等：三维泊松方程数值模拟的多重网格方法
复杂实际模型、更细密网格剖分以及更快速度要求，
难于有进一步提高．要有所突破，需借助计算数学最
新进展，引入新的高效算法．
多重网格（ＭＧ）方法在２０世纪６０年代由苏联
计算数学家Ｆｅｄｏｒｅｎｋｏ［７３提出最初思想，当时并未引起注意，其真正兴起是在近２０多年，Ｂｒａｎｄｔ［８］和Ｈａｃｋｂｕｓｃｈ［９１的工作才标志多重网格方法的研究全面开始．该方法的最大特点是，当离散化更加精细时，收敛速度并不变慢，即其收敛速度与网格尺度大小无关，且在求解更加复杂的问题时，这种有效性也不丧失．Ｈａｃｋｂｕｓｃｈ［９３等在数学上证明了至少对于
线形椭圆型偏微方程，多重网格算法是最优的，计算
量仅与网格节点数的一次方成正比，把解方程所需
计算量降到最低数量级，因而其它方法在计算效率
上很难与多重网格法相媲美．该方法是现今数值计
算最活跃的分支之一，在计算流体力学、结构力学、波动方程等诸多领域得到广泛应用，许多成果总结
在Ｈａｃｋｂｕｓｃｈ［９。、Ｂｒｉｇｇｓ［１０］、Ｗｅｓｓｅｌｉｎｇ［儿］、刘超
群［和Ｔｒｏｔｔｅｎｂｅｒｇ［１３３等专著中．在计算地球物理领域，该方法的应用也渐渐得到关注［１４￣１９Ｉ，不过，还
多限于二维问题．本文利用ＭＧ方法求解三维泊松方程，并较详细地比较了ＭＧ方法和ＩＣＣＧ、Ｇａｕｓｓ
法的内存需求、计算速度及精度，希望为多重网格
（ＭＧ）方法应用于地球物理三维快速正演计算取到抛砖引玉的作用．
１
多重网格（ＭＧ）方法［－９～１３］
多重网格法最初动力来源于对网格方程迭代求
解时，误差的各个Ｆｏｕｒｉｅｒ分量的不同衰减程度．分析线性方程组收敛特征时，把全部误差分量分为两类：光滑分量和高频分量．把频率低、变化慢、看上去
较光滑的分量称为光滑分量；而把频率高、变化快的
分量称为高频分量．对于某一固定网格，误差分量由光滑分量和高频分量组成．高频分量在传统迭代（如快，而光滑分量在迭代过程中衰减很慢，这就是为什么迭代法求解线性方程组在开始迭代时收敛速度相当快，随着迭代的进行，收敛速度越来越慢的原因．多重网格法为了克服固定网格的缺点，先在较细网
格上进行迭代，把高频分量衰减掉，然后在较粗网格
上迭代，把次高频分量衰减掉，逐层变粗直到最粗一层网格，把各种频率分量衰减掉．再由粗网格开始，依次返回到各级细网格，最后在最细网格上获得所万方数据
格Ｖ循环（图ｌａ）过程简述如下：
令最细网格层的线性方程表示为Ｌ＾“＾一ｆｈ，
（１）
其中Ｌ。为离散微分算子，“。为方程精确解，ｈ为细
网格单元长度．在细网格上进行ｅ（一般ｅ一１或２）次松弛光滑迭代（本文采用红黑Ｇａｕｓｓ—Ｓｉｅｄｅｌ迭代法［９?２０３），使迭代近似解霞。与“。间误差的高频分量衰减而变得光滑．其残差方程为
Ｌ＾口＾一ｄ＾，
（２）
其中仇一“。一霞。为误差，可视为对面。的校正量，
ｄ。一ｆｈ—Ｌ。葫。为残差向量．正常情况下迭代法解（２）式得到饥的近似解乞，。，然后磊：…一霞。＋饥形成迭代．而多重网格法完全是另一种新的思路，将（２）式的残差方程近似到更粗网格上解，即
ＬＨ训Ｈ—ｄＨ，
（３）
通常粗网格单元长度Ｈ一２ｈ，则对于三维网格，单元数减少至上一级细网格的１／８，在粗网格上
解（３）式较解（２）式容易得多．从ｄ。到ｄＨ通过所谓
限制（Ｒｅｓｔｒｉｃｔｉｏｎ）算子Ｒ，
ｄＨ—Ｒｄ＾，
（４）
解（３）式得到一个解ｏＨ，再通过插值（Ｐｒｏｌｏｎｇａｔｉｏｎ）算子Ｐ，得到％的近似解ｏ。，
讥一心Ｈ，
（５）
限制算子Ｒ和插值算子Ｐ都是线性算子［９］，（４）和
（５）式容易计算．修正近似解葫。为
“－－＾ｎ。”一面＾＋矛＾，
（６）
称为粗网格校正．若残差方程（３）在次一级粗网格上
就精确求解，则形成二层网格法，见图１中的二层网格（２一ｇｒｉｄ），Ｓ表示松弛光滑迭代，斜杠＼表示限制算子，Ｅ表示精确求解，斜杠／表示插值算子．对于三层网格，就是在次一级粗网格仍不精确求解，而是到更次一层进行粗网格校正，即二层网格中的精确解Ｅ由下一个二层网格代替，形成三层网格法．多层网格法以此类推（见图ｌａ），形成多重网格算法的Ｖ循环，这一递归性质则称为套迭代技术．到最粗一层网格时，节点已经很少，直接精确求解的计算量几乎可以忽略．
图１ａ的Ｖ循环是从随机给定的初始值开始松
弛迭代的．图１ｂ是另一种称为完全多重网格的循环
模式，先在最粗网格上直接求解，然后将其逐级延拓
到最细一层的网格上，得到更接近真值的初始解，有利于之后Ｖ循环的松弛迭代．
多重网格法正是利用细网格上的松弛光滑特性消除迭代误差的高频分量，利用粗网格上的残差校
Ｊａｃｏｂｉ，Ｇａｕｓｓ—Ｓｅｉｄｅｌ，ＳＯＲ迭代等）过程中衰减很需方程的解，详见Ｈａｃｋｂｕｓｃｈ［９］．一个简单的多重网
地球物理学进展
正特性消除迭代误差的低频分量，套迭代技术通过
限制和插值算子连接所有网格层共同求解，如此消除了不同频率范围的误差分量，加速了迭代的收敛
速度．
２
图１ａ多重网格算法的Ｖ循环
Ｆｉｇ．１
ａ
ＲｅｃｕｒｒｅｎｃｅＶ
ｏｆｍｕｈｉｇｒｉｄａｌｇｏｒｉｔｈｍ
图１ｂ完全多重网格算法的循环
ｓ为松弛迭代，斜杠＼为限制算子，斜杠／为延拓算子，Ｅ为精确求解
Ｆｉｇ．１ｂ
Ｒｅｃｕｒｒｅｎｃｅｏｆ
ｆｕｌｌｍｕｌｔｉｇｒｉｄａｌｇｏｒｉｔｈｍ
Ｓ－ｒｅｌａｘａｔｉｏｎ
ｉｔｅｒａｔｉｏｎ，＼ＩＲｅｓｔｒｉｃｔｉｏｎ，／－Ｐｒｏｌｏｎｇａｔｉｏｎ，Ｅ—Ｅｘａｃｔｓｏｌｖｅｒ
２
三维泊松方程的多重网格算法
考虑在区域Ｑ及其边界ａＱ上的如下三维泊松
方程边值问题：
ｆＬｕ一３７ｒ２ｓｉｎＯｒｘ）ｓｉｎＯｒｙ）ｓｉｎ（＿ｒｚ）
１甜Ｉ加一０
’
（ｚ，Ｙ，ｚ）∈Ｑ—Ｌ０，１Ｊ
３
该问题有理论解ｕ（ｘ，Ｙ，ｚ）一ｓｉｎ（ｚｒｘ）ｓｉｎ（ｚｒｙ）ｓｉｎ（ｚｒｚ）．
在区域ｎ上对泊松方程采用中心有限差分法进行均匀离散，分别用ＭＧ法、ＩＣＣＧ法以及Ｇａｕｓｓ法求解离散所得的大型线性方程组，其中ＩＣＣＧ及
ＭＧ算法先在最粗网格３×３×３（即最粗网格
步长为０．５）上直接求解，实际只是一个方程解一个１ｂ的完全多重网格的循环模式求解，当迭代误差ａ小于阈值ｄ。一２×１０＿７时，计算终止，ａ表示为
万方数据
图２随三维网格的加密，ＭＧ的收敛曲线
Ｆｉｇ．２
Ｃｏｎｖｅｒｇｅｎｃｅ
ｃｕｒｖｅ
ｏｆｍｕｌｔｉｇｒｉｄ
ａｌｇｏｒｉｔｈｍｆｏｒｖａｒｉｏｕｓｇｒｉｄｓｉｚｅ
图３
随三维网格的加密，ＩＣＣＧ的收敛曲线
Ｆｉｇ．３
Ｃｏｎｖｅｒｇｅｎｃｅ
ｃｕｒｖｅ
ｏｆ１ＣＣＧｍｅｔｈｏｄ
ｆｏｒｖａｒｉｏｕｓｇｒｉｄｓｉｚｅｓ
图４ＭＧ、ＩＣＣＧ、Ｇａｕｓｓ法的计算时间随网格
节点数的变化曲线
Ｆｉｇ．４
Ｒｕｎｎｉｎｇ
ｔｉｍｅ
ｔ
Ｖｓ．ｔｈｅ
ｎｕｍｂｅｒｏｆｇｒｉｄｎｏｄｅｓＮ
ｆｏｒ
ＭＧ、ＩＣＣＧａｎｄＧａｕｓｓ
ｍｅｔｈｏｄｓ
Ｇａｕｓｓ方法本文不讨论，详见文献［４，５］．
未知数，将其逐级延拓到更细一层网格上，其步长ｈ是上一层粗网格步长Ｈ的一半（Ｈ一２ｈ），然后按图
１期鲁晶津，等：三维泊松方程数值模拟的多重网格方法
ａ一．／≥：ｒ２（ｉ，ｊ，走）／≥：ｒｊ（ｉ，Ｊ，忌），
。Ｖ万■
ｚ：ｊ：；
其中ｒ。（ｉ，Ｊ，是）、，．（ｉ，Ｊ，ｋ）分别代表一次迭代前后网格节点（ｉ，Ｊ，ｋ）处的残差．
图２和图３分别是ＭＧ和ＩＣＣＧ迭代的收敛曲线，从１７×１７×１７到２５７×２５７×２５７的不同网格尺度，ＭＧ算法均是７次迭代收敛，证明了ＭＧ的收敛速度与网格尺度无关．而ＩＣＣＧ随着网格逐次加密，迭代次数明显增多，收敛曲线明显变慢．
表１是在不同网格尺度上ＭＧ法、ＩＣＣＧ法以及Ｇａｕｓｓ法的ＰＣ（ＣＰＵ２．８０ＧＨＺ，内存１ＧＢ）运行时间，在相同尺度的网格上，ｔＭ。最短，ｔ。。。Ｇ为ｔＭ。的数倍，ｔＧ。。。，最大可达１０００ｔＭＧ．图４显示了各算法运行时间随节点数Ｎ的变化，ｔＭＧ随Ｎ线性增加，而ｔｌ（ｃ。随Ｎ明显是非线性增长，在３３×３３×３３网格上ｔｌｃｃＧ／ｔＭＧ一２，６５×６５×６５网格上ｔＩｃｃＧ／ｔＭＧ一３，到１２９×１２９×１２９网格上时ｔＴｃ∞／ｔＭ。一７，随着网格加密ＭＧ法的优势越明显．这主要是由于ＭＧ的收敛速度与网格尺度无关，而ＩＣＣＧ的迭代次数随着网格加密明显增多造成的．对于２５７×２５７×２５７网格，
由于ＭＧ方法不需存储系数矩阵，只需存储各层网
格节点上的场值及右端项，所需内存最少，因而多重网格法在内存１ＧＢ的ＰＣ上仍能运行，费时仅４１２
Ｓ．
ＩＣＣＧ和Ｇａｕｓｓ法的内存需求详见文献［４］，尤其Ｇａｕｓｓ直接解法所需内存巨大，对于６５×６５×６５网格机器内存即不能满足要求，计算时间更是呈级数增长，难于解决大型的三维数值模拟．
表１不同三维网格尺度上ＭＧ、ＩＣＣＧ
和Ｇａｕｓｓ算法的计算时间
Ｔａｂｌｅ１
Ｒｕｎｎｉｎｇ
ｔｉｍｅｏｆＭＧ、ＩＣＣＧ、Ｇａｕｓｓ
ｍｅｔｈｏｄｓｆｏｒｖａｒｉｏｕｓ３Ｄｇｒｉｄｓ
表２是各方法数值解与理论解之间的最大相对误差，表示为
ｅｔ＇ｒ—ｍａｘ（ＪＵ。ｈ一“。。。Ｉ／ｕ。ｈ），
ｔ，，，ｋ
其中‰代表理论解，Ｒｎｕｍ代表数值解．结果表明ＭＧ
法达到Ｇａｕｓｓ直接求解的精确，并随网格加密精度
越高，对于３３×３３×３３网格计算得到的数值解误差
万方数据
不大于０．０８０４％．
表２
ＭＧ、ＩＣＣＧ、Ｇａｕｓｓ数值解的误差
Ｔａｂｌｅ２
Ｅｒｒｏｒｓｉｎｎｕｍｅｒｉｃａｌｓｏｌｕｔｉｏｎｓｏｆ
ＭＧ、ＩＣＣＧ、Ｇａｕｓｓｍｅｔｈｏｄｓ
需要指出的是，以上三维泊松问题是在正方形区域的均匀网格下实施的，也没有物性不均匀性问
题，多重网格算法是最快捷的，ＩＣＣＧ方法也能达到其最佳效率．然而在地球物理（如地电）三维模拟中，
地下存在复杂的电性不均匀结构，电性差异甚至可
达几个数量级；同时地球物理模拟通常是无界问题，
因此非均匀网格剖分也不可避免，这些将导致多重
网格法和ＩＣＣＧ算法的效率下降．不同的是，预条件共扼梯度法对此没有更进一步解决方案，而多重网格方法通过松驰过程、限制算子、插值算子等相应改
进凹“１’”一“２３］，使得在各种复杂情况之下，其有效性并不丧失，多重网格方法仍可得到实现，并重新获得其固有的高效性［９’１１‘１“．当然，这将导致整个多重网格算法技术深度和复杂程度的上升，而这些深层研究，往往属于计算数学界的专家们［９’１１’”＇２ｚ－ｚｓ３，
其它领域似乎难于问津．正因为如此，多重网格方法
迄今还极少应用于（地球）电磁三维数值模拟．
３
结论
（１）由于不需存储系数矩阵，多重网格法所耗内存最少，因而能够用细密的网格剖分模拟更复杂模型．
（２）在同样的收敛准则下，ＭＧ法收敛速度与网格尺度无关，而ＩＣＣＧ法收敛速度则随网格加密而迅速减小．
（３）同样网格尺度下，ＭＧ法计算速度比ＩＣＣＧ
法快数倍，并且随着网格加密优势越明显．Ｇａｕｓｓ法
的计算量明显太大，不宜于大规模数值计算．
（４）误差分析也表明，多重网格法是精确、高效的三维数值模拟算法．
可见，多重网格方法有其固有的优越性，但针对
实际地球物理正演计算，许多问题还有待进一步深
１５８
地球物理学进展
２４卷
入研究．
参考文献（Ｒｅｆｅｒｅｎｃｅｓ）：
Ｅ１］
周熙襄，钟本善．电法勘探数值模拟技术［Ｍ］．四川科学技术
出版社，１９８６．
Ｚｈｏｕ
ＸＸ，ＺｈｏｎｇＢＳ．ＮｕｍｅｒｉｃａｌＳｉｍｕｌａｔｉｏｎ
ｔｅｃｈｎｉｑｕｅｓｉｎ
ｅｌｅｃｔｒｉｃａｌ
ｐｒｏｓｐｅｃｔｉｎｇ［Ｍ］．Ｃｈｅｎｇｄｕ：Ｓｉｃｈｕａｎｐｕｂｌｉｓｈｉｎｇｈｏｕｓｅ
ｏｆＳｃｉｅｎｃｅａｎｄＴｅｃｈｎｏｌｏｇｙ，１９８６．
Ｌ２］徐世浙．地球物理中的有限单元法［Ｍ］．北京；科学出版社，
１９９４．
ＸｕＳＺ．Ｆｅｍｉｎ
Ｇｅｏｐｈｙｓｉｃｓ［Ｍ］．Ｂｅｉｊｉｎｇ：Ｓｃｉｅｎｃｅ
Ｐｒｅｓｓ，１９９４．
［３］阮百尧，熊彬，徐世浙．三维地电断面电阻率测深有限元数值模拟［Ｊ］．地球科学，２００１，２６（１）：７３～７７．
Ｒｕａｎ
ＢＹ，Ｘｉｏｎｇ
Ｂ，ＸｕＳＺ．Ｆｉｎｉｔｅｅｌｅｍｅｎｔｍｅｔｈｏｄ
ｆｏｒ
ｍｏｄｅｌｉｎｇｒｅｓｉｓｔｉｖｉｔｙｓｏｕｎｄｉｎｇｏｎ
３－Ｄｇｅｏｅｌｅｃｔｒｉｃ
ｓｅｃｔｉｏｎ［Ｊ］．
ＥａｒｔｈＳｃｉｅｎｃｅｓ，２００１，２６（１）：７３～７７．
［４］吴小平，徐果明，李时灿．利用不完全Ｃｈｏｌｅｓｋｙ共轭梯度法求解
点源三维地电场口］．地球物理学报，１９９８，４１（６）：８４８～８５５．
ＷｕＸＰ，ＸｕＧＭ，ＬｉＳＣ．Ｔｈｅｃａｌｃｕｌａｔｉｏｎｏｆ３－Ｄｇｅｏｅｌｅｃｔｒｉｃ
ｆｉｅｌｄ
ｙｉｅｌｄｅｄ
ｂｙ
ｐｏｉｎｔ
ｓｏｕｒｃｅ
ｕｓｉｎｇｉｎｐｌｅｔｅｃｈｏｌｅｓｋｙ
ｃｏｎｊｕｇａｔｅｇｒａｄｉｅｎｔｍｅｔｈｏｄ［Ｊ］．ＣｈｉｎｅｓｅＪ．Ｇｅｏｐｈｙｓ．，１９９８，４１
（６）：８４８～８５５．
［５］
ＷｕＸＰ，ＸｉａｏＹ，ＱｉＣ，ＷａｎｇＴ．Ｃｏｍｐｕｔａｔｉｏｎｓｏｆｓｅｃｏｎｄａｒｙ
ｐｏｔｅｎｔｉａｌｆｏｒ３－ＤＤＣｒｅｓｉｓｔｉｖｉｔｙｍｏｄｅｌｌｉｎｇｕｓｉｎｇａｎ
ｉｎｐｌｅｔｅ
Ｃｈｏｌｅｓｋｙ
ｃｏｎｊｕｇａｔｅｇｒａｄｉｅｎｔｍｅｔｈｏｄ［Ｊ］．Ｇｅｏｐｈｙｓ．Ｐｒｏｓｐｅｃｔ．，
２００３。５１（６）：５６７～５７７．
［６］
ＷｕＸＰ．Ａ３－Ｄｆｉｎｉｔｅ—ｅｌｅｍｅｎｔａｌｇｏｒｉｔｈｍｆｏｒＤＣｒｅｓｉｓｔｉｖｉｔｙ
ｍｏｄｅｌｉｎｇｕｓｉｎｇｓｈｉｆｔｅｄｉｎｐｌｅｔｅｅｈｏｌｅｓｋｙ
ｃｏｎｊｕｇａｔｅｇｒａｄｉｅｎｔ
ｍｅｔｈｏｄ［Ｊ］．Ｇｅｏｐｈｙｓ．Ｊ．Ｉｎｔ，２００３，１５４（３）：９４７～９５６．
［７］
ＦｅｄｏｒｅｎｋｏＲＰ．Ｔｈｅ
ｓｐｅｅｄ
ｏｆ
ｃｏｎｖｅｒｇｅｎｃｅ
ｏｆ
ｏｎｅ
ｉｔｅｒａｔｉｖｅ
ｐｒｏｃｅｓｓ［Ｊ］．Ｚｈ．Ｖｙｃｈ．Ｍａｔ，１９６４，４（５）：５５９—５６４．
［８１
Ｂｒａｎｄｔ
Ａ．Ｍｕｌｔｉ—ｌｅｖｅｌａｄａｐｔｉｖｅｓｏｌｕｔｉｏｎｓ
ｔｏ
ｂｏｕｎｄａｒｙ－ｖａｌｕｅ
ｐｒｏｂｌｅｍｓ［Ｊ］．Ｍａｔｈ．Ｃｏｍｐ，１９７７，３１（１３８）：３３３—３９０．
［９］
Ｈａｃｋｂｕｓｅｈ
Ｗ．Ｍｕｈｉｇｒｉｄ
ｍｅｔｈｏｄｓａｎｄ
ａｐｐｌｉｃａｔｉｏｎ［Ｍ］．
ＢｅｒｌｉｎＩＳｐｒｉｎｇｅｒ—Ｖｅｒｌａｇ，１９８５．
［ｉｏ］ＢｒｉｇｇｓＷＬ．Ａｍｕｌｔｉｇｒｉｄｔｕｔｏｒｉａｌ［Ｍ］．ＳＩＡＭ，Ｐｈｉｌａｄｅｌｐｈｉａ，
１９８７．
［１１］
ＷｅｓｓｅｌｉｎｇＰＷ．Ａｎｉｎｔｒｏｄｕｃｔｉｏｎｔｏ
ｍｕｌｔｉｇｒｉｄ
ｍｅｔｈｏｄｓ［Ｍ］．
Ｃｈｉｃｈｅｓｔｅｒ：Ｗｉｌｅｙ＆Ｓｏｎｓ，１９９２．
［１２］刘超群．多重网格法及其在计算流体力学中的应用ｌ－Ｍ］．北京：清华大学出版社，ｉ９９５．
ＬｉｕＣＱ．Ｍｕｌｔｉｇｒｉｄｍｅｔｈｏｄａｎｄａｐｐｌｉｃａｔｉｏｎｉｎｐｕｔａｔｉｏｎａｌｆｌｕｉｄ
ｄｙｎａｍｉｃｓ［Ｍ］．Ｂｅｉｊｉｎｇ；ＴｓｉｎｇｈｕａＵｎｉｖｅｒｓｉｔｙＰｒｅｓｓ，
１９９５．（ｉｎＣｈｉｎｅｓｅ）
［１３］ＴｒｏｔｔｅｎｂｅｒｇＵ，ＯｏｓｔｅｒｌｅｅＣＷ，ＳｃｈｕｌｌｅｒＡ．Ｍｕｌｔｉｇｒｉｄ［Ｍ］．
ＡｃａｄｅｒｎｉｃＰｒｅｓｓ，Ｉｎｃ，２００１．
［１４］ＢｕｎｋｓＣ．ＳａｌｅｅｋＦＭ，ＺａｌｅｓｋｉＳ，ＣｈａｖｅｎｔＧ．Ｍｕｌｔｉｓｃａｌｅ
ｓｅｉｓｍｉｃｗａｖｅｆｏｒｍｉｎｖｅｒｓｉｏｎ［Ｊ］．Ｇｅｏｐｈｙｓｉｃｓ．，１９９５，６０（５）：
１４５７～１４７３．
［１５］刘伊克，常旭．盆地模拟水动力油气二次运移隐式多重网格法［Ｊ］．地球物理学报，１９９８，４１（３）：３４２～３４８．
Ｌｉｕ万方数据
ＹＫ，ＣｈａｎｇＸ．Ｔｈｅｉｍｐｌｉｃｉｔ
ｍｕｈｉｇｒｉｄｍｅｔｈｏｄｆｏｒｔｈｅ
ｓｅｃｏｎｄａｒｙｈｙｄｒｏｃａｒｂｏｎｍｉｇｒａｔｉｏｎ
ｄｒｉｖｅｄ
ｂｙ
ｗａｔｅｒ［Ｊ］．
ＣｈｉｎｅｓｅＪ．Ｇｅｏｐｈｙｓ．，１９９８，４１（３）：３４２～３４８．
［１６］ＰｌａｋｓＡ，ＴｓｕｋｅｒｍａｎＩ，ＰａｉｎｃｈａｕｄＳ。ＴａｂａｒｏｖｓｋｙＬ．Ｍｕｌｔｉｇｒｉｄ
ｍｅｔｈｏｄｓｆｏｒｏｐｅｎ
ｂｏｕｎｄａｒｙｐｒｏｂｌｅｍｓ
ｉｎ
ｇｅｏｐｈｙｓｉｃｓ［Ｊ］．ＩＥＥＥ
Ｔｒａｎｓａｃｔｉｏｎｓ
ｏｎ
Ｍａｇｎｅｔｉｃｓ，２０００．３６（４）：６３３～６３８．
［１７］
ＭｏｕｃｈａＲ，ＢａｉｌｅｙＲＣ．Ａｎ
ａｃｃｕｒａｔｅ
ａｎｄ
ｒｏｂｕｓｔｒｎｕｈｉｇｒｉｄ
ａｌｇｏｒｉｔｈｍｆｏｒ２Ｄｆｏｒｗａｒｄｒｅｓｉｓｔｉｖｉｔｙｍｏｄｅｌｉｎｇ［Ｊ］
Ｇｅｏｐｈｙｓ．
Ｐｒｏｓｐｅｃｔ，２００４．５２（３）：１９７～２１２．
［１８］
田小波，吴庆举，曾融生．弹性波场数值模拟的隐式差分多重网格算法ＦＪ］．地球物理学报，２００４，４７（１）：８１～８７．
ＴｉａｎＸ
Ｂ，ＷｕＱＪ，ＺｅｎｇＲＳ．Ｍｕｌｔｉ—ｇｒｉｄａｌｇｏｒｉｔｈｍｆｏｒ
ｎｕｍｅｒｉｃａｌｍｏｄｅｌｉｎｇｏｆｅｌａｓｔｉｃｗａｖｅｆｉｅｌｄｕｓｉｎｇｆｉｎｉｔｅ—ｄｉｆｆｅｒｅｎｃｅ
ｍｅｔｈｏｄ［Ｊ］．ＣｈｉｎｅｓｅＪ．Ｇｅｏｐｈｙｓ．，４７（１）：８１～８７，２００４．［１９］
ＭｕｌｄｅｒＷ
Ａ．Ａ
ｍｕｈｉｇｒｉｄ
ｓｏｌｖｅｒ
ｆｏｒ
３Ｄ
ｅｌｅｃｔｒｏｍａｇｎｅｔｉｃ
ｄｉｆｆｕｓｉｏｎ［Ｊ］．Ｇｅｏｐｈｙｓ．Ｐｒｏｓｐｅｃｔ．，２００６，５４（５）：６３３～６４９．
［２０］
ＺｈａｎｇＪ．Ｆａｓｔａｎｄ
ｈｉｇｈ
ａｃｃｕｒａｃｙ
ｍｕｈｉｇｒｉｄ
ｓｏｌｕｔｉｏｎｏｆｔｈｅ
ｔｈｒｅｅｄｉｍｅｎｓｉｏｎａｌｐｏｉｓｓｏｎ
ｅｑｕａｔｉｏｎ［Ｊ］．Ｊ．Ｃｏｍｐ．Ｐｈｙｓ．，
１９９８，１４３：４４９～４６１．
［２１］
张智诠，周立伟，黄应清，等．轴对称静电场不等距多重网格
法的研究［Ｊ］．装甲兵工程学院学报，２０００，１４（２）：３５～３９．
ＺｈａｎｇＺＱ，ＺｈｏｕＬＷ，ＨｕａｎｇＹＱ．Ａｍｕｈｉｇｒｉｄｍｅｔｈｏｄ
ｗｉｔｈ
ｎｏｎ—－ｕｎｉｆｏｒｍ
ｍｅｓｈｅｓ
ｆｏｒ
ｔｈｅ
ａｘｉａｌ—－ｓｙｍｍｅｔｒｉｃａｌ
ｅｌｅｃｔｒｏｓｔａｔｉｃ
ｆｉｅｌｄｓ［Ｊ］．ＪｏｕｒｎａｌｏｆＡｃａｄｅｍｙｏｆＡｒｍｏｒｄｅＦｏｒｃｅ
Ｅｎｇｉｎｅｅｒｉｎｇ，２０００，１４（２）：３５～３９．
［２２３
黄朝晖，常谦顺．多重网格插值算子的改进算法［Ｊ］．应用数
学学报，２００３，２６（３）．
ＨｕａｎｇＣＨ，ＣｈａｎｇＱＳ．ＡｎｉｍｐｒｏｖｅｄａｌｇｏｒｉｔｈｍｆｏｒｔｈｅＭＧ
ｉｎｔｅｒｐｏｌａｔｉｏｎ
ｏｐｅｒａｔｏｒ［Ｊ］．ＡｃｔａＭａｔｈｅｒｎａｔｉｃａｅ
Ａｐｐｌｉｃａｔａｅ
Ｓｉｎｉｃａ。２００３，２６（３）．
［２３］
刘之行，封卫兵．三维区域上的多重网格算法［Ｊ］．西安交通大学学报，１９９８，３２（１２）：９４～９７．
ＬｉｕＺＸ．Ｆｅｎｇ
ＷＢ．Ｍｕｌｔｉｇｒｉｄｍｅｔｈｏｄａｐｐｌｉｅｄ
ｔｏｔｈｒｅｅ—
ｄｉｍｅｎｓｉｏｎａｌ
ｄｏｍａｉｎｓ［Ｊ］．
Ｊｏｕｒｎａｌ
ｏｆ
Ｘｉ’ａｎ
Ｊｉａｏｔｏｎｇ
Ｕｎｉｖｅｒｓｉｔｙ．１９９８。３２（１２）：９４～９７．
［２４］
ＺｈａｎｇＪ．Ａｃｃｅｌｅｒａｔｅｄｍｕｈｉｇｒｉｄｈｉｇｈ
ａｃｃｕｒａｃｙ
ｓｏｌｕｔｉｏｎｏｆｔｈｅ
ｃｏｎｖｅｃｔｉｏｎ—ｄｉｆｆｕｓｉｏｎｅｑｕａｔｉｏｎｗｉｔｈｈｉｇｈＲｅｙｎｏｌｄｓｎｕｍｂｅｒ［Ｊ］．
Ｎｕｍｅｒｉｃａｌｍｅｔｈｏｄｓｆｏｒｐａｒｔｉａｌｄｉｆｆｅｒｅｎｔｉａｌｅｑｕａｔｉｏｎｓ，１９９７，
１３：７７～９２．
［２５］
Ａｉ
ＭＨ，ＥｓｓｅｌａｏｕｉＤ，ＨａｋｉｍＡ，ＲａｇｈａｙＳ．Ｆｉｎｉｔｅｖｏｌｕｍｅ
ｍｕｈｉｇｒｉｄ
ｍｅｔｈｏｄ
ｏｆ
ｔｈｅ
ｐｌａｎａｒ
ｃｏｎｔｒａｃｔｉｏｎ
ｆｌｏｗｏｆ
ａ
ｖｉｓｃｏｅｌａｓｔｉｃ
ｆｌｕｉｄ［Ｊ］．Ｉｎｔｅｒｎａｔｉｏｎａｌ
Ｊｏｕｒｎａｌ
ｆｏｒＮｕｍｅｒｉｃａｌ
ＭｅｔｈｏｄｓｉｎＦｌｕｉｄｓ，２００１，３６（８）：８８５～９０２．
［２６］
Ａｒｍｅｎｉｏ
Ｖ．Ａｎｉｍｐｒｏｖｅｄ
ＭＡＣ
ｍｅｔｈｏｄ（ＳＩＭＡＣ）ｆｏｒｕｎｓｔｅａｄｙｈｉｇｈ—ｒｅｙｎｏｌｄｓ
ｆｒｅｅ
ｓｕｒｆａｃｅ
ｆｌｏｗｓ［Ｊ］．Ｉｎｔｅｒｎａｔｉｏｎａｌ
ＪｏｕｒｎａｌｆｏｒＮｕｍｅｒｉｃａｌＭｅｔｈｏｄｓｉｎＦｌｕｉｄｓ，１９９７，２４（２）：１８５～２１４．
［２７］
ＡｒｕｌｉａｈＤ，Ａｓｃｈｅｒ
Ｕ．Ｍｕｌｔｉｇｒｉｄｐｒｅｃｏｎｄｉｔｉｏｎｉｎｇｆｏｒｔｉｍｅ
ｈａｒｍｏｎｉｃ
ｍａｘｗｅｌｌ’ｓ
ｅｑｕａｔｉｏｎｓｉｎ３Ｄ［Ｊ］．ＳＩＡＭＪｏｕｒｎａｌｏｆ
ＳｃｉｅｎｔｉｆｉｃＣｏｍｐｕｔｉｎｇ，２００２，２４：７０２～７１８．
［２８］
ＫａｎｓｃｈａｔＧ．Ｍｕｌｔｉｌｅｖｅｌｍｅｔｈｏｄｓｆｏｒｄｉｓｃｏ力ｔｉｎｕｏｕｓＧａｌｅｒｋｉｎ
ＦＥＭ
ｏｎ
ｌｏｃａｌｌｙｒｅｆｉｎｅｄ
ｍｅｓｈｅｓ［Ｊ］．Ｃｏｍｐｕｔｅｒｓ＆Ｓｔｒｕｃｔｕｒｅｓ，
２００４。８２（２８）：２４３７～２４４５．
范文五：线性方程组的数值解法
《大学数学实验》作业
线性方程组的数值解法
班级： 姓名： 学号： 日期：
目录
目录 .................................................................................................................................................. 2 【实验目的】 ................................................................................................................................... 3 【实验内容】 ................................................................................................................................... 3
【题目1】（课本习题第五章第3题） .................................................................................. 3
【第(1)问求解】 .............................................................................................................. 3 【第(2)问求解】 ............................................................................................................ 11 【本题小结】 ................................................................................................................. 12 【题目2】（课本习题第五章第5题） ................................................................................ 13
【第(1)问求解】 ............................................................................................................ 14 【第(2)问求解】 ............................................................................................................ 15 【本题小结】 ................................................................................................................. 16 【题目3】（课本习题第五章第7题） ................................................................................ 16
【第（1）问求解】 ....................................................................................................... 16 【第（2）问求解】 ....................................................................................................... 18 【本题小结】 ................................................................................................................. 19
【实验心得、体会】 ..................................................................................................................... 19
注：
本实验作业matlab脚本文件均以ex5_3形式命名，其中ex代表作业，5_3表示第五章第三小题
【实验目的】
1.学会用MATLAB软件数值求解线性代数方程组，对迭代法的收敛性和解的稳定性作初步分析；
2.通过实例学习用线性代数方程组解决简化的实际问题。
【实验内容】
【题目1】（课本习题第五章第3题）
已知方程组Ax=b,其中A?R
20?20
，定义为
?1/2?1/4?3
??1/23?1/2?1/4?
??1/4?1/23?1/2?1/4?
1/4?1/23??
A???1/4??
?
???
??????
?????
????????
??
3?1/2?1/4?
?
?1/23?1/2??1/4?1/23??
试通过迭代法求解此方程组，认识迭代法收敛的含义以及迭代初值和方程组系数
矩阵性质对收敛速度的影响。实验要求：
(1)选取不同初始向量x和不同的方程组右端向量b,给定迭代误差要求，用雅（0）
克比和高斯——赛德尔迭代法计算，观测得到的迭代向量序列是否均收敛？若收敛，记录迭代次数，分析计算结果并得出你的结论；
(2)取定右端向量b和初始向量x(0)，将A的主对角元素成倍增长若干次，非主
(k?1)(k)?5|x?x||?10对角元素不变，每次用雅克比迭代法，要求迭代误差满足,
比较收敛速度，分析现象并得出你的结论。
【第(1)问求解】
【模型分析】
对于矩阵A，将其分解为：
A=D-L-U,
其中：
?3?0?D=??
??0??0
0?
0?0??0
?30?0?0?0?
0?00?, L=?
?????030?
?0??
0?03???0???
?0
0????0?0?
??，U=?
00??
??
?00?
?0?0?0???
0?00
0?
?
?0?????
??
?0?
?00?? 0
其中D、L、U均为20×20的方阵。
则雅可比迭代公式为
x(k?1)?D?1(L?U)x(k)?D?1b, k?0,1,2...
高斯—赛德尔迭代公式为：
x(k?1)?D?1(Lx(k?1)?Ux(k))?D?1b, k?0,1,2...
选取初始向量x
(0)
=[0,0,0?0]?，b=[1,1,1?1]?，迭代误差为：x(k?1)?x(k)
?
?10?5，
用雅可比迭代公法和高斯—赛德尔迭代法对方程组Ax=b进行迭代计算，在Matlab
中编写如下的程序代码： 【matlab程序】
%------------------------------作业题5_3脚本M文件源程序---------------------------- clear all;clc; n=20;
A1=sparse(1:n,1:n,3,n,n); % 输入A的对角元素 A2=sparse(1:n-1,2:n,-0.5,n,n); % 输入A的(上)次对角元素 A3=sparse(1:n-2,3:n,-0.25,n,n); % 输入(上)次次对角元素 A=A1+A2+A3+A2'+A3'; % 得到A矩阵 D=diag(diag(A)); % 提取D矩阵 L=-tril(A,-1); % 提取L矩阵 U=-triu(A,1); % 提取U矩阵
b= ones(20,1); % 设定b初值（进行变化） x(:,1)=zeros(n,1); % 设定初始向量（进行变化） for i=2:1000
x(:,i)=D\ (L+U)*x(:,i-1)+D\b; % 按雅克比迭代公式进行计算 if norm(x(:,i)-x(:,i-1),inf)
Bg_s= (D-L)\U; Fg_s= (D-L)\b;
y(:,1)=x(:,1); % 设定初始向量
for j=2:1000
y(:,j)=Bg_s*y(:,j-1)+Fg_s; % 按高斯—赛德尔迭代公式进行计算 if norm(y(:,j)-y(:,j-1),inf)
[i-1,j-1] % 输出两次迭代的次数
[x(:,i),y(:,j), A\b] % 输出迭代得到的方程的解及精确解 x,y % 输出两次迭代的所有结果
plot(1:i,x(1,:),'r',1:j,y(1,:),'b'); % 作出x1在两种迭代过程中的变化图像 xlabel('迭代次数');ylabel('x1');title('x1在两种迭代过程中的变化图像'); % 加入X轴、Y轴标记和标题
legend('雅克比迭代法','高斯—赛德尔迭代法', 'location', 'south') % 加入图例 gtext('x1');
可以得到，所得的向量序列均是收敛的。具体的原因见后面的理论分析部分。 用雅可比迭代公式计算，共迭代16次，得到的迭代结果如下表一：
表1：取x
(0)
=[0,0,0?0]?，b=[1,1,1?1]?雅可比迭代公式计算部分结果
用高斯—赛德尔迭代法计算，共迭代11次，得到的部分迭代结果如下表2：
表2：取x
(0)
=[0,0,0?0]?，b=[1,1,1?1]?高斯—赛德尔迭代法计算部分结果
?
x(0)=[0,0,0?0]?，b=[1,1,1?1]?，x(k?1)?x(k)
?10?5时，为了更清晰比较两种迭代
方法收敛快慢，作出解x1随迭代次数变化的图像：
x1在两种迭代过程中的变化图像
x1
迭代次数
图1： 解x1随迭代次数变化图像
【结果分析】
由以上数据图表可知，x
(0)
=[0,0,0?0]?，b=[1,1,1?1]?，x(k?1)?x(k)
?
?10?5时，解
向量序列均是收敛的，都收敛到精确解，不过收敛速度不同。雅克比迭代法共迭代16次，高斯—赛德尔迭代法共迭代11次，因此可以得出高斯—赛德尔迭代法收敛速度更快。两种迭代法的收敛情况类似。迭代前几步变化较大，后面均在小范围内变化，与真实解已十分接近，直到满足相应的迭代误差条件时停止迭代过程。这从图1可以很直观地看出。
选取不同初始向量x和不同的方程组右端向量b,给定不同迭代误差要求，比较雅克比和高斯——赛德尔迭代法
① 选取相同初始向量x和不同的方程组右端向量b,给定相同迭代误差要求，用雅克比和高斯——赛德尔迭代法计算，比较结果：
表3：
、迭代误差相同，不同右端向量b时两种迭代方法的计算比较
(0)(0)
【结果分析】
由表格可知，仅右端向量b不同时，解向量序列均是收敛的，在误差允许范围内，
对相同b，两种方法得到解向量相同；不同b之间，解向量不同。对相同b，两种方法收敛速度不同。雅克比迭代法比高斯—赛德尔迭代法迭代次数多，收敛速
度慢，高斯—赛德尔迭代法收敛速度快。还可以看出，b向量越复杂，值越大，两种迭代方法的迭代次数均增加。 分析：具体的原因见后面的理论分析部分。
② 选取不同初始向量x给定相同的方程组右端向量b,迭代误差要求，用雅克比和高斯——赛德尔迭代法计算，比较结果：
表4：
不同，向量b、迭代误差相同时两种迭代方法的计算比较
(0)
【结果分析】
由表格可知，仅初始向量x不同时，解向量序列均是收敛的，在误差允许范围内解向量均相同，不过收敛速度不同。雅克比迭代法比高斯—赛德尔迭代法迭代次数多，收敛速度慢，高斯—赛德尔迭代法收敛速度快。还可知，初始向量x越大，偏离解越多，所需要的迭代次数越多，两种迭代方法的迭代次数均增加。 分析：这从直观上也是很容易理解的：当初始
距离解x越远时，所需要的修
(0)
(0)
正次数就越多，故迭代次数越多。因为对于同一个b，方程形式是相同的，故在收敛的前提下，不同的
最终将趋于同一个值，也就是方程的真实解。
③ 选取相同初始向量x、方程组右端向量b,给定不同迭代误差要求，用雅克比和高斯——赛德尔迭代法计算，比较结果：
表5：
、向量b相同、迭代误差不同时两种迭代方法的计算比较
(0)
【结果分析】
由表格可知，仅迭代误差要求不同时，解向量序列均是收敛的，即便给定的条件
相同，两种方法得到的解向量也略有不同，如表格中红色区域所示，说明迭代精度要求不高时，解略有偏差。对于相同条件，两种迭代法收敛速度也不同。雅克比迭代法比高斯—赛德尔迭代法迭代次数多，收敛速度慢，高斯—赛德尔迭代法收敛速度快。还可以看出，随着迭代误差要求提高，两种迭代方法的迭代次数均增加。
分析：误差条件变严格后，原来的迭代次数所对应的解已无法满足新的误差条件，故需要更多的迭代次数以使解满足当前的误差条件。
【理论分析】
①收敛性的分析
两种迭代方法在求解线性方程组时，都能写成下列形式：
从而有
进而得到下列结果：
（1）对于雅可比迭代法法，B=D?1(L?U)。在Matlab中运行eig(D_inv*(L+U)语句，得到B=D?1(L?U)的全部特征值为-0.4893，-0.4579，-0.4079，-0.3425，-0.2658，-0.1825，-0.0978，-0.0162，0.0577，0.1205，0.1677，0.1700，0.1819，0.1914，0.2088，0.2106，0.2298，0.2317，0.2445，0.2454。其中绝对值的最大的值为：0.4893＜1，即B的谱半径小于1，故对方程Ax=b进行雅克比迭代是收敛的。
（2）对于高斯—赛德尔迭代法，
。同样可以求得B的全部特征值
绝对值最大的值为0.1210＜1，即B的谱半径小于1，故对方程Ax=b进行高斯—赛德尔迭代是收敛的。
可见，上述分析与之前的运行结果是相一致的。 ②收敛速度分析
由于矩阵的谱半径不超过它的任一种范数，所以若敛，且有误差估计式
,则迭代公式收
可见，
越小，序列
收敛越快。
对于雅可比迭代法，B=D?1(L?U)。利用Matlab的norm命令可以求得B的∞-范数为0.5000。对于高斯—赛德尔迭代法，
。同样可以求得B的
∞-范数为0.2802。即高斯—赛德尔迭代法对应的序列
值更小，故与之对应的
收敛更快，从而对于一定的迭代误差条件，高斯—赛德尔迭代法比雅可
比迭代法法的迭代次数要少。
可见，上述分析与之前的运行结果是相一致的。
【拓展对比】
本题采用的是迭代法求解线性方程组。也可以通过直接法求解线性方程组，
如高斯消元法，LU分解等。可以在Matlab中利用x=A\b或者利用递推关系式进行求解。求出来的结果与用迭代法求出来的解略有出入。这是由求解过程中的误差所致。
【第(2)问求解】
取定方程组右端向量b=[1,2,3?20]?，同时取定初始向量x(0)=[0,0,0?0]?，并将迭代误差设为为：x(k?1)?x(k)
?
?10?5，将A的主对角线元素分别增长为6，
12，24，48，非主对角元素不变，每次用雅可比迭代法进行迭代计算，观察收敛速度的变化，得到结果如下面的表6：
表6： A主对角元素加倍时收敛速度的变化(b=
，x(0)=
【结果分析】
由表6可以明显看到，随着矩阵A对角元素的增大，对于雅可比迭代法，与之对应的
的值越来越小。在b和
相同以及迭代误差条件相同的前提
下，迭代次数越来越少，即收敛速度越来越快。
分析： 由前一问对收敛速度的分析可得， 由于矩阵的谱半径不超过它的任一种范数，所以若收敛，且有误差估计式
,则迭代公式
可见，
越小，序列
收敛越快。
对于雅可比迭代法，B=D?1(L?U)。利用Matlab的norm命令可以求得当 A的主对角元素分别为3，6，12，24，48时，相应的B的∞-范数分别为为0.5000，0.2500，0.1250，0.0625，0.03125。即随着A的主对角元素越来越大，相应的的值越来越小，从而序列
收敛越快，故迭代次数越来越少（在相同的b和
，
以及相同的误差条件的前提下）。
【本题小结】
(1)对于本题中的方程组Ax=b，选取不同的初始向量和不同的方程组右端向量b，所得到的迭代向量序列均是收敛的。这是由两种迭代法所对应的B的
谱半径均小于1所决定的。
(2)对于固定的和b，以及相同的迭代误差条件，两种迭代方法的迭代次数不同。高斯—赛德尔迭代法对应的迭代次数要比雅可比迭代法对应的迭代次数少。这是因为高斯—赛德尔迭代法对应的收敛更快。
(3)对于同一个b，不同的对应的迭代次数不同，且距离解x越远，所需的迭代次数越多。但最后得到的解是一样的。而对于不同的b，最后得到的解一般是不同的。
(3)对于同一个，相同b，给定不同的迭代误差要求，最后得到的解可能略有偏差。误差精度越高，所需的迭代次数越多。误差条件变严格后需要更多的迭代次数以使解满足当前的误差条件。
(5)通过数据表和图像还可以观察到，两种迭代法的收敛情况类似。迭代前几步变化较大，后面均在小范围内变化，与真实解已十分接近，直到满足相应的迭代误差条件时停止迭代过程。
(6)随着矩阵A对角元素的增大，在b和相同以及迭代误差条件相同的前提下，迭代次数越来越少，即收敛速度越来越快。这是由与之对应的的值越来越小所决定的。
值更小，
故与之对应的序列
【题目2】（课本习题第五章第5题）
设国民经济由农业、制造业和服务业三个部门构成，已知某年它们之间的投入产
出关系、外部需求、初始投入等如表7所示。
表7
（1） 如果今年对农业、制造业和服务业的外部需求分别为50、150、100亿元，
问这部门的总产出分别应为多少？
（2） 如果三个部门的外部需求分别增加一个单位，问它们的总产出分别增加多
少？
【第(1)问求解】
由表7得到三个部门的投入产出表如表8所示：
表8
xi
设有n个部门，记一定时期内第i个部门的总产出为投入为
xij
，其中对第j个部门的
，外部需求为
di
，则
xi??xij?di,i?1,2,...,n.
j?1n
根据投入系数
aij
的定义，有
xij?aijxj,i,j?1,2,...,n,
即
aij
是第j个部门的单位产出所需的第i个部门的投入。由此可得
xi??aijxj?di,i?1,2,...,n.
j?1n
，需求向量
记投入系数矩阵
A?(aij)n?n
，产出向量
x?(x1,x2,...,xn)T
d?(d1,d2,...,dn)T
，则上式可写作x?Ax?d，或
(I?A)x?d.
已知今年对农业、制造业和服务业的外部需求分别为50亿元，150亿元，100亿元，求这三个部门的总产出： 【matlab程序】
%------------------------------作业题5_5脚本M文件源程序---------------------------- clear all;clc;
A=[0.15 0.1 0.2;0.3 0.05 0.3;0.2 0.3 0]; % 设定系数矩阵A d=[50 150 100]'; % 给定外部需求 b=eye(3)-A;
x=b\d % 求解方程
dx=inv(b)
得到结果如下：
x ??139.2801, 267.6056, 208.1377?
T
?1.34590.25040.3443?
?dx??0.56341.26760.4930??
??0.43820.43041.2167??
即三个部门的总产出应分别为：139.2801、267.6056、208.1377（亿元）。
【第(2)问求解】
当三个部门的外部需求分别增加1个单位时，求它们总产出相应增加的量： 由(I?A)x?d可得
x?(I?A)?1d.
表明总产出x对外部需求d是线性的，所以当d增加1个单位（记作?d）时，x
?1
?x?(I?A)?d. 的增量为
于是从dx=inv(b)矩阵中便可得三个部门的外部需求分别增加一个单位时，各部
门总产出分别增加量。
?1.34590.25040.3443?
?dx??0.56341.26760.4930??
??0.43820.43041.2167??
结果如下：
T?1
?d?(1,0,0)(I?A)?x若农业的外部需求增加1个单位，，为的第1列
?x=?1.3459，0.5634,0.4382?；即当农业的需求增加1个单位时，农业、制造业和服务业的总产出应分别增加1.3459, 0.5634, 0.4382单位。
若制造业的外部需求增加1个单位，?d?(0,1,0)，?x为(I?A)的第2列
T
?1
T
?x=?0.2504,1.2676,0.4304?；即当制造业的需求增加1个单位时，农业、制造业和服务业的总产出应分别增加0.2504,1.2676，0.4304单位。
若服务业的外部需求增加1个单位，?d?(0,0,1)，?x为(I?A)的第3列
T
?1
T
?x=?0.3443,0.4930,1.2167?；即当服务业的需求增加1个单位时，农业、制造业和服务业的总产出应分别增加0.3443,0.4930，1.2167单位。
T
【本题小结】
本题是对实际问题的数学建模，模型较简单，问题最后转化为线性方程组求解问题。因为方程并不复杂，我采用直接法求解。很快得到了答案。
【题目3】（课本习题第五章第7题）
输电网络：一种大型输电网络可简化为图5.5所示电路，其中R1，R2，…，Rn表示负载电阻，r1，r2，…，rn表示线路内阻，I1，I2，…，In表示负载上的电流。设电源电压为V。
(1)列出求各负载电流I1，I2，…，In的方程；
(2)设R1=R2=…=Rn=R，r1=r2=…=rn=r，在r=1，R=6，V=18，V=18，n=10的情况下求I1，I2，…，In及总电流I0。
图5.5
【第（1）问求解】
问题分析及模型建立：
I?Ii?Ii?1??Inr
通过对以上电路的分析可得，流过电阻i的电流ri。对于每一个
由Ri,ri+1,Ri?1组成的闭合回路，根据基尔霍夫回路电压法，有
-IiRi+Iri+1ri+1?Ii?1Ri?1?0。对于由电源V,R1,r1组成的回路，有
-V+I1R1?Ir1r1?0
。
因此，可列方程如下：
?I1R1?(I1+I2+?In)r1?V?0
?-IR?IR?(I?I??I)r?0112223n2??
?-I2R2?I3R3?(I3?I4??In)r3?0 ?????-In-1Rn?1?InRn?Inrn?0
整理后可得，
R11?(1+)I+I+?I?V12n?r1r1?
R2?R1
-I?(1?)I2?I3??In=0?r1
r?22
?
RR?-2I?(1?3)I?I??I=0
34n
?r32r3?
?-Rn?1I?(1?Rn)I?0
n-1n
??rn
rn变成矩阵的形式，即AX?b，其中，
??1+R1 1 1 … … 1 ?r1
???R1R2?r 1?2r 1 … ? 12 A=?? ?R
2?r ? 3
?? ? … ? ? ? ??? ?Rn?1r 1?Rnnrn?I?V??1?
?I?
?2?r?X???
1??I?0?3?
???b??0?
??，???。 ??In-1??????0??I?n?????
?
最后列出求各负载电流I1，I2，…，In的方程如下：
??
?
?
?????????
??
?R1?I1+ 1 1 … … 1 ?r??1??V?
1
?????r1?
R2?R1?????
??r 1?r 1 … ? 1 ??I2??0?
2?2?????
??????R2
? ? ???I3???0? r3
????????? ? … ? ? ????????In-1??0?? ?????
?????Rn?1Rn?????
1? ??In??0?? ?rnrn??
【第（2）问求解】
由题意可得， R1=R2=…=Rn=R，r1=r2=…=rn=r，r=1，R=6，V=18， n=10因此，
矩阵A变作
?7 1 1 … 1 ?
?0??6 7 1 … 1 ?
???
? ?6 ? ? 矩阵b变作 b=?0 A=??? ? ? ????? ?
????0
? ?6 7?
?7 1 1 … 1 ??I1??18???6 7 1 … 1 ??I??0????2???
整理方程得到 ? ?6 ? ??I3?=?0?
?????? ? ? ????????? ??In-1??0??????? ?6 7????In???0?
?18?
??? ????
其中n=10。
【Matlab程序】
%------------------------------作业题5_7脚本M文件源程序---------------------------- clear all;clc; n=10;
A1=sparse(1:n,1:n,7,n,n); A2=sparse(2:n,1:n-1,-6,n,n); A3=triu(ones(10,10),1);
A=A1+A2+A3; % 设定系数矩阵A b=[18,zeros(1,9)]'; % 设定列向量b x=[A\b]' % 求解各支路电流值 I0=sum(x) % 求解电流I0
得到如下结果：
表9：r=1，R=6，V=18，V=18，n=10时I1，I2，…，In，I0的值
所以要求的I1，I2，…，In及总电流I0值如表9所示。
【本题小结】
本题是电路问题的求解。可以采用差分方程求解，也可以用线性方程组求解。我采用线性方程组求解得到答案。因为问题中矩阵A并不复杂，因此采用直接法求解方程组，模型建好后很快便得到了答案。
【实验心得、体会】
本次实验作为数学实验课程的第三次作业，让我们掌握了线性方程组的数值解法相关基础知识，同时我进一步了解了matlab。本次实验第一个题目难度较大其余两题目都比较轻松。通过第一题，我深刻理解了雅克比和高斯——赛德尔迭代法，增长了许多知识。实验中最令我困惑的是为什么很简单的一个矩阵求解（线性方程组求解）问题要用雅克比和高斯——赛德尔迭代法求解，比如第一问中已知矩阵A与b，直接用x=A\b便可得解，我们却要绕一个大弯得到x的解，后来我逐渐明白了，直接法适用于求解中小型（n
我尝试着对每个问题用学到的不同的方法去解决（如对线性方程组分别用直接法中的高斯消元法和LU分解以进行求解，并和两种迭代方法进行比较），同时尝试着从不同的方面和角度去分析，并尽可能地发掘出更多值得去思考和分析的东西，在这个过程中不仅巩固了此次所学到的了线性方程组的数值解等相关知识，也不仅进一步熟练了Matlab的操作，更重要的是通过这次作业进行了一次思维上的极大锻炼，进一步熟悉了通过数学建模思想和数学软件解决一些和现实密切相关的问题。
总之，本次实验收获挺大的。
转载请注明出处范文大全网 » 金属—绝缘层—半导体结构泊松
