Ejercicio St

Stochastic Simulation Methods Rebeca de la Fuente Mara˜ no´n November 3, 2015

√ We have he differential equation x˙ = −ax + Dξ(t) where ξ(t) is a Gaussian white noise of zero mean and correlation hξ(t)ξ(t0 )i = δ(t − t0 ) For a=2, D=0.05 and the initial condition x(0)=4 we integrate the equation. To do this we use the Milsthein algorithm discretising the time: x(ti+1 ) = x(ti )(1 − ha) +

√

DhUi

where Ui is a Gaussian random variable that measures the related noise for every step in our algorithm, that is, we calculate a new random number every time, which have a Gaussian distribution. We find trajectories with fluctuation because of the white noise, that si the derivative of the Wiener process. At first, we plot only one trajectory to realized how the fluctuations behave:

In this plot we can see ten diferent trajectories:

1

The deterministic term of the differential equation is x˙ = −ax. For the same conditions than before, we find as a solution x(t) = 4 exp−2t , and defines the deterministic dynamics. We can see this green trajectory without fluctuations in this plot:

To identify an initial dynamical regime dominated by deterministic dynamics, I’ve imposed in my algorithm different conditions which break the generation of every trajectory and gives me the values at which the deterministic equation domains. I’ve done it for three different contitions to search the values: For the condition | 4e−2t − x(t) |< 0.001 7.81 12.19 12.2 10.25 10.26 13.15 5.99 7.01 7.39 1.67 1.68 1.69 1.59 2.61 2.62 6.17 11.22 13.84 6.44 6.45 6.75 We can realized that the error is very small. Thus, the times at which this condition is accomplished is very different among trajectories.Furthermore, we can see in one arbitrary trajectory that sometimes the condition occurs for times that are away from each other. Our condition is too much strict. Anyway, if we do the average between trajectories at the first time step in which the condition is accomplished we find that the initial dynamical regime searched is at time: 5.70. This is not very much realistic. For the condition | 4e−2t − x(t) |< 0.01 2.98

2.45

1.7

2.67

3.55

2.02

2.93 2

2.22

2.65

2.55

Every time corresponds to the first time step of different trajectories at which the condition is acomplished. Now we can see that this is more realistic, that is, out condition gives to us a more definited regime, from 1.7 to 3.55. The average is 2.512. If we compare it with the plot, it has sense. For the condition | 4e−2t − x(t) |< 0.1 1.64 1.83 1.67 1.43 1.88 2.07 1.77 Here we find values smaller than in the previous condition. That is because of the error. Now we have an error bigger, and then the distance between the times at which this occurs is smaller. Conclusion We find the second condition as the better one, so I identify an initial dynamical regime in [2.02,2.98] in which the deterministic dynamics dominate because the times 1.7 and 3.55 are too much away from the mean. The stationary regime in which the system fluctuates around the steady state is [-0.5,0.5] such that x(t) ∈ [−0.5, 0.5]. We can see it in the plot. How does the time depend on the parameter a? We can see now the plot of trajectories for different values of the parameter a and for different initial conditions in the same algorithm implemented by c++ language. For a=0.25 and x(0)=4

For a=0.5 and x(0)=4

3

For a=2 and x(0)=4

For a=4 and x(0)=4

4

For a=2 and x(0)=0

For a=2 and x(0)=0.5

5

For a=2 and x(0)=1

Since a is the parameter in the deterministic term, it dominates more with high values, and the fluctuations becomes smaller. So the trajectories with small values of a fluctuates very much, and the stationary regime is not so defined. We have seen before the time it takes to reach the stationary state and for higher values of a this time is smaller. Altought for enough high values of h the differential equation becomes more deterministic (the deterministic term dominates), the noise will always perturbate the system. The equation 6

decays exponentially with time, and it do it faster and the white noise becomes smoother when a is bigger. From numerical integration we evaluate hx(t)i, where h...i stands for averages over trajectories, with initial condition x(0)=1 Over 10 trajectories

Over 100 trajectories

7

Over 1000 trajectories

From numerical integration we evaluate hx(t)2 i, with the same conditions than before Over 10 trajectories

Over 100 trajectories 8

Over 1000 trajectories

Conclusions As we can see hx(t)i becomes constant at zero (because of the negative and positive values) with the number of trajectories and the noise becomes irrelevant in the equation. So it becomes deterministic. And hx(t)i can only be positive, so we can see the same behaviour but for values grather than zero: the limit to this term remains displazed of the last one and more deterministic since it seems more a straight line, a constant once reaches the stationary state. 9

Correlation function as function of s C(t, x) = hx(t)x(t + s)i for a=2, D=0.05 and x(0)=1 averaging over 50000 trajectories. s ∈ [0, 10]. We plot C(t,s) for t=1,t=2,t=5,t=10 and t=20

For t grather or equal to 2 we find C(t,s)=C(s), so the stationary regime in which this happens is [2, ∞) because C(t,s) is equal in t=2, t=5, t=10 and t=20, but not in t=1. This is related with the previous conclusions of our differential equation, where a dominates the equation for values sufficiently high. The correlation function C(t, s) = hx(t)x(t + s)it over time in a single long trajectory, with a=2 and D=0.05 from x=0 We plot this function together with C(s)

10

The correlation functions indicates that the diferences between two values at diffeent times of one trajectory behaves as the diferences between values at diferents times in all trajectories. Since we find that both correlation functions behaves equal, this function does not depend on t. Altought there are some diferences in the plot, this happens because it must been ploted with more precision of s. The discretitzacion could be better. PROGRAMS Integration of the differential equation

#include #include #include #include #include #define pi 3.14159265358979323846 using namespace std; int main(void) { double s,k,u; double a=2; double D; double h=0.01; double x[100000],y[100000]; D=0.05; 11

y[0]=4; x[0]=0; srand(time(NULL)); for (int i=1; i