3

確率微分方程式の解を計算するための Matlab の使用について質問があります。方程式は、この論文(PDF)の 2.2a、b、ページ 3です。

私の教授は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 は大きく異なり、SDEode45を解くために などの ODE 用のツールを使用しないでください。リンク先の論文を見ると、彼らは基本的なオイラー・マルヤマ方式を使用してシステムを統合していました。これは、自分で実装するための非常に単純なソルバーです。

先に進む前に、あなた (そしてあなたの教授も!) は、SDE とそれらを数値的に解く方法について読む必要があります。多くの Matlab の例が含まれているこの論文をお勧めします。

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

論文の Matlab ファイルへの URL は機能しません。これを使用してください。これは 15 年前の論文であるため、乱数生成に関連するコードの一部が古くなっていることに注意してください (ジェネレーターをシードするrng(1)代わりに使用してください)。randn('state',1)

よく知っている場合は、GitHubode45にある私の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 に答える