简单随机微分方程数值解

news2024/9/21 20:38:58

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求解

考虑一个由下列方程支配的随机激励的单自由度振子：
$\ddot{X}+h(X,\dot{X})=XW_{1}(t)+\dot XW_{2}(t)+W_{3}(t)$
其中， $h(X,\dot X)$ 表示阻尼力和恢复力， $W_{l}(t),l=1,2,3$ 是独立的高斯白噪声，其相关函数为 $E[W_{l}(t)W_{s}(t+\tau)]=2\pi K_{ls} \delta_{ls}\delta(\tau)$
这里， $\delta_{ls}$ 与 $\delta(\tau)$ 不一样，前者是克罗内克(Kronecker) $\delta$ ，即 $\delta_{ls}=1,l=s;\delta_{ls}=0,l\neq s.$ ， $\delta(\tau)$ 是狄拉克函数， $\delta=\infty, \tau=0;\delta=0,\tau \neq 0.$

记 $X_{1}=X,X_{2}=\dot{X}$ 可转换状态空间的两个方程
$\dot{X}_{1}=X_{2}$
$\dot{X}_{2}=-h(X_{1},X_{2})+X_{1}W_{1}(t)+X_{2}W_{2}(t)+W_{3}(t)$

对于 $\frac{d}{dt}X(t)=f(X,t)+g(X,t)W(t)$
相应的FPK方程为 $\frac{\partial }{\partial t}p+\sum_{j=1}^{n}\frac{\partial }{\partial x _{j}}(a_{j}p)-\frac{1}{2}\sum_{j,k=1}^{n}\frac{\partial ^{2}}{\partial x_{j}\partial x_{k}}(b_{jk}p)=0$
其中，一、二阶导数矩为 $a_{j}(\mathbf{x},t)=f_{j}(\mathbf{x},t)+\pi\sum_{r=1}^{n}\sum_{l,s=1}^{m}K_{ls}g_{rs}(\mathbf{x},t)\frac{\partial }{\partial x_{r}}g_{jl}(\mathbf{x},t)$ ,
$b_{jk}(\mathbf{x},t)=2\pi\sum_{l,s=1}^{m}K_{ls}g_{jl}(\mathbf{x},t)g_{ks}(\mathbf{x},t)$

由此可以得出, $f_{1}=x_{2},f_{2}=-h(x_{1},x_{2})$ , $g_{1j}=0(j=1,2,3)$
$g_{21}=x_{1},g_{22}=x_{2},g_{23}=1$ , $n = 2, m = 3$ .
则一、二阶导数矩：
$a_{1}=x_{2},a_{2}=-h(x_{1},x_{2})+\pi K_{22}x_{2}$
$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}.$
从而得到FPK方程：
$\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$ .

相应的也可以得到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'); % 添加标题