ode45 для уравнения Ланжевена

Question

У меня вопрос об использовании Matlab для вычисления решения стохастических дифференциальных уравнений. Уравнения - это 2.2a, b, стр. 3, в 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. Глядя на бумагу, с которой вы связались, они использовали базовую схему для интеграции системы. Это очень простой решатель, чтобы реализовать себя.

Прежде чем продолжить, вы (и ваш профессор!) Должны занять некоторое время, чтобы прочитать SDE и как их решить численно. Я рекомендую этот документ, который включает в себя множество примеров Matlab: (.. Образа Sect)

Desmond J. Higham, 2001, Алгоритмическое Введение в Численное моделирование стохастических дифференциальных уравнений, SIAM Rev., 43 525- 46. http://dx.doi.org/10.1137/S0036144500378302

URL-адрес файлов Matlab в документе не будет работать; use this one. Обратите внимание, что поскольку этот 15-летний документ, некоторые из кода, связанные с генерацией случайных чисел, устарели (используйте для использования в качестве генератора генератор rng(1) вместо randn('state',1)).

Если вы знакомы с ode45, вы можете посмотреть мой SDETools Matlab toolbox на GitHub. Он был разработан так, чтобы быть быстрым и имеет интерфейс, который очень похож на набор ODE от Matlab. Вот как вы можете закодировать ваш пример с использованием решателя Эйлера-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, чтобы попытаться справиться с этим, или вам, возможно, потребуется создать свой собственный решатель. Я также не мог найти, какой промежуток времени интеграции используется в документе. Вы также должны обратиться к авторам статьи напрямую.

UPDATE

Вот реализация Эйлера-Maruyama, который соответствует sde_euler выше код:

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');