matlab - 朗之万方程的 ode45

Question

我有一个关于使用 Matlab 计算随机微分方程解的问题。方程式是第 3 页的 2.2a、b，在 this paper 中(PDF)。

我的教授建议使用时间步长较小的ode45，但结果与文章中的不符。特别是时间序列和pdf。我对函数中白噪声的定义也有疑问。

这里是集成函数的代码:

function dVdt = R_Lang( t,V )

global  sigma lambda alpha

W1=sigma*randn(1,1); 
W2=sigma*randn(1,1);
dVdt=[alpha*V(1)+lambda*V(1)^3+1/V(1)*0.5*sigma^2+W1;
        sigma/V(1)*W2];

end

主要脚本:

clear variables
close all
global  sigma lambda alpha
sigma=sqrt(2*0.0028);
alpha=3.81;
lambda=-5604;

tspan=[0,10];
options = odeset('RelTol',1E-6,'AbsTol',1E-6,'MaxStep',0.05);

A0=random('norm',0,0.5,[2,1]);
[t,L]=ode45(@(t,L) R_Lang(t,L),tspan,A0,options);

如果您有任何建议，我将不胜感激。

这里是对抗我的 EM 方法和“sde_euler”的新代码。

lambda = -5604; 
sigma=sqrt(2*0.0028) ; 
Rzero = 0.03;    % problem parameters
phizero=-1;
dt=1e-5;
T = 0:dt:10;
N=length(T);
Xi1 = sigma*randn(1,N);         % Gaussian Noise with variance=sigma^2
 Xi2 = sigma*randn(1,N);
alpha=3.81;

Rem = zeros(1,N);                 % preallocate for efficiency
Rtemp = Rzero;
phiem = zeros(1,N);                 % preallocate for efficiency
phitemp = phizero;
for j = 1:N
     Rtemp = Rtemp + dt*(alpha*Rtemp+lambda*Rtemp^3+sigma^2/(2*Rtemp)) + sigma*Xi1(j);
   phitemp=phitemp+sigma/Rtemp*Xi2(j);
   phiem(j)=phitemp;
   Rem(j) = Rtemp;

end

f = @(t,V)[alpha*V(1)+lambda*V(1)^3+0.5*sigma^2/V(1)/2;
           0];                 % Drift function
g = @(t,V)[sigma;
           sigma/V(1)];        % Diffusion function
A0 = [0.03;0];                % 2-by-1 initial condition
opts = sdeset('RandSeed',1,'SDEType','Ito'); % Set random seed, use Ito formulation
L = sde_euler(f,g,T,A0,opts);       

plot(T,Rem,'r')
hold on
plot(T,L(:,1),'b')

再次感谢您的帮助!

Answer 1

ODE 和 SDE 非常不同，不应使用 ODE 工具(例如 ode45)来尝试求解 SDE。查看您链接到的论文，他们使用了基本的 Euler-Maruyama系统集成方案。这是一个非常简单的求解器，可以自己实现。

在继续之前，您(和您的教授!)应该花一些时间阅读 SDE 以及如何用数值求解它们。我推荐这篇论文，其中包含许多 Matlab 示例:

Desmond J. Higham, 2001, An Algorithmic Introduction to Numerical Simulation of Stochastic Differential Equations, SIAM Rev. (Educ. Sect.), 43 525–46. http://dx.doi.org/10.1137/S0036144500378302

论文中 Matlab 文件的 URL 无效； use this one .请注意，由于这是一篇已有 15 年历史的论文，一些与随机数生成相关的代码已过时(使用 rng(1) 而不是 randn('state',1 ) 为生成器播种)。

如果您熟悉 ode45，您可能会看看我的 SDETools Matlab toolbox在 GitHub 上。它的设计速度很快，并且具有与 Matlab 的 ODE 套件非常相似的界面。以下是使用 Euler-Maruyma 求解器编写示例代码的方法:

sigma = 1e-1*sqrt(2*0.0028);
lambda = -5604;
alpha = 3.81;
f = @(t,V)[alpha*V(1)+lambda*V(1)^3+0.5*sigma^2/V(1);
           0];                 % Drift function
g = @(t,V)[sigma;
           sigma/V(1)];        % Diffusion function
dt = 1e-3;                     % Time step
t = 0:dt:10;                   % Time vector
A0 = [0.03;-2];                % 2-by-1 initial condition
opts = sdeset('RandSeed',1,'SDEType','Ito'); % Set random seed, use Ito formulation
L = sde_euler(f,g,t,A0,opts);                % Integrate

figure;
subplot(211);
plot(t,L(:,2));
ylabel('\phi');
subplot(212);
plot(t,L(:,1));
ylabel('r');
xlabel('t');

我不得不减小 sigma 的大小，否则噪声太大以至于可能导致半径变量变为负数。我不确定这篇论文是否讨论了他们如何处理这种奇点。您可以尝试 sdeset 中的 'NonNegative' 选项来尝试处理此问题，或者您可能需要构建自己的求解器。我也找不到论文使用的积分时间步长。您还应该考虑直接联系论文的作者。

更新

这是一个与上面的 sde_euler 代码匹配的 Euler-Maruyama 实现:

sigma = 1e-1*sqrt(2*0.0028);
lambda = -5604;
alpha = 3.81;
f = @(t,V)[alpha*V(1)+lambda*V(1)^3+0.5*sigma^2/V(1);
           0];                 % Drift function
g = @(t,V)[sigma;
           sigma/V(1)];        % Diffusion function
dt = 1e-3;                     % Time step
t = 0:dt:10;                   % Time vector
A0 = [0.03;-2];                % 2-by-1 initial condition

% Create and initialize state vector (L here is transposed relative to sde_euler output)
lt = length(t);
n = length(A0);
L = zeros(n,lt);
L(:,1) = A0;

% Set seed and pre-calculate Wiener increments with order matching sde_euler
rng(1);
r = sqrt(dt)*randn(lt-1,n).';

% General Euler-Maruyama integration loop
for i = 1:lt-1
    L(:,i+1) = L(:,i)+f(t(i),L(:,i))*dt+r(:,i).*g(t(i),L(:,i));
end

figure;
subplot(211);
plot(t,L(2,:));
ylabel('\phi');
subplot(212);
plot(t,L(1,:));
ylabel('r');
xlabel('t');