在数学 上,以皮埃尔-西蒙·拉普拉斯 命名的拉普拉斯方法 是用于得出下列积分 形式的近似解的方法:
∫
a
b
e
M
f
(
x
)
d
x
{\displaystyle \int _{a}^{b}\!e^{Mf(x)}\,dx}
其中的 ƒ (x ) 是一个二次可微 函数 , M 是一个很大的数,而积分边界点 a 与 b 则允许为无限大。此外,函数 ƒ (x ) 在此积分范围内的 全域极大值 所在处必须是唯一的并且不在边界点上。则它的近似解可以写为:
∫
a
b
e
M
f
(
x
)
d
x
≈
2
π
M
|
f
″
(
x
0
)
|
e
M
f
(
x
0
)
as
M
→
∞
.
{\displaystyle \int _{a}^{b}\!e^{Mf(x)}\,dx\approx {\sqrt {\frac {2\pi }{M\left|f''(x_{0})\right|}}}e^{Mf(x_{0})}{\text{ as }}M\to \infty .\,}
其中的 x 0 为极大值所在处。这方法最早是拉普拉斯在 (1774, pp. 366–367) 所发表。(待考查)
拉普拉斯方法的想法概论
1. 相对误差
首先,我们得先知道的是,这里所谓的近似解是以 相对误差 在讲,而不是指 绝对误差 。因此,令
s
=
2
π
M
|
f
″
(
x
0
)
|
{\displaystyle s={\sqrt {\frac {2\pi }{M\left|f''(x_{0})\right|}}}}
,
此 s 很显然的当 M 很大时为一个微小的数,则上述的积分可以写为
∫
a
b
e
M
f
(
x
)
d
x
=
s
e
M
f
(
x
0
)
1
s
∫
a
b
e
M
(
f
(
x
)
−
f
(
x
0
)
)
d
x
=
s
e
M
f
(
x
0
)
∫
(
a
−
x
0
)
/
s
(
b
−
x
0
)
/
s
e
M
(
f
(
s
y
+
x
0
)
−
f
(
x
0
)
)
d
y
{\displaystyle {\begin{aligned}\int _{a}^{b}\!e^{Mf(x)}\,dx&=se^{Mf(x_{0})}{\frac {1}{s}}\int _{a}^{b}\!e^{M(f(x)-f(x_{0}))}\,dx\\&=se^{Mf(x_{0})}\int _{(a-x_{0})/s}^{(b-x_{0})/s}\!e^{M(f(sy+x_{0})-f(x_{0}))}\,dy\end{aligned}}}
因此,我们的相对误差很显然的为
|
∫
(
a
−
x
0
)
/
s
(
b
−
x
0
)
/
s
e
M
(
f
(
s
y
+
x
0
)
−
f
(
x
0
)
)
d
y
−
1
|
{\displaystyle \left|\int _{(a-x_{0})/s}^{(b-x_{0})/s}e^{M(f(sy+x_{0})-f(x_{0}))}dy-1\right|}
现在,让我们把积分分为
y
∈
[
−
D
y
,
D
y
]
{\displaystyle y\in [-D_{y},D_{y}]}
的与余下的部分。
基于上述四点,就有办法证明拉普拉斯方法的可靠性。而 Fog(2008) 又将此方法推广到任意精确。 ***待考查***
此方法的正式表述与证明:
假设
f
(
x
)
{\displaystyle f(x)}
是一个在
x
0
∈
(
a
,
b
)
{\displaystyle x_{0}\in (a,b)}
这点满足 (1)
f
(
x
0
)
=
max
[
a
,
b
]
f
(
x
)
{\displaystyle f(x_{0})=\max _{[a,b]}f(x)}
,(2)唯一全域最大,(3)
x
0
{\displaystyle x_{0}}
附近为二阶可微且
f
″
(
x
0
)
<
0
{\displaystyle f''(x_{0})<0}
(4) 当拉普拉斯方法的积分范围为无限大时,此积分会收敛,
则,
lim
n
→
+
∞
(
∫
a
b
e
n
f
(
x
)
d
x
(
e
n
f
(
x
0
)
2
π
n
(
−
f
″
(
x
0
)
)
)
)
=
1
{\displaystyle \lim _{n\to +\infty }\left({\frac {\int _{a}^{b}e^{nf(x)}\,dx}{\left(e^{nf(x_{0})}{\sqrt {\frac {2\pi }{n(-f''(x_{0}))}}}\right)}}\right)=1}
证明:
下限 :
在
f
″
{\displaystyle f''}
要求连续的条件底下,给定任意大于0的
ε
{\displaystyle \varepsilon }
就能找到一个
δ
>
0
{\displaystyle \delta >0}
使得在
|
x
0
−
x
|
<
δ
{\displaystyle |x_{0}-x|<\delta }
这区间内所有的
f
″
(
x
)
≥
f
″
(
x
0
)
−
ε
.
{\displaystyle f''(x)\geq f''(x_{0})-\varepsilon .}
。 由 泰勒展开 我们可以得知,在
x
∈
(
x
0
−
δ
,
x
0
+
δ
)
{\displaystyle x\in (x_{0}-\delta ,x_{0}+\delta )}
区间内,
f
(
x
)
≥
f
(
x
0
)
+
1
2
(
f
″
(
x
0
)
−
ε
)
(
x
−
x
0
)
2
{\displaystyle f(x)\geq f(x_{0})+{\frac {1}{2}}(f''(x_{0})-\varepsilon )(x-x_{0})^{2}}
。
这样,我们就得到本方法的积分下限了:
∫
a
b
e
n
f
(
x
)
d
x
≥
∫
x
0
−
δ
x
0
+
δ
e
n
f
(
x
)
d
x
≥
e
n
f
(
x
0
)
∫
x
0
−
δ
x
0
+
δ
e
n
2
(
f
″
(
x
0
)
−
ε
)
(
x
−
x
0
)
2
d
x
=
e
n
f
(
x
0
)
1
n
(
−
f
″
(
x
0
)
+
ε
)
∫
−
δ
n
(
−
f
″
(
x
0
)
+
ε
)
δ
n
(
−
f
″
(
x
0
)
+
ε
)
e
−
1
2
y
2
d
y
{\displaystyle \int _{a}^{b}e^{nf(x)}\,dx\geq \int _{x_{0}-\delta }^{x_{0}+\delta }e^{nf(x)}\,dx\geq e^{nf(x_{0})}\int _{x_{0}-\delta }^{x_{0}+\delta }e^{{\frac {n}{2}}(f''(x_{0})-\varepsilon )(x-x_{0})^{2}}\,dx=e^{nf(x_{0})}{\sqrt {\frac {1}{n(-f''(x_{0})+\varepsilon )}}}\int _{-\delta {\sqrt {n(-f''(x_{0})+\varepsilon )}}}^{\delta {\sqrt {n(-f''(x_{0})+\varepsilon )}}}e^{-{\frac {1}{2}}y^{2}}\,dy}
其中最后一个等号来自于变数变换:
y
=
n
(
−
f
″
(
x
0
)
+
ε
)
(
x
−
x
0
)
{\displaystyle y={\sqrt {n(-f''(x_{0})+\varepsilon )}}(x-x_{0})}
。请记得
f
″
(
x
0
)
<
0
{\displaystyle f''(x_{0})<0}
,所以我们才会对它的负号取开根号。
接着让我们对上面的不等式两边同除以
e
n
f
(
x
0
)
2
π
n
(
−
f
″
(
x
0
)
)
{\displaystyle e^{nf(x_{0})}{\sqrt {\frac {2\pi }{n(-f''(x_{0}))}}}}
并且对
n
{\displaystyle n}
取极限,则
lim
n
→
+
∞
(
∫
a
b
e
n
f
(
x
)
d
x
(
e
n
f
(
x
0
)
2
π
n
(
−
f
″
(
x
0
)
)
)
)
≥
lim
n
→
+
∞
1
2
π
∫
−
δ
n
(
−
f
″
(
x
0
)
+
ε
)
δ
n
(
−
f
″
(
x
0
)
+
ε
)
e
−
1
2
y
2
d
y
−
f
″
(
x
0
)
−
f
″
(
x
0
)
+
ε
=
−
f
″
(
x
0
)
−
f
″
(
x
0
)
+
ε
{\displaystyle \lim _{n\to +\infty }\left({\frac {\int _{a}^{b}e^{nf(x)}\,dx}{\left(e^{nf(x_{0})}{\sqrt {\frac {2\pi }{n(-f''(x_{0}))}}}\right)}}\right)\geq \lim _{n\to +\infty }{\frac {1}{\sqrt {2\pi }}}\int _{-\delta {\sqrt {n(-f''(x_{0})+\varepsilon )}}}^{\delta {\sqrt {n(-f''(x_{0})+\varepsilon )}}}e^{-{\frac {1}{2}}y^{2}}\,dy{\sqrt {\frac {-f''(x_{0})}{-f''(x_{0})+\varepsilon }}}={\sqrt {\frac {-f''(x_{0})}{-f''(x_{0})+\varepsilon }}}}
既然上式是对任意大于0的
ε
{\displaystyle \varepsilon }
都对,所以我们得到:
lim
n
→
+
∞
(
∫
a
b
e
n
f
(
x
)
d
x
(
e
n
f
(
x
0
)
2
π
n
(
−
f
″
(
x
0
)
)
)
)
≥
1
{\displaystyle \lim _{n\to +\infty }\left({\frac {\int _{a}^{b}e^{nf(x)}\,dx}{\left(e^{nf(x_{0})}{\sqrt {\frac {2\pi }{n(-f''(x_{0}))}}}\right)}}\right)\geq 1}
请注意,上面的证明也适用于
a
=
−
∞
{\displaystyle a=-\infty }
或
b
=
∞
{\displaystyle b=\infty }
或双双跑到正负无限大的情况
上限 :
证明上限的部分其实和证明下限的部分很像,但是会较麻烦。再一次,我们取
ε
>
0
{\displaystyle \varepsilon >0}
,不过,不再是任意而是多了一个要求
ε
{\displaystyle \varepsilon }
得够小以致于
f
″
(
x
0
)
+
ε
<
0
{\displaystyle f''(x_{0})+\varepsilon <0}
。接着,就如同之前的证明,因着
f
″
{\displaystyle f''}
被要求连续,并且根据 泰勒定理 我们会发现总存在一个
δ
>
0
{\displaystyle \delta >0}
以至于在
|
x
−
x
0
|
<
δ
{\displaystyle |x-x_{0}|<\delta }
区间里,
f
(
x
)
≤
f
(
x
0
)
+
1
2
(
f
″
(
x
0
)
+
ε
)
(
x
−
x
0
)
2
{\displaystyle f(x)\leq f(x_{0})+{\frac {1}{2}}(f''(x_{0})+\varepsilon )(x-x_{0})^{2}}
总是可以成立。 最后,在我们的假设里 (假设
a
,
b
{\displaystyle a,b}
都是有限值) ,由于
x
0
{\displaystyle x_{0}}
是全域最大所在处,我们总可以找到一个
η
>
0
{\displaystyle \eta >0}
使得所有在
|
x
−
x
0
|
≥
δ
{\displaystyle |x-x_{0}|\geq \delta }
这区间里,
f
(
x
)
≤
f
(
x
0
)
−
η
{\displaystyle f(x)\leq f(x_{0})-\eta }
总是成立。
现在,万事俱备,东风就在下面啦:
∫
a
b
e
n
f
(
x
)
d
x
≤
∫
a
x
0
−
δ
e
n
f
(
x
)
d
x
+
∫
x
0
−
δ
x
0
+
δ
e
n
f
(
x
)
d
x
+
∫
x
0
+
δ
b
e
n
f
(
x
)
d
x
≤
(
b
−
a
)
e
n
(
f
(
x
0
)
−
η
)
+
∫
x
0
−
δ
x
0
+
δ
e
n
f
(
x
)
d
x
{\displaystyle \int _{a}^{b}e^{nf(x)}\,dx\leq \int _{a}^{x_{0}-\delta }e^{nf(x)}\,dx+\int _{x_{0}-\delta }^{x_{0}+\delta }e^{nf(x)}\,dx+\int _{x_{0}+\delta }^{b}e^{nf(x)}\,dx\leq (b-a)e^{n(f(x_{0})-\eta )}+\int _{x_{0}-\delta }^{x_{0}+\delta }e^{nf(x)}\,dx}
≤
(
b
−
a
)
e
n
(
f
(
x
0
)
−
η
)
+
e
n
f
(
x
0
)
∫
x
0
−
δ
x
0
+
δ
e
n
2
(
f
″
(
x
0
)
+
ε
)
(
x
−
x
0
)
2
d
x
≤
(
b
−
a
)
e
n
(
f
(
x
0
)
−
η
)
+
e
n
f
(
x
0
)
∫
−
∞
+
∞
e
n
2
(
f
″
(
x
0
)
+
ε
)
(
x
−
x
0
)
2
d
x
{\displaystyle \leq (b-a)e^{n(f(x_{0})-\eta )}+e^{nf(x_{0})}\int _{x_{0}-\delta }^{x_{0}+\delta }e^{{\frac {n}{2}}(f''(x_{0})+\varepsilon )(x-x_{0})^{2}}\,dx\leq (b-a)e^{n(f(x_{0})-\eta )}+e^{nf(x_{0})}\int _{-\infty }^{+\infty }e^{{\frac {n}{2}}(f''(x_{0})+\varepsilon )(x-x_{0})^{2}}\,dx}
≤
(
b
−
a
)
e
n
(
f
(
x
0
)
−
η
)
+
e
n
f
(
x
0
)
2
π
n
(
−
f
″
(
x
0
)
−
ε
)
{\displaystyle \leq (b-a)e^{n(f(x_{0})-\eta )}+e^{nf(x_{0})}{\sqrt {\frac {2\pi }{n(-f''(x_{0})-\varepsilon )}}}}
如果我们对上面的不等式两边皆除以
e
n
f
(
x
0
)
2
π
n
(
−
f
″
(
x
0
)
)
{\displaystyle e^{nf(x_{0})}{\sqrt {\frac {2\pi }{n(-f''(x_{0}))}}}}
并且顺便取极限的话,会得到:
lim
n
→
+
∞
(
∫
a
b
e
n
f
(
x
)
d
x
(
e
n
f
(
x
0
)
2
π
n
(
−
f
″
(
x
0
)
)
)
)
≤
lim
n
→
+
∞
(
(
b
−
a
)
e
−
η
n
n
(
−
f
″
(
x
0
)
)
2
π
+
−
f
″
(
x
0
)
−
f
″
(
x
0
)
−
ε
)
=
−
f
″
(
x
0
)
−
f
″
(
x
0
)
−
ε
{\displaystyle \lim _{n\to +\infty }\left({\frac {\int _{a}^{b}e^{nf(x)}\,dx}{\left(e^{nf(x_{0})}{\sqrt {\frac {2\pi }{n(-f''(x_{0}))}}}\right)}}\right)\leq \lim _{n\to +\infty }\left((b-a)e^{-\eta n}{\sqrt {\frac {n(-f''(x_{0}))}{2\pi }}}+{\sqrt {\frac {-f''(x_{0})}{-f''(x_{0})-\varepsilon }}}\right)={\sqrt {\frac {-f''(x_{0})}{-f''(x_{0})-\varepsilon }}}}
再一次,因为
ε
{\displaystyle \varepsilon }
可以取任意大于0的值,所以我们得到了上限了:
lim
n
→
+
∞
(
∫
a
b
e
n
f
(
x
)
d
x
(
e
n
f
(
x
0
)
2
π
n
(
−
f
″
(
x
0
)
)
)
)
≤
1
{\displaystyle \lim _{n\to +\infty }\left({\frac {\int _{a}^{b}e^{nf(x)}\,dx}{\left(e^{nf(x_{0})}{\sqrt {\frac {2\pi }{n(-f''(x_{0}))}}}\right)}}\right)\leq 1}
把上限与下限两个证明同时考虑,整个证明就完成了。
注意,关于上限的证明,很明显的当我们把它应用在
a
{\displaystyle a}
或
b
{\displaystyle b}
为正负无限大时,该上限证明会失败。那怎么办呢?我们需要再多假设一些东西。一个充分但非必要的假设是:当
n
=
1
{\displaystyle n=1}
时,此积分
∫
a
b
e
n
f
(
x
)
d
x
{\displaystyle \int _{a}^{b}e^{nf(x)}\,dx}
为有限值,并且上面所说的
η
{\displaystyle \eta }
是存在的 (注意,当
[
a
,
b
]
{\displaystyle [a,b]}
区间是无限的时候,这假设是必要的) 。整个证明过程就如同先前所显示的那样,只不过下列的积分部分要做点改变:
∫
a
x
0
−
δ
e
n
f
(
x
)
d
x
+
∫
x
0
+
δ
b
e
n
f
(
x
)
d
x
{\displaystyle \int _{a}^{x_{0}-\delta }e^{nf(x)}\,dx+\int _{x_{0}+\delta }^{b}e^{nf(x)}\,dx}
必须利用上述的假设,而改为:
∫
a
x
0
−
δ
e
n
f
(
x
)
d
x
+
∫
x
0
+
δ
b
e
n
f
(
x
)
d
x
≤
∫
a
b
e
f
(
x
)
e
(
n
−
1
)
(
f
(
x
0
)
−
η
)
d
x
=
e
(
n
−
1
)
(
f
(
x
0
)
−
η
)
∫
a
b
e
f
(
x
)
d
x
{\displaystyle \int _{a}^{x_{0}-\delta }e^{nf(x)}\,dx+\int _{x_{0}+\delta }^{b}e^{nf(x)}\,dx\leq \int _{a}^{b}e^{f(x)}e^{(n-1)(f(x_{0})-\eta )}\,dx=e^{(n-1)(f(x_{0})-\eta )}\int _{a}^{b}e^{f(x)}\,dx}
以取代先前会得到的
(
b
−
a
)
e
n
(
f
(
x
0
)
−
η
)
{\displaystyle (b-a)e^{n(f(x_{0})-\eta )}}
,这样的话,当我们除以
e
n
f
(
x
0
)
2
π
n
(
−
f
″
(
x
0
)
)
{\displaystyle e^{nf(x_{0})}{\sqrt {\frac {2\pi }{n(-f''(x_{0}))}}}}
,就会改得到如下的结果:
e
(
n
−
1
)
(
f
(
x
0
)
−
η
)
∫
a
b
e
f
(
x
)
d
x
e
n
f
(
x
0
)
2
π
n
(
−
f
″
(
x
0
)
)
=
e
−
(
n
−
1
)
η
n
e
−
f
(
x
0
)
∫
a
b
e
f
(
x
)
d
x
−
f
″
(
x
0
)
2
π
{\displaystyle {\frac {e^{(n-1)(f(x_{0})-\eta )}\int _{a}^{b}e^{f(x)}\,dx}{e^{nf(x_{0})}{\sqrt {\frac {2\pi }{n(-f''(x_{0}))}}}}}=e^{-(n-1)\eta }{\sqrt {n}}e^{-f(x_{0})}\int _{a}^{b}e^{f(x)}\,dx{\sqrt {\frac {-f''(x_{0})}{2\pi }}}}
这样的话,当我们取
n
→
∞
{\displaystyle n\rightarrow \infty }
时,上式的值就会趋近于
0
{\displaystyle 0}
。而剩下的部分的证明就还是如同原先的证明,不做改变。
再强调一次,这里我们多加给无限大积分范围的情形的条件,是充分,但非必要。不过,这样的条件已经可以适用在许多情形了(但非全部)。 这考虑条件简单来讲就是积分区间得是被良好定义的(即不能是无限大的),并且被积函数在
x
0
{\displaystyle x_{0}}
必须是真的极大 (意即
η
>
0
{\displaystyle \eta >0}
必须真的存在) ;如果这积分区间是无限大的话,要求
n
=
1
{\displaystyle n=1}
时的此拉普拉斯方法所用的积分值要为有限并非必要的,其实只要当
n
{\displaystyle n}
大于某数时,此积分值会是有限的即可。
其他形式
有时拉普拉斯方法也会被写成其他形式,如:
∫
a
b
h
(
x
)
e
M
g
(
x
)
d
x
≈
2
π
M
|
g
″
(
x
0
)
|
h
(
x
0
)
e
M
g
(
x
0
)
as
M
→
∞
{\displaystyle \int _{a}^{b}\!h(x)e^{Mg(x)}\,dx\approx {\sqrt {\frac {2\pi }{M|g''(x_{0})|}}}h(x_{0})e^{Mg(x_{0})}{\text{ as }}M\to \infty \,}
其中
h
{\displaystyle h}
为正 (好像不必要)。
重要的是,这方法精确度与函数
g
(
x
)
{\displaystyle g(x)}
和
h
(
x
)
{\displaystyle h(x)}
有关。 [ 1] ***待考查***
相对误差的推导
首先,由于极大值所在设为
x
0
=
0
{\displaystyle x_{0}=0}
并不影响计算结果,所以以下的推导皆如此假设。此外,我们想要的是相对误差
|
R
|
{\displaystyle \left|R\right|}
,所以
∫
a
b
h
(
x
)
e
M
g
(
x
)
d
x
=
h
(
0
)
e
M
g
(
0
)
s
∫
a
/
s
b
/
s
h
(
x
)
h
(
0
)
e
M
[
g
(
s
y
)
−
g
(
0
)
]
d
y
⏟
1
+
R
{\displaystyle \int _{a}^{b}\!h(x)e^{Mg(x)}\,dx=h(0)e^{Mg(0)}s\underbrace {\int _{a/s}^{b/s}{\frac {h(x)}{h(0)}}e^{M\left[g(sy)-g(0)\right]}dy} _{1+R}}
其中
s
≡
2
π
M
|
g
″
(
0
)
|
{\displaystyle s\equiv {\sqrt {\frac {2\pi }{M\left|g''(0)\right|}}}}
,所以,若我们令
A
≡
h
(
s
y
)
h
(
0
)
e
M
[
g
(
s
y
)
−
g
(
0
)
]
{\displaystyle A\equiv {\frac {h(sy)}{h(0)}}e^{M\left[g(sy)-g(0)\right]}}
且
A
0
≡
e
−
π
y
2
{\displaystyle A_{0}\equiv e^{-\pi y^{2}}}
,则
|
R
|
=
|
∫
a
/
s
b
/
s
A
d
y
−
∫
−
∞
∞
A
0
d
y
|
{\displaystyle \left|R\right|=\left|\int _{a/s}^{b/s}A\,dy-\int _{-\infty }^{\infty }A_{0}\,dy\right|}
所以,只要我们求得上式的上限为何,即可得相对误差。注意,这里的推导还可以再优化,不过,由于此处我只想显示它会收敛到0,因此,不再细推。
由于
|
A
+
B
|
≤
|
A
|
+
|
B
|
{\displaystyle \left|A+B\right|\leq |A|+|B|}
,因此
|
R
|
<
|
∫
D
y
∞
A
0
d
y
|
⏟
(
a
1
)
+
|
∫
D
y
b
/
s
A
d
y
|
⏟
(
b
1
)
+
|
∫
−
D
y
D
y
(
A
−
A
0
)
d
y
|
⏟
(
c
)
+
|
∫
a
/
s
−
D
y
A
d
y
|
⏟
(
b
2
)
+
|
∫
−
∞
−
D
y
A
0
d
y
|
⏟
(
a
2
)
{\displaystyle |R|<\underbrace {\left|\int _{D_{y}}^{\infty }A_{0}dy\right|} _{(a_{1})}+\underbrace {\left|\int _{D_{y}}^{b/s}Ady\right|} _{(b_{1})}+\underbrace {\left|\int _{-D_{y}}^{D_{y}}\left(A-A_{0}\right)dy\right|} _{(c)}+\underbrace {\left|\int _{a/s}^{-D_{y}}Ady\right|} _{(b_{2})}+\underbrace {\left|\int _{-\infty }^{-D_{y}}A_{0}dy\right|} _{(a_{2})}}
其中
(
a
1
)
{\displaystyle (a_{1})}
与
(
a
2
)
{\displaystyle (a_{2})}
雷同,所以只算
(
a
1
)
{\displaystyle (a_{1})}
,经过
z
≡
π
y
2
{\displaystyle z\equiv \pi y^{2}}
的变换后,得到
(
a
1
)
=
|
1
2
π
∫
π
D
y
2
∞
e
−
z
z
−
1
/
2
d
z
|
<
e
−
π
D
y
2
2
π
D
y
{\displaystyle (a_{1})=\left|{\frac {1}{2{\sqrt {\pi }}}}\int _{\pi D_{y}^{2}}^{\infty }e^{-z}z^{-1/2}dz\right|<{\frac {e^{-\pi D_{y}^{2}}}{2\pi D_{y}}}}
也就是说,只要
D
y
{\displaystyle D_{y}}
取得够大,它就会趋近于0。
而
(
b
1
)
{\displaystyle (b_{1})}
与
(
b
2
)
{\displaystyle (b_{2})}
也雷同,所以也只算一个:
(
b
1
)
≤
|
∫
D
y
b
/
s
[
h
(
s
y
)
h
(
0
)
]
max
e
M
m
(
s
y
)
d
y
|
{\displaystyle (b_{1})\leq \left|\int _{D_{y}}^{b/s}\left[{\frac {h(sy)}{h(0)}}\right]_{\text{max}}e^{Mm(sy)}dy\right|}
其中
m
(
x
)
≥
1
M
ln
h
(
x
)
h
(
0
)
+
g
(
x
)
−
g
(
0
)
as
x
∈
[
s
D
y
,
b
]
{\displaystyle m(x)\geq {\frac {1}{M}}\ln {\frac {h(x)}{h(0)}}+g(x)-g(0)\,\,{\text{as}}\,\,x\in [sD_{y},b]}
并且
h
(
x
)
{\displaystyle h(x)}
在此范围内要与
h
(
0
)
{\displaystyle h(0)}
同号。
这里让我们选过
x
=
s
D
y
{\displaystyle x=sD_{y}}
的切线为
m
(
x
)
{\displaystyle m(x)}
吧!即
m
(
s
y
)
=
g
(
s
D
y
)
−
g
(
0
)
+
g
′
(
s
D
y
)
(
s
y
−
s
D
y
)
{\displaystyle m(sy)=g(sD_{y})-g(0)+g'(sD_{y})\left(sy-sD_{y}\right)}
,如图所示:
m
(
x
)
{\displaystyle m(x)}
以斜直线代表,为通过位于
x
=
s
D
y
{\displaystyle x=sD_{y}}
的切线
由图可以看出,当
s
{\displaystyle s}
或者
D
y
{\displaystyle D_{y}}
变小时,满足上列不等式的区间会变大,因此,为了涵盖整个区间,
D
y
{\displaystyle D_{y}}
有了上限。此外,
e
−
α
x
{\displaystyle e^{-\alpha x}}
的积分也相对容易,因此,很适合用来预估误差。
由泰勒展开我们得知,
M
[
g
(
s
D
y
)
−
g
(
0
)
]
=
M
[
g
″
(
0
)
2
s
2
D
y
2
+
g
‴
(
ξ
)
6
s
3
D
y
3
]
as
ξ
∈
[
0
,
s
D
y
]
=
−
π
D
y
2
+
(
2
π
)
3
/
2
g
‴
(
ξ
)
D
y
3
6
M
|
g
″
(
0
)
|
3
/
2
,
{\displaystyle {\begin{aligned}M\left[g(sD_{y})-g(0)\right]&=M\left[{\frac {g''(0)}{2}}s^{2}D_{y}^{2}+{\frac {g'''(\xi )}{6}}s^{3}D_{y}^{3}\right]\,\,{\text{as}}\,\,\xi \in [0,sD_{y}]\\&=-\pi D_{y}^{2}+{\frac {(2\pi )^{3/2}g'''(\xi )D_{y}^{3}}{6{\sqrt {M}}|g''(0)|^{3/2}}},\end{aligned}}}
且
M
s
g
′
(
s
D
y
)
=
M
s
(
g
″
(
0
)
s
D
y
+
g
‴
(
ζ
)
2
s
2
D
y
2
)
,
as
ζ
∈
[
0
,
s
D
y
]
=
−
2
π
D
y
+
2
M
(
π
|
g
″
(
0
)
|
)
3
/
2
g
‴
(
ζ
)
D
y
2
{\displaystyle {\begin{aligned}Msg'(sD_{y})&=Ms\left(g''(0)sD_{y}+{\frac {g'''(\zeta )}{2}}s^{2}D_{y}^{2}\right),\,\,{\text{as}}\,\,\zeta \in [0,sD_{y}]\\&=-2\pi D_{y}+{\sqrt {\frac {2}{M}}}\left({\frac {\pi }{|g''(0)|}}\right)^{3/2}g'''(\zeta )D_{y}^{2}\end{aligned}}}
然后将它们代回
(
b
1
)
{\displaystyle (b_{1})}
的计算里,不过,您可以看到,上述的两个余项皆与
M
{\displaystyle M}
的开根号成反比,为了式子的漂亮,让我省略它们吧!不省略只是看起来较丑陋而已,不过,那样子较严谨便是。
(
b
1
)
≤
|
[
h
(
s
y
)
h
(
0
)
]
max
e
−
π
D
y
2
∫
0
b
/
s
−
D
y
e
−
2
π
D
y
y
d
y
|
≤
|
[
h
(
s
y
)
h
(
0
)
]
max
e
−
π
D
y
2
1
2
π
D
y
|
.
{\displaystyle {\begin{aligned}(b_{1})&\leq \left|\left[{\frac {h(sy)}{h(0)}}\right]_{\text{max}}e^{-\pi D_{y}^{2}}\int _{0}^{b/s-D_{y}}e^{-2\pi D_{y}y}dy\right|\\&\leq \left|\left[{\frac {h(sy)}{h(0)}}\right]_{\text{max}}e^{-\pi D_{y}^{2}}{\frac {1}{2\pi D_{y}}}\right|.\end{aligned}}}
所以,原则上也是
D
y
{\displaystyle D_{y}}
越大则越趋近于0。只不过,在决定
m
(
x
)
{\displaystyle m(x)}
的过程,设下了
D
y
{\displaystyle D_{y}}
的上限。
至于靠近极大值的点的积分,我们一样可以利用泰勒展开来求,当
h
(
x
)
{\displaystyle h(x)}
在此点的一阶微分不为0时,
(
c
)
≤
∫
−
D
y
D
y
e
−
π
y
2
|
s
h
′
(
ξ
)
h
(
0
)
y
|
d
y
<
2
π
M
|
g
″
(
0
)
|
|
h
′
(
ξ
)
h
(
0
)
y
|
max
(
1
−
e
−
π
D
y
2
)
{\displaystyle {\begin{aligned}(c)&\leq \int _{-D_{y}}^{D_{y}}e^{-\pi y^{2}}\left|{\frac {sh'(\xi )}{h(0)}}y\right|\,dy\\&<{\sqrt {\frac {2}{\pi M|g''(0)|}}}\left|{\frac {h'(\xi )}{h(0)}}y\right|_{\text{max}}\left(1-e^{-\pi D_{y}^{2}}\right)\end{aligned}}}
会看到它原则上与
M
{\displaystyle M}
的开根号成反比,其实,当
h
(
x
)
{\displaystyle h(x)}
为常数时,积分出来的也会有如此的行为。
所以,概括的讲,靠近极大点的附近的积分原则上会随着
M
{\displaystyle {\sqrt {M}}}
的增大而变小,而其余的部分,只要
D
y
{\displaystyle D_{y}}
够大,也会趋近于0,只是这
D
y
{\displaystyle D_{y}}
是有上限的,取决于我们所找的上限函数
m
(
x
)
{\displaystyle m(x)}
是否在这区间内总是大于
g
(
x
)
−
g
(
0
)
{\displaystyle g(x)-g(0)}
;不过,一旦有一个满足条件的
m
(
x
)
{\displaystyle m(x)}
被找到,由于切线的起点为
x
=
s
D
y
{\displaystyle x=sD_{y}}
,因此,
D
y
{\displaystyle D_{y}}
可以取正比于
M
{\displaystyle {\sqrt {M}}}
即可,所以,
M
{\displaystyle M}
越大,
D
y
{\displaystyle D_{y}}
也可以设越大。
至于多维的情形,让我们令
x
{\displaystyle \mathbf {x} }
是一个
n
{\displaystyle n}
维向量,而
f
(
x
)
{\displaystyle f(\mathbf {x} )}
则是一个标量函数,则此拉普拉斯方法可以写成
∫
e
M
f
(
x
)
d
n
x
≈
(
2
π
M
)
n
/
2
|
−
H
(
f
)
(
x
0
)
|
−
1
/
2
e
M
f
(
x
0
)
as
M
→
∞
{\displaystyle \int e^{Mf(\mathbf {x} )}\,d^{n}x\approx \left({\frac {2\pi }{M}}\right)^{n/2}|-H(f)(\mathbf {x} _{0})|^{-1/2}e^{Mf(\mathbf {x} _{0})}{\text{ as }}M\to \infty \,}
其中的
H
(
f
)
(
x
0
)
{\displaystyle H(f)(\mathbf {x} _{0})}
是一个在
x
0
{\displaystyle \mathbf {x} _{0}}
取值的 海森矩阵 ,而
|
⋅
|
{\displaystyle |\cdot |}
则是指矩阵的 行列式 ;此外,与单变数的拉普拉斯方法类似,这里的 海森矩阵 必须为 负定矩阵 ,即该矩阵的所有本征值皆小于0,这样才会是极大值所在。.[ 2]
拉普拉斯方法的推广:最速下降法
此拉普拉斯方法可以被推广到 复分析 里头使用,搭配 留数定理 ,以找出一个过最速下降点的 contour (翻路径的话,会与path integral 相冲,所以,这里还是以英文原字称呼) 的 曲线积分 ,用来取代原有的复数空间的 contour积分。因为有时
f
(
z
)
{\displaystyle f(z)}
的一阶微分为0的点未必就在实数轴上,而就算真在实数轴上,也未必二阶微分在
x
{\displaystyle x}
方向上为小于0 的实数;此种情况下,就得使用最速下降法了。由于最速下降法中,已经利用另一条通过最速下降的鞍点来取代原有的 contour 积分,经过变数变换后就会变得有如拉普拉斯方法,因此,我们可以透过这新的 contour ,找到原本的积分的渐进近似解,而这将大大的简化整个计算。就好像原本的路径像是在蜿蜒的山路开车,而新的路径就像干脆绕过这座山开,反正目的只是到达目的地而已,留数定理已经帮我们把中间的差都算好了。请读 Erdelyi (2012)与 Arfken & Weber (2005) 的书里有关 steepest descents 的章节。
以下就是该方法在z 平面下的形式:
∫
a
b
e
M
f
(
z
)
d
z
≈
2
π
−
M
f
″
(
z
0
)
e
M
f
(
z
0
)
as
M
→
∞
.
{\displaystyle \int _{a}^{b}\!e^{Mf(z)}\,dz\approx {\sqrt {\frac {2\pi }{-Mf''(z_{0})}}}e^{Mf(z_{0})}{\text{ as }}M\to \infty .\,}
其中 z 0 就是新的路径通过的鞍点。
注意,开根号里的负号是用来指定最速下降的方向,千万别认为取
f
″
(
z
0
)
{\displaystyle f''(z_{0})}
的绝对值来取代这个负号,若然,那就大错特错了。
另外要注意的是,如果该被积函数是 在手动字词转换规则中检测到错误 ,就有必要将被新旧 contour 包到的极点所贡献的留数给加入,范例请参考 Okounkov 的文章 arXiv:math/0309074 的第三章。
更进一步一般化
最速下降法还可以更进一步的推广到所谓的 非线性稳定相位/最速下降法 (nonlinear stationary phase/steepest descent method)。
这方法主要用在解非线性偏微分方程,透过将非线性偏微分方程转换为求解柯西变换(Cauchy transform)的积分形式,就可以借由最速下降法的想法来得到非线性解的渐进解。
以 艾里方程 (线性)
y
″
=
x
y
{\displaystyle y''=xy}
为例,它可以写成积分形式如下:
y
j
(
x
)
=
1
2
π
i
∫
γ
j
e
1
3
t
3
−
x
t
d
t
{\displaystyle y_{j}(x)={\frac {1}{2\pi i}}\int _{\gamma _{j}}e^{{\frac {1}{3}}t^{3}-xt}dt}
由这条积分式,我们就可以借由最速下降法(若
γ
j
{\displaystyle \gamma _{j}}
指的是负实数轴,那么就回到此拉普拉斯方法了)来得到它的渐进解了。
然而,若方程式如 KdV方程
y
t
+
6
y
y
x
+
y
x
x
x
=
0
{\displaystyle y_{t}+6yy_{x}+y_{xxx}=0}
是个非线性偏微分方程,想要找到它的解相对应的一个复数 contour 积分的话,就没那么简单,在非线性稳定相位/最速下降法里所用到的概念主要是基于散射逆转换 (inverse scattering transform) 的处理方式,即先将原本的非线性偏微分方程变成 Lax 对 ,其中一个像是线性的 薛定谔方程式 ,其位能障为我们要找的
y
{\displaystyle y}
,本征值为
λ
{\displaystyle \lambda }
,波函数为
ϕ
λ
{\displaystyle \phi _{\lambda }}
(不过,它并非我们所要的
y
{\displaystyle y}
);因此可以解它的散射矩阵,若利用解析延拓将原本的波函数由实数
λ
{\displaystyle \lambda }
延拓到复数空间时,就可以得到黎曼希尔伯特问题(RHP)的形式。利用这个黎曼希尔伯特问题(RHP) ,我们可以解得
ϕ
{\displaystyle \phi }
的柯西变换的积分形式,再利用此线性薛定谔方程的特性,就可以反推出
y
{\displaystyle y}
的复数 contour 积分 形式了。
而 Lax 对 的另一个偏微分方程则是决定每个
ϕ
λ
{\displaystyle \phi _{\lambda }}
随时间变化的行为,由于
y
{\displaystyle y}
在
x
→
±
∞
{\displaystyle x\rightarrow \pm \infty }
时被要求为0 ,会发现整个偏微分方程会变得十分简单,并且只决定
c
(
λ
,
t
)
ϕ
λ
{\displaystyle c(\lambda ,t)\phi _{\lambda }}
里的
c
(
λ
,
t
)
{\displaystyle c(\lambda ,t)}
的值,不过,条件是
ϕ
λ
{\displaystyle \phi _{\lambda }}
必须是指由正负无限远入射或反射波的解。这样,我们所得到的这个只与时间与本征值有关的系数
c
(
λ
,
t
)
{\displaystyle c(\lambda ,t)}
就可以直接被应用在上述的黎曼希尔伯特问题(RHP)里的跃变矩阵里了。
接着就是非线性稳定相位/最速下降法所要做的工作,即找出 鞍点 来,在该点附近基于最速下降法的精神做近似。不过,这近似因着考虑到收敛性,需要将原本的 contour 变形,与将原本的黎曼希尔伯特问题(RHP)作转换,所以有再比原本的最速下降法多出一些步骤来。
这整个方法最早由 Deift 与 Zhou 在 1993 基于 Its 之前的工作所提出的,后续又有许多人加以改进,主要的应用则有 孤波 理论,可积模型等。
复数积分
以下的积分常用在 拉普拉斯变换#拉普拉斯逆变换 里,
1
2
π
i
∫
c
−
i
∞
c
+
i
∞
g
(
s
)
e
s
t
d
s
{\displaystyle {\frac {1}{2\pi i}}\int _{c-i\infty }^{c+i\infty }g(s)e^{st}\,ds}
假定我们想要得到该积分在
t
>>
1
{\displaystyle t>>1}
时的结果(若
t
{\displaystyle t}
为时间,通常就是在找经过够久时间后达稳定的结果),我们可以透过 解析延拓 的概念,先将这时间换成虚数,如 t = iu 并且一并做
s
=
c
+
i
x
{\displaystyle s=c+ix}
的变换,则我们可以将上式转换为如下的 拉普拉斯变换#双边拉普拉斯变换
1
2
π
∫
−
∞
∞
g
(
c
+
i
x
)
e
−
u
x
e
i
c
u
d
x
.
{\displaystyle {\frac {1}{2\pi }}\int _{-\infty }^{\infty }g(c+ix)e^{-ux}e^{icu}\,dx.}
这里就可以套用此拉普拉斯方法求渐进解,最后,再利用 u = t / i 把 t 换回来,就可以得到该拉普拉斯逆变换的渐进解了。
例子1:斯特灵公式
拉普拉斯方法可以用在推导 斯特灵公式 上;当 N 很大时,
N
!
≈
2
π
N
N
N
e
−
N
{\displaystyle N!\approx {\sqrt {2\pi N}}N^{N}e^{-N}\,}
证明:
由 Γ函数 的积分定义,我们可以得到
N
!
=
Γ
(
N
+
1
)
=
∫
0
∞
e
−
x
x
N
d
x
.
{\displaystyle N!=\Gamma (N+1)=\int _{0}^{\infty }e^{-x}x^{N}\,dx.}
接着让我们做变数变换,
x
=
N
z
{\displaystyle x=Nz\,}
因此
d
x
=
N
d
z
.
{\displaystyle dx=N\,dz.}
将这些代回 Γ函数 的积分定义里,我们可以得到
N
!
=
∫
0
∞
e
−
N
z
(
N
z
)
N
N
d
z
=
N
N
+
1
∫
0
∞
e
−
N
z
z
N
d
z
=
N
N
+
1
∫
0
∞
e
−
N
z
e
N
ln
z
d
z
=
N
N
+
1
∫
0
∞
e
N
(
ln
z
−
z
)
d
z
.
{\displaystyle {\begin{aligned}N!&=\int _{0}^{\infty }e^{-Nz}\left(Nz\right)^{N}N\,dz\\&=N^{N+1}\int _{0}^{\infty }e^{-Nz}z^{N}\,dz\\&=N^{N+1}\int _{0}^{\infty }e^{-Nz}e^{N\ln z}\,dz\\&=N^{N+1}\int _{0}^{\infty }e^{N(\ln z-z)}\,dz.\end{aligned}}}
经由此变数变换后,我们有了拉普拉斯方法所需要的
f
(
x
)
{\displaystyle f(x)}
为
f
(
z
)
=
ln
z
−
z
{\displaystyle f\left(z\right)=\ln {z}-z}
而它乃为二次可微函数,且
f
′
(
z
)
=
1
z
−
1
,
{\displaystyle f'(z)={\frac {1}{z}}-1,\,}
f
″
(
z
)
=
−
1
z
2
.
{\displaystyle f''(z)=-{\frac {1}{z^{2}}}.\,}
因此, ƒ (z ) 的极大值出现在 z 0 = 1 而且在该点的二次微分为
−
1
{\displaystyle -1}
。因此,我们得到
N
!
≈
N
N
+
1
2
π
N
e
−
N
=
2
π
N
N
N
e
−
N
.
{\displaystyle N!\approx N^{N+1}{\sqrt {\frac {2\pi }{N}}}e^{-N}={\sqrt {2\pi N}}N^{N}e^{-N}.\,}
关于概率推理的简介,请参考 http://doc.baidu.com/view/e9c1c086b9d528ea81c77989.html 。
而在 Azevedo-Filho & Shachter 1994 的文章里,则回顾了如何应用此拉普拉斯方法 (无论是单变数或者多变数) 如何应用在 概率推理 上,以加速得到系统的 后验矩 (数学) (posterior moment) , 贝氏参数 等,并举医学诊断上的应用为例。
相关维基百科文章
参考文献
Deift, P.; Zhou, X., A steepest descent method for oscillatory Riemann–Hilbert problems. Asymptotics for the MKdV equation, Ann. of Math. 137 (2), 1993, 137 (2): 295–368, doi:10.2307/2946540 .
Azevedo-Filho, A.; Shachter, R., Laplace's Method Approximations for Probabilistic Inference in Belief Networks with Continuous Variables, Mantaras, R.; Poole, D. (编), Uncertainty in Artificial Intelligence, San Francisco, CA: Morgan Kauffman, 1994, CiteSeerX : 10.1.1.91.2064 .
Erdelyi, A., Asymptotic Expansions , Dover Publications, 2012, ISBN 9780486155050 .
Arfken, G.B.; Weber, H.J., Mathematical Methods For Physicists International Student Edition , Elsevier Science, 2005, ISBN 9780080470696 .