3

我有一个关于使用 Matlab 计算随机微分方程解的问题。方程是本文(PDF)中第 3 页的 2.2a,b 。

我的教授建议使用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')

再次感谢您的帮助 ！

4

1 回答 1

1

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

在继续之前，您（和您的教授！）应该花一些时间阅读 SDE 以及如何以数字方式求解它们。我推荐这篇论文，其中包含许多 Matlab 示例：

Desmond J. Higham，2001，随机微分方程数值模拟的算法介绍，SIAM Rev.（Educ. Sect.），43 525–46。http://dx.doi.org/10.1137/S0036144500378302

论文中 Matlab 文件的 URL 不起作用；用这个。请注意，由于这是一篇已有 15 年历史的论文，一些与随机数生成相关的代码已经过时（使用rng(1)而不是randn('state',1)生成生成器的种子）。

如果你熟悉ode45你可以看看我在 GitHub 上的SDETools Matlab 工具箱。它的设计速度很快，并且具有与 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否则噪音太大，可能导致半径变量变为负数。我不确定这篇论文是否讨论了他们如何处理这个奇点。您可以尝试其中的'NonNegative'选项sdeset来尝试处理此问题，或者您可能需要构建自己的求解器。我也找不到论文使用的积分时间步长。您还应该考虑直接联系论文的作者。

更新

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');
于 2016-06-08T15:08:44.573 回答