输入参数化
当系统的随机输入是系统参数时,参数化过程将很简单,因为输入已经是参数形式。我们关心的问题是这些系统参数之间是否独立, 具体而言:
令
Y
=
(
Y
1
,
⋯
,
Y
n
)
,
n
>
1
Y=(Y_1,\cdots, Y_n), n>1
Y=(Y1,⋯,Yn),n>1 是系统参数, 其对应的分布函数为
F
Y
(
y
)
=
P
(
Y
≤
y
)
F_Y(y)=P(Y\leq y)
FY(y)=P(Y≤y), 我们希望寻找一组相互独立的随机变量
Z
=
(
Z
1
,
⋯
,
Z
n
)
∈
R
d
,
1
≤
d
≤
n
Z=(Z_1,\cdots,Z_n)\in\mathbb{R}^d,1\leq d\leq n
Z=(Z1,⋯,Zn)∈Rd,1≤d≤n, 使得
Y
=
T
(
Z
)
Y=T(Z)
Y=T(Z), 其中
T
T
T 是一个合适的变换.
我们用一个简单的例子来解释上述事情。 考虑一个带两个随机参数的常微分方程:
d
u
d
t
(
t
,
ω
)
=
−
α
(
ω
)
u
,
u
(
0
,
ω
)
=
β
(
ω
)
,
\frac{du}{dt}(t,\omega)=-\alpha(\omega)u,\quad u(0,\omega)=\beta(\omega),
dtdu(t,ω)=−α(ω)u,u(0,ω)=β(ω),
其中
α
,
β
\alpha,\beta
α,β 是随机的, 也就是说: 输入变量
Y
(
ω
)
=
(
α
(
ω
)
,
β
(
ω
)
)
∈
R
2
Y(\omega)=(\alpha(\omega),\beta(\omega))\in \mathbb{R}^2
Y(ω)=(α(ω),β(ω))∈R2.
如果 α , β \alpha,\beta α,β 是相互独立的, 则我们令 Z ( ω ) = Y ( ω ) Z(\omega)=Y(\omega) Z(ω)=Y(ω). 则解 u ( t , ω ) : [ 0 , T ] × Ω → R u(t,\omega):[0,T]\times \Omega\to \mathbb{R} u(t,ω):[0,T]×Ω→R 可以被表述为 u ( T , Z ) : [ 0 , T ] × R 2 → R , u(T,Z):[0,T]\times \mathbb{R}^2\to\mathbb{R}, u(T,Z):[0,T]×R2→R, 即方程有一个时间维度和两个随机维度。
如果
α
,
β
\alpha,\beta
α,β 不是相互独立的, 则存在函数
f
f
f 使得
f
(
α
,
β
)
=
0
f(\alpha,\beta)=0
f(α,β)=0. 则我们可以找到一个随机变量
Z
(
ω
)
Z(\omega)
Z(ω) 使得
α
(
ω
)
=
a
(
Z
(
ω
)
)
,
β
(
ω
)
=
b
(
Z
(
ω
)
)
,
f
(
a
,
b
)
=
0.
\alpha(\omega)=a(Z(\omega)),\quad \beta(\omega)=b(Z(\omega)), \quad f(a,b)=0.
α(ω)=a(Z(ω)),β(ω)=b(Z(ω)),f(a,b)=0.
或者, 等价地,
α
,
β
\alpha,\beta
α,β 相关可以推出存在函数
g
g
g 使得
β
=
g
(
α
)
.
\beta=g(\alpha).
β=g(α). 接着我们令
Z
(
ω
)
=
α
(
ω
)
,
β
(
ω
)
=
g
(
Z
(
ω
)
)
.
Z(\omega)=\alpha(\omega),\beta(\omega)=g(Z(\omega)).
Z(ω)=α(ω),β(ω)=g(Z(ω)). 此时, 问题的解便变为
u
(
t
,
Z
)
:
[
0
,
T
]
×
R
→
R
,
u(t,Z):[0,T]\times \mathbb{R}\to \mathbb{R},
u(t,Z):[0,T]×R→R,
即只有一个随机维度。
实际问题中, 系统中有许多随机参数, 找到这些随机参数依赖的精确函数形式(即
f
,
g
f,g
f,g) 是很困难的。我们能知道的信息仅仅是他们对应的(联合)概率分布函数, 于是我们的目标便是使用分布函数来找到变换
T
T
T.
高斯参数
对于高斯参数, 这件事情是比较容易的。 假设
Y
=
(
Y
1
,
⋯
,
Y
n
)
∼
N
(
0
,
C
)
,
Y=(Y_1,\cdots,Y_n)\sim \mathcal{N}(0,\bm{C}),
Y=(Y1,⋯,Yn)∼N(0,C), 其中
C
∈
R
n
×
n
\bm{C}\in\mathbb{R}^{n\times n}
C∈Rn×n 是协方差矩阵。 令
Z
∼
N
(
0
,
I
)
,
I
∈
R
n
×
n
Z\sim \mathcal{N}(0,\bm{I}),\bm{I}\in \mathbb{R}^{n\times n}
Z∼N(0,I),I∈Rn×n.
A
\bm{A}
A 是一个
n
×
n
n\times n
n×n 的矩阵, 由高斯变换的性质我们知
A
Z
∼
N
(
0
,
A
A
T
)
\bm{A}Z\sim \mathcal{N}(0, \bm{AA^T})
AZ∼N(0,AAT). 因此, 存在矩阵
A
\bm{A}
A 使得
A
A
T
=
C
,
\bm{AA^T}=\bm{C},
AAT=C, 则
Y
=
A
Z
∼
N
(
0
,
C
)
.
Y=\bm{A}Z\sim\mathcal{N}(0,\bm{C}).
Y=AZ∼N(0,C). 实际计算时, 我们可以采用
C
h
o
l
e
s
k
y
Cholesky
Cholesky 分解来计算矩阵
A
\bm{A}
A. 具体有
a
i
1
=
c
i
1
/
c
11
,
1
≤
i
≤
n
,
a
i
i
=
c
i
i
−
∑
k
=
1
i
−
1
a
i
k
2
,
1
<
i
≤
n
,
a
i
j
=
(
c
i
j
−
∑
k
=
1
j
−
1
a
i
k
a
j
k
)
/
a
j
j
,
1
<
j
<
i
≤
n
,
a
i
j
=
0
,
i
<
j
<
n
.
\begin{array}{ll} a_{i1}=c_{i1}/\sqrt{c_{11}}, &1\leq i\leq n,\\ a_{ii}=\sqrt{c_{ii}-\sum_{k=1}^{i-1}a^2_{ik}},& 1<i\leq n,\\ a_{ij} = (c_{ij}-\sum_{k=1}^{j-1}a_{ik}a_{jk})/a_{jj}, &1<j<i\leq n,\\ a_{ij}=0, &i<j<n. \end{array}
ai1=ci1/c11,aii=cii−∑k=1i−1aik2,aij=(cij−∑k=1j−1aikajk)/ajj,aij=0,1≤i≤n,1<i≤n,1<j<i≤n,i<j<n.
非高斯参数
当系统参数是非高斯分布的时候, 直接参数化变变得很困难了。 Rosenblatt 给出了如下结果来理论上解决这一问题:
令
Y
=
(
Y
1
,
⋯
,
Y
n
)
Y=(Y_1,\cdots,Y_n)
Y=(Y1,⋯,Yn) 是一个随机向量, 其分布函数为
F
Y
(
y
)
=
P
(
Y
≤
y
)
,
F_Y(y)=P(Y\leq y),
FY(y)=P(Y≤y), 令
Z
=
(
z
1
,
⋯
,
z
n
)
=
T
y
=
T
(
y
1
,
⋯
,
y
n
)
Z=(z_1,\cdots,z_n)=Ty=T(y_1,\cdots,y_n)
Z=(z1,⋯,zn)=Ty=T(y1,⋯,yn) 为
z
1
=
P
(
Y
1
≤
y
1
)
=
F
1
(
y
1
)
,
z
2
=
P
(
Y
2
≤
y
2
∣
Y
1
=
y
1
)
=
F
2
(
y
2
∣
y
1
)
,
⋯
,
z
n
=
P
(
Y
n
≤
y
n
∣
Y
n
−
1
=
y
n
−
1
,
⋯
,
Y
1
=
y
1
)
=
F
n
(
y
n
∣
,
y
n
−
1
,
⋯
,
y
1
)
.
\begin{array}{l} z_1=P(Y_1\leq y_1)=F_1(y_1),\\ z_2 = P(Y_2\leq y_2\mid Y_1=y_1)=F_2(y_2\mid y_1),\\ \cdots,\\ z_n = P(Y_n\leq y_n\mid Y_{n-1}=y_{n-1},\cdots,Y_1=y_1)=F_n(y_n\mid,y_{n-1},\cdots,y_1). \end{array}
z1=P(Y1≤y1)=F1(y1),z2=P(Y2≤y2∣Y1=y1)=F2(y2∣y1),⋯,zn=P(Yn≤yn∣Yn−1=yn−1,⋯,Y1=y1)=Fn(yn∣,yn−1,⋯,y1).
则简单计算有
P
(
Z
i
≤
z
i
;
i
=
1
,
2
,
⋯
,
n
)
=
∫
{
Z
∣
Z
i
≤
z
i
}
⋯
∫
d
y
n
F
n
(
y
n
∣
y
n
−
1
,
⋯
,
y
1
)
⋯
d
y
1
F
1
(
y
1
)
=
∫
0
z
n
⋯
∫
0
z
1
d
z
1
⋯
d
z
n
=
∏
i
=
1
n
z
i
,
\begin{array}{rl} P(Z_i&\leq z_i;i=1,2,\cdots,n)\\ &=\int_{\{Z\mid Z_i\leq z_i\}}\cdots \int d{y_n} F_n(y_n\mid y_{n-1},\cdots,y_1)\cdots d_{y_1} F_1(y_1)\\ &= \int _0^{z_n} \cdots \int_{0}^{z_1}dz_1\cdots dz_n =\prod_{i=1}^{n}z_i, \end{array}
P(Zi≤zi;i=1,2,⋯,n)=∫{Z∣Zi≤zi}⋯∫dynFn(yn∣yn−1,⋯,y1)⋯dy1F1(y1)=∫0zn⋯∫0z1dz1⋯dzn=∏i=1nzi,
即
Z
=
(
Z
1
,
⋯
,
Z
n
)
∼
[
0
,
1
]
n
Z=(Z_1,\cdots,Z_n)\sim [0,1]^n
Z=(Z1,⋯,Zn)∼[0,1]n 独立同分布.
尽管数学上此定理简单有力, 但在实际中基本上不可以使用,原因在于它依赖随机变量的条件概率分布,实际中这样的信息很少完全知道。况且,即使我们知道一些条件概率分布,也没有精确地表示形式。 故实际中我们会采用 Rosenblatt 的一些数值形式, 但这不属于我们讨论的范围。
随机过程与降维
在许多情况下, 随机输入往往是随机过程。 例如, 输入可以是随时间变化的随机强迫项, 此时它为时间上的随机过程; 或者是不确定的材料性质(导电率), 此时它为空间上的随机过程。 参数化的问题可以被描述为下述形式:
令
(
Y
t
,
t
∈
T
)
(Y_t, t\in T)
(Yt,t∈T) 是一个随机过程, 求变换
R
R
R 使得
Y
t
=
R
(
Z
)
Y_t = R(Z)
Yt=R(Z), 其中
Z
=
(
Z
1
,
⋯
,
Z
d
)
,
d
≥
1
Z=(Z_1,\cdots,Z_d), d\geq 1
Z=(Z1,⋯,Zd),d≥1 是相互独立的。
注意到, 指标集
T
T
T 可以是时间区域或者空间区域,因而其通常是一个无限维的对象。 因为我们要求
d
d
d 是有限的, 故上述变换不是精确的, 即
Y
t
≈
R
(
Z
)
Y_t\approx R(Z)
Yt≈R(Z).
一种直接的方式便是考虑
Y
t
Y_t
Yt 的有限维版本, 这便要求将
T
T
T 离散为有限维的指标进而考虑
(
Y
t
1
,
Y
t
2
,
⋯
,
Y
t
n
)
,
t
1
,
t
2
,
⋯
,
t
n
∈
T
.
(Y_{t_1},Y_{t_2},\cdots,Y_{t_n}),\quad t_1,t_2,\cdots, t_n\in T.
(Yt1,Yt2,⋯,Ytn),t1,t2,⋯,tn∈T.
Y
t
Y_t
Yt 的离散是一个逼近问题,显然,离散越密, 逼近效果越好。但是, 这会导致我们问题的维数变大,计算代价变大。 此处我们讨论三种方式:
Karhunen-Loeve Expansion
KL 展开是表达随机过程中最常用模型降维方法。具体而言, 对于任意一个随机函数
κ
(
x
,
ω
)
\kappa(x,\omega)
κ(x,ω), 我们定义其均值函数, 即数学期望:
κ
ˉ
(
x
)
=
E
[
κ
(
x
,
⋅
)
]
=
∫
Ω
κ
(
x
,
ω
)
d
P
,
x
∈
D
.
\bar{\kappa}(x)=\mathbb{E}[\kappa(x,\cdot)]=\int_\Omega \kappa(x,\omega)dP,\quad x\in D.
κˉ(x)=E[κ(x,⋅)]=∫Ωκ(x,ω)dP,x∈D.
类似地, 我们可以定义其 $\tau
阶
矩
(
阶矩(
阶矩(\tau\in \mathbb
N$)
κ
ˉ
τ
(
x
)
=
E
[
κ
τ
(
x
,
⋅
)
]
=
∫
Ω
κ
τ
(
x
,
ω
)
d
P
,
x
∈
D
.
\bar{\kappa}_\tau (x)=\mathbb{E}[\kappa^\tau (x,\cdot)]=\int_\Omega \kappa^\tau (x,\omega)dP,\quad x\in D.
κˉτ(x)=E[κτ(x,⋅)]=∫Ωκτ(x,ω)dP,x∈D.
随机函数
κ
(
x
,
ω
)
\kappa(x,\omega)
κ(x,ω) 的协方差函数定义为
C
o
v
κ
(
x
,
y
)
=
E
[
(
κ
(
x
,
⋅
)
−
κ
ˉ
(
x
)
)
(
κ
(
y
,
⋅
)
−
κ
ˉ
(
y
)
)
]
,
x
,
y
∈
D
.
Cov_\kappa (x,y) = \mathbb{E}[(\kappa(x,\cdot)-\bar{\kappa}(x))(\kappa(y,\cdot)-\bar{\kappa}(y))],\quad x,y\in D.
Covκ(x,y)=E[(κ(x,⋅)−κˉ(x))(κ(y,⋅)−κˉ(y))],x,y∈D.
根据定义,
C
o
v
κ
(
x
,
y
)
:
D
×
D
‾
→
R
Cov_\kappa(x,y):\overline{D\times D}\to \mathbb{R}
Covκ(x,y):D×D→R 是对称正定的协方差函数。 我们引入自共轭紧算子
T
k
:
L
2
(
D
)
→
L
2
(
D
)
:
T_k:L^2(D)\to L^2(D):
Tk:L2(D)→L2(D):
T
k
v
(
⋅
)
=
∫
D
C
o
v
κ
(
x
,
⋅
)
v
(
x
)
d
x
,
∀
v
∈
L
2
(
D
)
.
T_kv(\cdot)=\int_D Cov_\kappa (x,\cdot)v(x)dx,\quad \forall v\in L^2(D).
Tkv(⋅)=∫DCovκ(x,⋅)v(x)dx,∀v∈L2(D).
令
{
λ
i
}
i
=
1
∞
\{\lambda_i\}_{i=1}^{\infty}
{λi}i=1∞ 和
{
ϕ
i
}
i
=
1
∞
\{\phi_i\}_{i=1}^\infty
{ϕi}i=1∞ 分别为算子
T
κ
T_\kappa
Tκ 的非负特征值序列和归一化正交特征函数, 即
T
κ
ϕ
i
=
λ
i
ϕ
i
,
λ
1
≥
λ
2
≥
⋯
>
0
,
⟨
ϕ
i
,
ϕ
j
⟩
=
δ
i
j
,
i
,
j
∈
N
+
.
T_\kappa \phi_i = \lambda_i \phi_i, \lambda_1\geq \lambda_2\geq \cdots >0,\quad \langle \phi_i,\phi_j\rangle =\delta_{ij},\quad i,j\in\mathbb{N}^+.
Tκϕi=λiϕi,λ1≥λ2≥⋯>0,⟨ϕi,ϕj⟩=δij,i,j∈N+.
进一步, 我们定义互不相关、且均值为0方差为1(i.i.d)的随机变量
Z
i
(
ω
)
≔
1
λ
i
∫
D
(
κ
(
x
,
ω
)
−
κ
ˉ
(
x
)
)
ϕ
i
(
x
)
d
x
,
i
=
1
,
2
,
⋯
,
E
[
Z
i
]
=
0
,
E
[
Z
i
Z
j
]
=
δ
i
j
,
i
,
j
∈
N
+
.
\begin{array}{l} Z_i(\omega)\coloneqq \frac{1}{\sqrt{\lambda_i}}\int_D (\kappa(x,\omega)-\bar{\kappa}(x))\phi_i(x)dx,\quad i=1,2,\cdots,\\ \mathbb{E}[Z_i]=0,\quad \mathbb{E}[Z_iZ_j]=\delta_{ij},\quad i,j\in \mathbb{N}^+. \end{array}
Zi(ω):=λi1∫D(κ(x,ω)−κˉ(x))ϕi(x)dx,i=1,2,⋯,E[Zi]=0,E[ZiZj]=δij,i,j∈N+.
那么,成立如下 Karhunen-Loeve 展开:
κ
(
x
,
ω
)
=
κ
ˉ
(
x
)
+
∑
i
=
1
∞
λ
i
Z
i
(
ω
)
ϕ
i
(
x
)
.
\kappa(x,\omega)=\bar{\kappa}(x)+\sum_{i=1}^{\infty}\sqrt{\lambda_i}Z_i(\omega)\phi_i(x).
κ(x,ω)=κˉ(x)+i=1∑∞λiZi(ω)ϕi(x).
对上述的Karhunen-Loeve 展开作有限项截断, 我们便可以得到关于随机函数
κ
(
x
,
ω
)
\kappa(x,\omega)
κ(x,ω) 的有限维逼近
κ
(
x
,
ω
)
≈
κ
d
(
x
,
ω
)
=
κ
ˉ
(
x
)
+
∑
i
=
1
d
λ
i
Z
i
(
ω
)
ϕ
i
(
x
)
.
\kappa(x,\omega)\approx \kappa_d(x,\omega)=\bar{\kappa}(x)+\sum_{i=1}^d\sqrt{\lambda_i}Z_i(\omega)\phi_i(x).
κ(x,ω)≈κd(x,ω)=κˉ(x)+i=1∑dλiZi(ω)ϕi(x).
并且, 由 Mercer 定理, 我们有
lim
d
→
∞
sup
x
∈
D
E
[
(
κ
(
x
,
⋅
)
−
κ
d
(
x
,
⋅
)
)
2
]
=
lim
d
→
∞
sup
x
∈
D
∑
i
=
d
+
1
∞
λ
i
ϕ
i
2
(
x
)
=
0
,
\lim\limits_{d\to \infty}\sup \limits_{x\in D}\mathbb{E}[(\kappa(x,\cdot)-\kappa_d(x,\cdot))^2]=\lim\limits_{d\to \infty}\sup\limits_{x\in D}\sum_{i=d+1}^{\infty}\lambda_i\phi_i^2(x)=0,
d→∞limx∈DsupE[(κ(x,⋅)−κd(x,⋅))2]=d→∞limx∈Dsupi=d+1∑∞λiϕi2(x)=0,
并且,在均方意义下, 有限阶截断的 Karhunen-Loeve 展开具有最优误差估计。如果你想较为严格全面的了解KL展开的具体细节, 请参考文献 A brief note on Karhunen-Loeve expansion。
高斯过程与非高斯过程
利用 KL 展开, 我们可以用有限个随机变量来逼近一个随机场, 并且 { Z i ( ω ) } \{Z_i(\omega)\} {Zi(ω)} 是不相关的。 对于高斯过程来说, 我们可以证明 Z i ( ω ) Z_i(\omega) Zi(ω) 也服从高斯分布, 进而我们便知 { Z i ( ω ) } \{Z_i(\omega)\} {Zi(ω)} 是相互独立的。 这是非常友好的。 对于非高斯过程而言, { Z i ( ω ) } \{Z_i(\omega)\} {Zi(ω)} 不是相互独立的, 因此 KL 展开不能提供独立参数化。 但是实际中, 我们仍然采用 KL 展开, 并且将 { Z i ( ω ) } \{Z_i(\omega)\} {Zi(ω)} 看作是相互独立的, 这当然时是不严格的。 此问题到目前仍然是个公开问题,有待进一步的讨论解决。
随机系统的参数化
考虑如下 PDE 系统:
{
u
t
(
x
,
t
,
ω
)
=
L
(
u
)
,
D
×
(
0
,
T
]
×
Ω
,
B
(
u
)
=
0
,
∂
D
×
[
0
,
T
]
×
Ω
,
u
=
u
0
,
D
×
{
t
=
0
}
×
Ω
,
\left\{\begin{array}{ll} u_t(x,t,\omega)=\mathcal{L}(u),& D\times (0,T]\times \Omega,\\ \mathcal{B}(u)=0,& \partial D\times [0,T]\times \Omega,\\ u=u_0,& D\times \{t=0\}\times \Omega, \end{array}\right.
⎩⎨⎧ut(x,t,ω)=L(u),B(u)=0,u=u0,D×(0,T]×Ω,∂D×[0,T]×Ω,D×{t=0}×Ω,
其中
L
\mathcal{L}
L 是(非线性)微分算子,
B
\mathcal{B}
B 是边界算子,
u
0
u_0
u0 是初始条件,
ω
∈
Ω
\omega\in \Omega
ω∈Ω 定义了系统的随机输入。 因此解便是一个随机量
u
(
x
,
t
,
ω
)
:
D
ˉ
×
[
0
,
T
]
×
Ω
→
R
n
u
,
u(x,t,\omega):\bar{D}\times [0,T]\times \Omega\to \mathbb{R}^{n_u},
u(x,t,ω):Dˉ×[0,T]×Ω→Rnu,
n
u
≥
1
n_u\geq 1
nu≥1 是
u
u
u 的维数。 由之前的知识我们知, 参数化后的系统为
{
u
t
(
x
,
t
,
Z
)
=
L
(
u
)
,
D
×
(
0
,
T
]
×
R
d
,
B
(
u
)
=
0
,
∂
D
×
[
0
,
T
]
×
R
d
,
u
=
u
0
,
D
×
{
t
=
0
}
×
R
d
,
\left\{\begin{array}{ll} u_t(x,t,Z)=\mathcal{L}(u),& D\times (0,T]\times \mathbb{R}^d,\\ \mathcal{B}(u)=0,& \partial D\times [0,T]\times \mathbb{R}^d,\\ u=u_0,& D\times \{t=0\}\times \mathbb{R}^d, \end{array}\right.
⎩⎨⎧ut(x,t,Z)=L(u),B(u)=0,u=u0,D×(0,T]×Rd,∂D×[0,T]×Rd,D×{t=0}×Rd,
其中
Z
=
(
Z
1
,
⋯
,
Z
d
)
Z=(Z_1,\cdots,Z_d)
Z=(Z1,⋯,Zd) 是一组独立的随机变量,用来参数化随机输入。我们举一个具体的一维随机椭圆问题:
{
∇
(
κ
(
x
,
ω
)
∇
u
)
=
f
(
x
,
ω
)
,
x
∈
(
−
1
,
1
)
,
u
(
−
1
,
ω
)
=
u
l
(
ω
)
,
u
(
1
,
ω
)
=
u
r
(
ω
)
.
\left\{\begin{array}{l} \nabla(\kappa(x,\omega)\nabla u )=f(x,\omega),\quad x\in(-1,1),\\ u(-1,\omega)=u_l(\omega),\quad u(1,\omega)=u_r(\omega). \end{array}\right.
{∇(κ(x,ω)∇u)=f(x,ω),x∈(−1,1),u(−1,ω)=ul(ω),u(1,ω)=ur(ω).
这里
κ
,
f
\kappa,f
κ,f 是随机场,
u
l
,
u
r
u_l,u_r
ul,ur 是随机变量。 假定随机场 $\kappa,f $ 有 KL 展开:
κ
(
x
,
ω
)
≈
κ
~
(
x
,
Z
κ
)
=
κ
ˉ
(
x
)
+
∑
i
=
1
d
κ
λ
i
κ
Z
i
κ
(
ω
)
ϕ
i
κ
(
x
)
,
f
(
x
,
ω
)
≈
f
~
(
x
,
Z
f
)
=
f
ˉ
(
x
)
+
∑
i
=
1
d
f
λ
i
f
Z
i
f
(
ω
)
ϕ
i
f
(
x
)
.
\begin{array}{l} \kappa(x,\omega)\approx \tilde{\kappa}(x,Z^\kappa)=\bar{\kappa}(x) + \sum_{i=1}^{d_{\kappa} }\sqrt{\lambda^{\kappa}_i}Z^{\kappa}_i(\omega)\phi^{\kappa}_i(x),\\ f(x,\omega)\approx \tilde{f}(x,Z^f)=\bar{f}(x)+\sum_{i=1}^{d_f} \sqrt{\lambda^f_i}Z^f_i(\omega)\phi^f_i(x). \end{array}
κ(x,ω)≈κ~(x,Zκ)=κˉ(x)+∑i=1dκλiκZiκ(ω)ϕiκ(x),f(x,ω)≈f~(x,Zf)=fˉ(x)+∑i=1dfλifZif(ω)ϕif(x).
我们假定
{
Z
f
}
,
{
Z
κ
}
,
u
l
,
u
r
\{Z^f\},\{Z^\kappa\}, u_l,u_r
{Zf},{Zκ},ul,ur 相互独立, 令
Z
=
(
Z
1
,
⋯
,
Z
d
)
=
Z
1
κ
,
⋯
,
Z
d
κ
κ
,
Z
1
f
,
⋯
,
Z
d
f
f
,
u
l
,
u
r
)
,
Z=(Z_1,\cdots,Z_d)= Z_1^\kappa,\cdots,Z^\kappa_{d_\kappa}, Z^f_1,\cdots,Z^f_{d_f},u_l,u_r),
Z=(Z1,⋯,Zd)=Z1κ,⋯,Zdκκ,Z1f,⋯,Zdff,ul,ur), 其中
d
=
d
κ
+
d
f
+
2
,
d=d^\kappa+d^f+2,
d=dκ+df+2, 则之前的问题可以表述为
{
∇
(
κ
~
(
x
,
Z
)
∇
u
)
=
f
~
(
x
,
Z
)
,
x
∈
(
−
1
,
1
)
,
u
(
−
1
,
Z
)
=
Z
d
−
1
,
u
(
1
,
Z
)
=
Z
d
.
\left\{ \begin{array}{l} \nabla(\tilde{\kappa}(x,Z)\nabla u)=\tilde{f}(x,Z),\quad x\in(-1,1),\\ u(-1,Z) = Z_{d-1}, u(1,Z)=Z_{d}. \end{array} \right.
{∇(κ~(x,Z)∇u)=f~(x,Z),x∈(−1,1),u(−1,Z)=Zd−1,u(1,Z)=Zd.
方程的解为
u
(
x
,
Z
)
:
[
−
1
,
1
]
×
R
d
→
R
.
u(x,Z):[-1,1]\times \mathbb{R}^d\to \mathbb{R}.
u(x,Z):[−1,1]×Rd→R.
几种传统的数值方法
我们以之前的常微分方程为例,即
d
u
d
t
(
t
,
ω
)
=
−
α
(
ω
)
u
,
u
(
0
,
ω
)
=
β
(
ω
)
.
\frac{du}{dt}(t,\omega)=-\alpha(\omega)u,\quad u(0,\omega)=\beta(\omega).
dtdu(t,ω)=−α(ω)u,u(0,ω)=β(ω).
上述方程的精确解为
u
(
t
,
ω
)
=
β
(
ω
)
e
−
α
(
ω
)
t
.
u(t,\omega)=\beta(\omega)e^{-\alpha(\omega)t}.
u(t,ω)=β(ω)e−α(ω)t.
如果
α
(
ω
)
,
β
(
ω
)
\alpha(\omega),\beta(\omega)
α(ω),β(ω) 是相互独立的, 则有
E
[
u
(
t
,
ω
)
]
=
E
[
β
]
E
[
e
−
α
t
]
.
\mathbb{E}[u(t,\omega)]=\mathbb{E}[\beta]\mathbb{E}[e^{-\alpha t}].
E[u(t,ω)]=E[β]E[e−αt].
Monte Carlo Sampling
Monte Carlo 方法. Monte Carlo 方法及其改进方法是简单常用的基于样本的方法. 在 Monte Carlo 方法中, 我们利用概率分布随机地产生一些样本, 对每一个样本而言, 所要解决的问题变成了一 个确定的问题. 通过求解这些确定问题, 可以得到精确解的一些统计量信息, 如均值或者方差. Monte Carlo 方法实施简单, 可以利用现成的程序, 利于并行计算. 但众所周知, Monte Carlo 方法的收敛率非 常低: 均值收敛率是 1 / K 1/\sqrt{K} 1/K ( K K K 是所使用样本的数量), 这意味着我们需要使用大量的样本才能得到较精确的数值结果. 如果所对应的确定问题求解困难, 这将是一个很大的挑战. 具体而言:
- 产生独立同分布的随机实现 { Z i = ( α i , β i ) } ; \{Z_i=(\alpha_i,\beta_i)\}; {Zi=(αi,βi)};
- 求解方程 u i ( t ) = u ( t , Z i ) ; u^i(t)=u(t,Z^i); ui(t)=u(t,Zi);
- 计算统计量信息 u ˉ ( t ) = 1 M ∑ i = 1 M u ( t , Z i ) ≈ E [ u ] . \bar{u}(t)=\frac{1}{M}\sum_{i=1}^M u(t,Z^i)\approx \mathbb{E}[u]. uˉ(t)=M1∑i=1Mu(t,Zi)≈E[u].
Moment Equation Approach
令
μ
(
t
)
=
E
[
u
]
,
\mu(t)=\mathbb{E}[u],
μ(t)=E[u], 在方程两边对
ω
\omega
ω 积分有,
d
μ
d
t
(
t
)
=
−
E
[
α
u
]
,
μ
(
0
)
=
E
[
β
]
.
\frac{d\mu}{dt}(t)=-\mathbb{E}[\alpha u],\quad \mu(0)=\mathbb{E}[\beta].
dtdμ(t)=−E[αu],μ(0)=E[β].
上述方程需要知道
E
[
α
u
]
\mathbb{E}[\alpha u]
E[αu] 的知识, 故我们两边乘
α
\alpha
α 再次积分有
d
d
t
E
[
α
u
]
(
t
)
=
−
E
[
α
2
u
]
,
E
[
α
u
]
(
0
)
=
E
[
α
β
]
,
\frac{d}{dt}\mathbb{E}[\alpha u](t)=-\mathbb{E}[\alpha^2 u],\quad \mathbb{E}[\alpha u](0)=\mathbb{E}[\alpha\beta],
dtdE[αu](t)=−E[α2u],E[αu](0)=E[αβ],
这便出现了一个新的量
E
[
α
2
u
]
\mathbb{E}[\alpha^2 u]
E[α2u], 我们对原始系统乘
α
2
\alpha^2
α2 再积分有
d
d
t
E
[
α
2
u
]
=
−
E
[
α
3
u
]
,
E
[
α
2
u
]
(
0
)
=
E
[
α
2
β
]
.
\frac{d}{dt}\mathbb{E}[\alpha^2u]=-\mathbb{E}[\alpha^3u],\quad \mathbb{E}[\alpha^2u](0)=\mathbb{E}[\alpha^2\beta].
dtdE[α2u]=−E[α3u],E[α2u](0)=E[α2β].
我们定义
μ
~
k
(
t
)
=
E
[
α
k
u
]
,
\tilde{\mu}_k(t)=\mathbb{E}[\alpha^k u],
μ~k(t)=E[αku], 接着我们有
d
d
t
μ
~
k
(
t
)
=
−
μ
~
k
+
1
(
t
)
,
μ
~
k
(
0
)
=
E
[
α
k
β
]
.
\frac{d}{dt}\tilde{\mu}_k(t)=-\tilde{\mu}_{k+1}(t), \quad \tilde{\mu}_k(0)=\mathbb{E}[\alpha^k \beta].
dtdμ~k(t)=−μ~k+1(t),μ~k(0)=E[αkβ].
即原始的系统是不封闭的,北京大学的李若老师在这方面做了很多杰出的工作。最简单的做法是固定一个
k
k
k, 假定高阶矩可以被低阶矩表示,即
μ
~
k
+
1
=
g
(
μ
~
0
,
⋯
,
μ
~
k
)
.
\tilde{\mu}_{k+1}=g(\tilde{\mu}_0,\cdots,\tilde{\mu}_k).
μ~k+1=g(μ~0,⋯,μ~k).
Perturbation method
当
ϵ
=
O
(
α
(
ω
)
)
∼
σ
α
<
<
1
,
\epsilon=O(\alpha(\omega))\sim \sigma_\alpha << 1,
ϵ=O(α(ω))∼σα<<1, 我们可以考虑解的Taylor展开:
u
(
t
,
ω
)
=
u
0
(
t
)
+
α
(
ω
)
u
1
(
t
)
+
α
2
(
ω
)
u
2
(
t
)
+
⋯
,
u(t,\omega)=u_0(t)+\alpha(\omega)u_1(t)+\alpha^2(\omega)u_2(t)+\cdots,
u(t,ω)=u0(t)+α(ω)u1(t)+α2(ω)u2(t)+⋯,
代入原方程, 我们有
d
u
0
d
t
+
α
d
u
1
d
t
+
α
2
d
u
2
d
t
+
⋯
=
−
α
(
u
0
+
α
u
1
+
α
2
u
2
+
⋯
)
.
\frac{du_0}{dt}+\alpha \frac{du_1}{dt}+\alpha^2\frac{du_2}{dt}+\cdots=-\alpha(u_0+\alpha u_1+\alpha^2 u_2+\cdots).
dtdu0+αdtdu1+α2dtdu2+⋯=−α(u0+αu1+α2u2+⋯).
因为
O
(
α
k
)
=
ϵ
k
O(\alpha^k)=\epsilon^k
O(αk)=ϵk, 于是我们有
O
(
1
)
:
d
u
0
d
t
=
0
,
O
(
ϵ
)
:
d
u
1
d
t
=
−
u
0
,
O
(
ϵ
2
)
:
d
u
2
d
t
=
−
u
1
,
⋯
⋯
O
(
ϵ
k
)
:
d
u
k
d
t
=
−
u
k
−
1
,
\begin{array}{ll} O(1): & \frac{du_0}{dt}=0,\\ O(\epsilon): & \frac{du_1}{dt}=-u_0,\\ O(\epsilon^2): & \frac{du_2}{dt}=-u_1,\\ \cdots&\cdots\\ O(\epsilon^k): & \frac{du_k}{dt}=-u_{k-1}, \end{array}
O(1):O(ϵ):O(ϵ2):⋯O(ϵk):dtdu0=0,dtdu1=−u0,dtdu2=−u1,⋯dtduk=−uk−1,
且
u
0
(
0
)
=
β
,
u
k
(
0
)
=
0
,
k
>
1.
u_0(0)=\beta,\quad u_k(0)=0,\quad k>1.
u0(0)=β,uk(0)=0,k>1.
简单计算有
u
0
(
t
)
=
β
,
u
1
(
t
)
=
−
β
t
,
⋯
,
u
k
=
β
(
−
1
)
k
t
k
k
!
.
u_0(t)=\beta, u_1(t)=-\beta t,\cdots, u_k=\beta (-1)^k \frac{t^k}{k!}.
u0(t)=β,u1(t)=−βt,⋯,uk=β(−1)kk!tk.
故
u
(
t
,
ω
)
=
∑
k
=
0
∞
β
−
t
k
k
!
α
k
(
ω
)
=
β
e
x
p
(
−
α
t
)
.
u(t,\omega)=\sum_{k=0}^{\infty}\beta \frac{-t^k}{k!}\alpha^k (\omega)=\beta exp(-\alpha t).
u(t,ω)=k=0∑∞βk!−tkαk(ω)=βexp(−αt).
通常情形下, 最多我们可以进行二 阶展开的截断, 因为对于更高阶的情形, 得到的求解系统将会变得非常复杂. 此方法被广泛地用于各 种工程领域的应用问题. 然而, 该方法还有一个固有缺陷, 即对于不确定性的放大, 因 此一般只应用于小尺度的的随机输入问题, 如小于 10% 的随机扰动.