matlab - 伊藤確率過程の符号化

Question

確率過程のルーチンを mathematica/matlab に実装しようとしています。ここに書かれたコードはすべて mathematica 用ですが、誰かがこれを matlab でエンコードするのを手伝ってくれれば (彼らがそれに慣れていれば)、それも問題ありません。ただし、可能であれば Mathematica が優先されます。

最初に方程式を調べた後、取得したいものを述べます。

以下は、私が興味を持っている伊藤確率過程です（ここでz(t)=[x(t),y(t)]:

次の量で:

また、次のものもあります (これから入力しますx(t),y(t) as x,y)。

ここで、は最初の退出時間であり 、 (最後に見てください)GT(x)の滑らかな関数ですx

目標: andだけQ0で取得したい。x

================================================== ======================================== _ 以下は mathematica のコードであることに注意してください。

終了時間の制約:

const[x_, y_] := And[10^-8 <= y <= 10^-3, 0.9*(Uc) <= a1*x^2/2 - a3*x^4/4 <= 1.1*(Uc)]

x0 = 0.1;   (* starting point for x[t] *) y0 = 0.1; (* starting point for y[t] *)

proc = ItoProcess[  {\[DifferentialD]x[t] == y[t] \[DifferentialD]t,   \[DifferentialD]y[t] == (-G*y[t] - (a1*x[t] - a3*x[t]^3) - eps*b3b*y[t]^3) \[DifferentialD]t + Sqrt[2*eps*G] \[DifferentialD]w[t]},    {t, x[t], y[t], Boole[const[x[t], y[t]]]},   {{x, y}, {x0, y0}}, {t, 0}, w \[Distributed] WienerProcess[]]

終了時刻は次のとおりです。

SeedRandom[3] 
sim = RandomFunction[proc, {0, 1, 0.001}];
First@Select[sim[[2, 1, 1]], #[[4]] == 0 &]

それが出力する{t, x[t], y[t], Boole[const[x[t], y[t]]]}

上記のコードを使用して終了時間を見つけ、それを で使用できるようにする必要がありますQ0。終了時間を見つけるための上記のスニペットは、適切に選択された に対してゼロ以外の時間を生成します。

================================================== ==============================

を見つけるために設定する必要がある数値タスクQ0(x):

--> この式の積分項から始めます。最初に、ドメイン内から始まる多くの初期条件 (たとえば、100 回の終了時間) を検索します (おそらく、この投稿で既に説明したスニペットを使用します)。これで、との関数として、100 回の終了時間ごとに積分を評価できます。x

--> での積分の期待値は、Q0以前に評価された 100 個の積分すべての標本平均 (大数の法則による) です。したがって、Q0を見つけることができ、 と の関数のみである必要がxあります。

================================================== ==============================

私の問題:

1. 上記の終了時間のコードは、適切に選択された初期条件に対してゼロ以外の終了時間を生成するように見えます。十分な終了時間が生成されるように適切なものを選択する方法について誰かが説明できれば、それはありがたいことです。上記の で指定した制約を維持したいのですがconst[x_, y_]、扱いやすい結果が得られない場合は、制約を緩和してもかまいません。

2. 以下のコードGT(x)は特異点をもたらします -- DSolve から不定式に関するエラー メッセージを受け取りました.... b0[xc] と DrhoDy[xc] の両方が不定になるため、初期条件 GT[xc] をそれも不確定になります...この問題を回避する方法は喜んでいただければ幸いです。

3. 最後に、Mathematica で 100 の積分を効率的に評価して (項が巨大であるため)、以前に見つけた各終了時間と期待値を取るために誰かの助けが本当に必要です。出口時間を正しく見つける方法がわかりません。

================================================== ============================== GT(x)を見つけるには:

b1b = 0.9; b3b = .8; a1b = 0.1; a3b = 0.2; G = (1/0.1^2)*
  b1b; a1 = (1/.1^2)*a1b; a3 = (1/.1^2)*a3b; xc = Sqrt[a1/a3]; Uc = 
 a1*xc^2/2 - a3*xc^4/4; L = (G + Sqrt[G^2 + 4*(-(a1 - 3*a3*xc^2))])/2;
y[x_] := (x *a1 - x^3*a3)/(G + L)

b0[x_] := (y[
      x]*(a1*x \[Minus] 
       a3*x ^3)*(1 \[Minus] (a1 \[Minus] 3*a3*x ^2)) \[Minus] 
    G*y[x]*(a1 \[Minus] 3*a3*x ^2)) /(y[
      x] ^2 + (\[Minus]G*y[x] + a1*x \[Minus] a3*x ^3) ^2)


DrhoDy [x_] := 
 Sqrt[y[x]^ 2 /(y[x] ^2 + (\[Minus]G*y[x] + a1*x \[Minus] a3*x ^3) ^2)]


linearequation = 
  y[x]*GT '[x] + b0[x]*GT[x] == G*DrhoDy [x]^2 *GT[x]^ 3;

GT[x] = 
 DSolve[{linearequation, GT[xc] == Sqrt[b0[xc]/( G*DrhoDy [xc]^2 )]}, 
  GT, x]

ODE:

次の初期条件を満たします。

どこ、

と