数值模拟导论 ——第十一讲
牛顿法实例学习 ——模拟图像滤波器
雅可比 .怀特
Thanks to Deepak Ramaswamy, Andrew Lumsdaine,
Jaime Peraire, MichalRewienski, and Karen Veroy
SMA-HPC ?2003 MIT
图像分割实例
牛顿迭代法
Gershgorin圆理论
弧长连续
概要：
简单的滤波器
SMA-HPC ?2003 MIT
电路图
光滑输
出
图像输入
非线性滤波器
SMA-HPC ?2003 MIT
电路图
光滑输出
出
图像输
入
非线性滤波器
SMA-HPC ?2003 MIT
电路图
非线性电阻器基本方程
SMA-HPC ?2003 MIT
2
()
()
1
v
v
iv
e
βγα
α
??
=
+
电
流
电流
改变β值
SMA-HPC ?2003 MIT
存在的问题
SMA-HPC ?2003 MIT
z用什么方程式表示？
——节点－分支法或节点法
z用什么牛顿法？
——标准的，阻尼的，或连续的？
——哪一种连续法？
z用什么线性求解器？
电流
——高斯消去法还是 Krylove法？
——Krylove法能快速收敛么？
——方程式的选择和牛顿法选择是否有冲突？
牛顿迭代法
SMA-HPC ?2003 MIT
基本算法
嵌套迭代
0
x
＝初始值，
0k =
重复 { ( ) ( )
() ()
11
11
,
1
kk
F
kk k k
F
kkkk
Fx J x
J xx Fx x
xx x
kk
α
+ +
++
?=? ?
=+?
=+
计算
求解（利用成组编码法）
解方程 求得
}
( )
11
,
kk
xFx
++
?直到 足够小为止
成组编码法应该取多大精度？
SMA-HPC ?2003 MIT
求解精度要求
成组编码法 l步之后
( )
N
N
1,
()
kk k kl
F
l
Jx x Fx γ
+
?
?=? +
步成组编码
成组编码
后的牛顿
余量
()
() ()
1
2
,
a) ( )
) ( Lipschitz Cont)
) ( ) ( )
k
F
FF
kl k
Jx
bJxJy lxy
cCFx
β
γ
?
≤
?≤?
≤
如果
反之不成立
导数为
更精确接近收敛
那么
牛顿迭代二次收敛
收敛性证明
SMA-HPC ?2003 MIT
通过定义牛顿迭代
( )()
( )
1
1,kk k kkl
F
xxJxFxγ
?
+
=? +
接近牛顿方向
多维均值迭代
() () ()( )
2
2
F
l
F xFyJyxy xy?? ?≤?
综合上述两式得
( ) () () () () () ()
2
11
1,,
2
kkkkkkl kkkl
FF F
l
Fx Fx J x J x Fx J x Fxγγ
??
+
??
?? +≤ +
??
??
消去上式中雅可比矩阵及其逆得
()()() () ()
2
1
1, ,
2
kkkkl kkkl
F
l
Fx Fx Fx J x Fxγγ
?
+
?++≤ +
化简并利用三角不等式得
() () ()
2
1
1,,
2
kkkklkl
F
l
Fx J x Fx γ γ
?
+
≤++
SMA-HPC ?2003 MIT
利用雅可比矩阵的极限和三角不等式
() ()
2,
2
2
1,
1
22
kl
kk kl
l
Fx Fx
βγ
β
γ
+
??
??
≤++
??
在迭代求解有误时利用范围限制
() () ()
2,
2
22
1
1
22
kl
kk k
l
Fx Fx CFx
βγ
β
+
??
??≤++
??
综合上述得
() ()
2,
2
2
1
1
22
kl
kk
l
Fx C Fx
βγ
β
+
??
??
??
??
≤++
??
易于限制起值的范围
牛顿迭代法
SMA-HPC ?2003 MIT
矩阵－释放思想
考虑应用成组编码法到牛顿迭代方程中
( )()
1kk k
F
J xx Fx
+
?=?
在每一次迭代中成组编码法形成一矩阵矢量值
()
( )
()
( )
1
l
kl k k
F
Jxp Fx p Fxε
ε
≈+?
可能采用牛顿成组编码法而无需雅可比矩阵
需要选出一个恰当的ε值
Gerschgorin 圆定理
SMA-HPC ?2003 MIT
定理陈述
给出矩阵
1,1 1,
,1 ,
N
NNN
mm
M
mm
? ?
? ?
=
? ?
? ?
? ?
"
# % #
"
对于每一M的特征值都存在一个i，1<i<N 使得
,,ii i j
ji
mmλ
≠
?≤
∑
可以得出特征值包含在Gerschgorin圆集合中
Gerschgrion图
SMA-HPC ?2003 MIT
接地电阻丝
SMA-HPC ?2003 MIT
节点矩阵
2.1 -1
-1 2.1
-1
-1 2.1
M
??
??
??
??
%
节点方程式
SMA-HPC ?2003 MIT
N
2 -1
-1 2
-1
-1 2
M
? ?
? ?
? ?
? ?
? ?
? ?
%
节点方程式
连续方案
SMA-HPC ?2003 MIT
基本概念
一般设置
求解 ( )( )
,0Fxλλ=
其中
( )( )
0,0 0Fx =
易于求解 开始时连续
( )( )
1,1 ( )F xFx=
终止时连续
( )
x λ
充分光滑 难以确保！
连续定律
SMA-HPC ?2003 MIT
基本概念
模板算法
( )( ) ()( )
() ()
()
()()
0
0,0 0, 0
0.01,
1 {
,0
,
=+ , 2
1
=+
2
}
F xxx
xx
Fx
xx
λ
δλ λ δλ
λ
λλ
λλ
λλ
λλδλδλ δλ
δλ δλ λ λ δλ
==
==
<
=
=
=
=
前一个
前一个
前一个
前一个
求解
当
试用牛顿法求解
若牛顿收敛 ＝
否则 ，
连续方案
SMA-HPC ?2003 MIT
雅克比迭代定律
定律描述
( )( ) ( )( ) ( ) ( )
,1Fx x xλ λλλ λλ=+?
观察
( )( ) ( )
0 0 ,0 0 0
((0),0)
Fx x
Fx
I
x
λ = ==
?
=
?
问题易于求解且雅克比矩阵确定非奇异。
( )( ) ( )( )
()()
1 1 ,1 1
1
((0),0)
Fx Fx
Fx
Fx
xx
λ ==
?
?
=
??
重新回到初始的问题和初始雅克比矩阵
基本算法
SMA-HPC ?2003 MIT
( )( ) ( ) ( )
() ()
()
()()
0
0 ,0 0, 0
0.01,
1 {
?
,0
, = + , 2
1
=+
2
}
Fx x x
xx
Fx
xx
λ
δλ λ δλ
λ
λλ
λλ
λ λ λλδλδλ δλ
δλ δλ λ λ δλ
==
==
<
=+
=
=
=
前一个
前一个
前一个
前一个
求解
当
试用牛顿法求解
若牛顿收敛 ＝
否则 ，
SMA-HPC ?2003 MIT
弧长步？
SMA-HPC ?2003 MIT
需求解λ
()()
2
2
2
2
(, ) 0
0
Fx
xx arc
λ
λλ λ
=
? +? ? =
前一个 前一个
弧长步牛顿法
SMA-HPC ?2003 MIT
()
()
()
( )
()
()
()()
1
1
22
2
2
,
,
2
2,
,
kk
kk
kk
kk
T
k
k
kk
kk
Fx
Fx
xx
x
xx
Fx
xx arc
λ
λ
λ
λλ
λλ
λ
λ
λλ λ
+
+
?
? ?
?
?
? ?
??
?
=
??
? ?
?
?
?
? ?
?
?
? ?
?
??
?
??
?
??? ?
??
前一个
前一个
前一个 前一个
弧长转折点
SMA-HPC ?2003 MIT
( )
()
()
( )
()
,
,
2
2,
kk
kk
T
k
k
Fx
Fx
x
xx
λ
λ
λ
λλ
λ
?
? ?
?
?
?
?
??
?
?
?
?
? ?
?
前一个
前一个
矩阵左上角模块为奇数
小结：
SMA-HPC ?2003 MIT
图像过滤实例
——大的非线性系统方程
——在选择数值方法中检验问题
牛顿迭代法
——无需精确求解迭代方程
Gershgorin圆定理
——有时给出特征值的有利范围
弧长连续
SMA-HPC ?2003 MIT
SMA-HPC ?2003 MIT
SMA-HPC ?2003 MIT
SMA-HPC ?2003 MIT
