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'); % 添加标题


 未完待续…










