Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I want to run an Ito stochastic process. I have the following parameters

b1b = 0.9;
b3b = .8;
a1b = 0.1;
a3b = 0.2;
eps = 0.1;
G = (1/eps^2)*b1b ; a1 = (1/eps^2)*a1b; a3 = (1/eps^2)*a3b;
xc = Sqrt[a1/a3];
Uc = a1*xc^2/2 - a3*xc^4/4

I want to stop the process when U[x[t]]=U[xc] and y[t]=0 (or at least very close for each of these parameters: say in the vicinity of 10^-4*U[xc] for U[x[t]] and 10^-8 for y[t]). Please note U[x] is the integral of U'[x].

I am trying to write code that simulates dz[t] for many realisations (say 100) and calculates the average of times at which the aforementioned constraint is satisfied (i.e., average of 100 times). However, I have no idea how to incorporate the constraints and extract the time at which the constraints are satisfied.

I tried to use the ItoProcess[] function, but had no fruitful outcome. Any help would be much appreciated.

share|improve this question
    
Have you seen this demonstration project? –  Rod Jul 30 '13 at 0:33
add comment

1 Answer

up vote 2 down vote accepted

Not sure I got your process right and not a very elegant solution : setup the process, do the simulations first, then look at the outcome. First off define your constraints :

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

I am outputting

{t, x[t], y[t], Boole[const[x[t], y[t]]]}

which contains all the information we need.

Your process :

x0 = 0.35;   (* starting point for x[t] *)
y0 = 0.0005; (* 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[]]

A sample run of the simulation : Q1

SeedRandom[3]
sim = RandomFunction[proc, {0, 1, 0.001}];

ListLinePlot[{sim[[2, 1, 1, All, {1, 2}]], sim[[2, 1, 1, All, {1, 3}]]}] 

plot

Find the first exit time :

First@Select[sim[[2, 1, 1]], #[[4]] == 0 &]
(* {0.001, 0.350001, 0.158739, 0} *)

So for this specific run your process exited at the first step past the initial condition. This is not very surprising since your constraints are quite tight; for instance : Q2

Reduce[const[x, y][[2]], x, Reals]/. Or -> List
(* { -1.37126 <= x <= -1.36069, 
    -0.385397 <= x <= -0.345918, 
     0.345918 <= x <= 0.385397, 
      1.36069 <= x <= 1.37126} *)
share|improve this answer
1  
I double checked and the code above seems to work; please try with a fresh kernel. –  b.gatessucks Jul 30 '13 at 10:40
    
If the above code works then yes, you can add the rest of your constraints. –  b.gatessucks Jul 30 '13 at 11:26
    
@b.gatessucks It could be that he is trying to run your code directly, without previously inputing his own values... Try to edit your code putting the assigned values. By the way, the code works perfectly for me. +1 ! –  Rod Jul 30 '13 at 18:51
    
Please see edit. –  b.gatessucks Jul 31 '13 at 6:57
    
@b.gatessucks can you kindly look at my other question here: mathematica.stackexchange.com/questions/29770/… I have explained the logic properly....please let me know if something is unclear...i'll try to make it clear ASAP....thank you once again for this and your help on this question would be much appreciated....Regards! –  Stoc Aug 5 '13 at 2:45
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.