1.随机微分方程求解：dX(t) =− αXtdt + σdWt
法一：Euler-Maruyama
%%
%O-U过程
%dX(t)=-alpha*Xt*dt+sigma*dWt,X|t=0=X0
%alpha=2,sigma=1,X0=1
% 设置初始参数
T = 1; % 时间区间长度
N = 1000; % 离散化的时间步数
dt = T/N; % 时间步长
X = zeros(1,N+1); % 存储解向量
X(1) = 1; % 初始条件
alpha = 2;
sigma = 1;
% 模拟数值解
for i = 1:N
dW = sqrt(dt)*randn(); % 标准正态分布增量
X(i+1) = X(i) - alpha*X(i)*dt + sigma*dW; % 欧拉方法更新
end
% 绘图
plot(linspace(0,T,N+1),X,'*-') % 根据时间步长将x轴离散化
xlabel('t')
ylabel('X(t)')
title('随机微分方程数值解')
legend('Euler数值解')
法二：Milstein
% 设置参数
alpha = 2;
sigma = 1;
X0 = 1;
T = 1;
N = 1000;
dt = T/N;
% 初始化向量和随机项
t = 0:dt:T;
W = [0,cumsum(randn(1,N).*sqrt(dt))];
%使用cumsum函数生成一个与t等长的Wiener过程（随机项）
% 初始化X
X = zeros(1,N+1);
X(1) = X0;
% Milstein方法计算X
for i = 1:N
dW = W(i+1) - W(i);
X(i+1) = X(i) - alpha*X(i)*dt + sigma*X(i)*dW ...
+ 0.5*sigma^2*X(i)*(dW^2-dt);
%在Milstein方法中，我们需要对二次变差数进行估计
%因此在计算时需要添加正交项0.5*sigma^2*X(i)*(dW^2-dt)
end
% 绘制图形
plot(t,X,'*-')
title('随机微分方程数值解')
xlabel('t')
ylabel('X(t)')
legend('Milstein数值解')
2.受高斯白噪声激励的系统，FPK求解
考虑一个由下列方程支配的随机激励的单自由度振子：
X
¨
+
h
(
X
,
X
˙
)
=
X
W
1
(
t
)
+
X
˙
W
2
(
t
)
+
W
3
(
t
)
\ddot{X}+h(X,\dot{X})=XW_{1}(t)+\dot XW_{2}(t)+W_{3}(t)
X¨+h(X,X˙)=XW1(t)+X˙W2(t)+W3(t)
其中，
h
(
X
,
X
˙
)
h(X,\dot X)
h(X,X˙)表示阻尼力和恢复力，
W
l
(
t
)
,
l
=
1
,
2
,
3
W_{l}(t),l=1,2,3
Wl(t),l=1,2,3是独立的高斯白噪声，其相关函数为
E
[
W
l
(
t
)
W
s
(
t
+
τ
)
]
=
2
π
K
l
s
δ
l
s
δ
(
τ
)
E[W_{l}(t)W_{s}(t+\tau)]=2\pi K_{ls} \delta_{ls}\delta(\tau)
E[Wl(t)Ws(t+τ)]=2πKlsδlsδ(τ)
这里，
δ
l
s
\delta_{ls}
δls与
δ
(
τ
)
\delta(\tau)
δ(τ)不一样，前者是克罗内克(Kronecker)
δ
\delta
δ，即
δ
l
s
=
1
,
l
=
s
;
δ
l
s
=
0
,
l
≠
s
.
\delta_{ls}=1,l=s;\delta_{ls}=0,l\neq s.
δls=1,l=s;δls=0,l=s.，
δ
(
τ
)
\delta(\tau)
δ(τ)是狄拉克函数，
δ
=
∞
,
τ
=
0
;
δ
=
0
,
τ
≠
0.
\delta=\infty, \tau=0;\delta=0,\tau \neq 0.
δ=∞,τ=0;δ=0,τ=0.
记
X
1
=
X
,
X
2
=
X
˙
X_{1}=X,X_{2}=\dot{X}
X1=X,X2=X˙可转换状态空间的两个方程
X
˙
1
=
X
2
\dot{X}_{1}=X_{2}
X˙1=X2
X
˙
2
=
−
h
(
X
1
,
X
2
)
+
X
1
W
1
(
t
)
+
X
2
W
2
(
t
)
+
W
3
(
t
)
\dot{X}_{2}=-h(X_{1},X_{2})+X_{1}W_{1}(t)+X_{2}W_{2}(t)+W_{3}(t)
X˙2=−h(X1,X2)+X1W1(t)+X2W2(t)+W3(t)
由此可以得出,
f
1
=
x
2
,
f
2
=
−
h
(
x
1
,
x
2
)
f_{1}=x_{2},f_{2}=-h(x_{1},x_{2})
f1=x2,f2=−h(x1,x2),
g
1
j
=
0
(
j
=
1
,
2
,
3
)
g_{1j}=0(j=1,2,3)
g1j=0(j=1,2,3)
g
21
=
x
1
,
g
22
=
x
2
,
g
23
=
1
g_{21}=x_{1},g_{22}=x_{2},g_{23}=1
g21=x1,g22=x2,g23=1,
n
=
2
,
m
=
3
n=2,m=3
n=2,m=3.
则一、二阶导数矩：
a
1
=
x
2
,
a
2
=
−
h
(
x
1
,
x
2
)
+
π
K
22
x
2
a_{1}=x_{2},a_{2}=-h(x_{1},x_{2})+\pi K_{22}x_{2}
a1=x2,a2=−h(x1,x2)+πK22x2
b
11
=
0
,
b
12
=
0
,
b
21
=
0.
b
22
=
2
π
K
11
x
1
2
+
2
π
K
22
x
2
2
+
2
π
K
33
.
b_{11}=0,b_{12}=0,b_{21}=0.b_{22}=2\pi K_{11}x_{1}^{2}+2\pi K_{22}x_{2}^{2}+2\pi K_{33}.
b11=0,b12=0,b21=0.b22=2πK11x12+2πK22x22+2πK33.
从而得到FPK方程：
∂
∂
t
p
+
∂
∂
x
1
(
x
2
p
)
+
∂
∂
x
2
{
[
(
−
h
(
x
1
,
x
2
+
π
K
22
x
2
]
p
}
−
π
∂
2
∂
x
2
2
[
(
K
11
x
1
2
+
K
22
x
2
2
+
K
33
)
p
]
=
0
\frac{\partial}{\partial t}p+\frac{\partial}{\partial x_{1}}(x_{2}p)+\frac{\partial}{\partial x_{2}}\left \{ [(-h(x_{1},x_{2}+\pi K_{22}x_{2}]p \right \}-\pi \frac{\partial ^{2}}{\partial x_{2}^{2}}[(K_{11}x_{1}^{2}+K_{22}x_{2}^{2}+K_{33})p]=0
∂t∂p+∂x1∂(x2p)+∂x2∂{[(−h(x1,x2+πK22x2]p}−π∂x22∂2[(K11x12+K22x22+K33)p]=0.
相应的也可以得到Stratonovich或Ito方程，它们得到的FPK方程都和上式一致。
3.绘制随机微分方程概率密度曲线
%(b)
clc;clear;
% 定义需要绘制的变量和参数
t_values = linspace(0.01, 5, 100); % 在t轴上均匀采样100个点，保证不会出现除数为0的情况
x_values = linspace(-5, 5, 200); % 在x轴上均匀采样200个点
[x_mesh, t_mesh] = meshgrid(x_values, t_values); % 创建网格
% 计算概率密度函数
p = (1./sqrt(2*pi.*t_mesh)).*cosh(x_mesh).*exp(-0.5./t_mesh).*exp(-(x_mesh.^2)./(2*t_mesh));
% 绘制演化图
mesh(x_mesh, t_mesh, p); % 使用surf函数绘制三维曲面图
xlabel('x'); ylabel('t'); zlabel('p(x,t)'); % 添加轴标签
title('Probability Density Evolution'); % 添加标题
未完待续…
