# Monte Carlo Simulation of a Double Well Potential (Edited)

I am looking for a possible way to simulate the temperature effect on a collection of atoms sitting within a double well potential, say something like

V[x_]:=-2 x^2 + 2 x^4 + 0.5


giving a potential looking like the following

I have written the following code to simulate the random population versus temperature :

sampleplot2D[n_, \[Lambda]_] :=
Graphics[Block[{samplex, sampley, s},
Table[samplex = RandomReal[{-1.1, 1.1}];
sampley = RandomVariate[ExponentialDistribution[\[Lambda]]];
If[sampley > 2 samplex^4 - 2 samplex^2 + 0.5, s = Red, s = Blue];
{s, Point[{samplex, sampley}]}, {n}]],
PlotRange -> {{-1.1, 1.1}, {-0, 1}}, AspectRatio -> 1, Frame -> True]

Histplot2D[n_, \[Lambda]_] :=
Block[{samplex, sampley, datax = {}, datay = {}},
Table[samplex = RandomReal[{-1.1, 1.1}];
sampley = RandomVariate[ExponentialDistribution[\[Lambda]]];
If[sampley >
2 samplex^4 - 2 samplex^2 + 0.5, {datax = Append[datax, samplex],
datay = Append[datay, sampley]}], {n}];
Histogram[datax, 50]]

Manipulate[
GraphicsRow[{sampleplot2D[5000, \[Lambda]],
Histplot2D[1000, \[Lambda]]}, ImageSize -> Large], {\[Lambda], 1,
25, 2}]


In this case I have assumed a random distribution in the x axis (can be a position), while along the y-axis (energy) I have assumed an exponentially decaying probability, where \[Lambda] is the parameter of the exponential decay and can be seen as the inverse of the temperature. For high \[Lambda] the population is localized in the bottom of the potential well; while for low \[Lambda] filling it more.

Q1 : I see that running the Manipulate function the code is running twice per every lambda selected, this is probably due to my limited understanding of local variables in the Block function and the fact that I have 2 separate functions, for the 2D plot and the histogram.

Is there a way to have a similar output using only 1 function? Can the code be optimized? Already now going beyond n = 10^4 makes the calculation very heavy.

Q2: Now assuming to fill only the right side of the potential (changing samplex = RandomReal[{-1.1, 1.1}]; into samplex = RandomReal[{0, 1.1}];, I am trying to apply an approach as discussed in How to set up a simple Monte Carlo Simulation, to randomize even further in time the position of the n atoms. In practice what I want to see is the amount of redistribution into the left side starting from the right side population. At the moment I was not able to successfully merge the codes into a working set.

• Matteo, StackExchange is not a free coding service. Do you have any concrete questions about what you've tried? May 22 '19 at 11:12
• Hi Roman, sure it wasn't intended as a request of "code for me please". I did not attach anything because at the moment I wasn't able to make anything working, not even remotely. I will try harder for a couple of days and update my post. Sorry for the inconvenience. Best May 22 '19 at 16:27
• I suggest to learn some functional programming reference.wolfram.com/language/howto/WorkWithPureFunctions.html and use functions like Gather , Select , Map and so on May 29 '19 at 9:45

samplot2D[n_, lam_]:=

RandomVariate[#,n]& /@ {
UniformDistribution@{-1.1, 1.1},
ExponentialDistribution@lam} //

, #2 > 2 #1^4 - 2 #1^2 + 0.5 & @@ # &] & //


There is no need for auxiliary variables and procedural commands like Table, If etc.
• @MatteoS sorry I do not understand what you mean. GatherBy splits the list of random points {{x1,y1},{x2,y2}, ..} , in two sublists; the first sublist consists of points such that their coordinates satisfy the inequality, the other sublist does not satisfy the inequality. May 29 '19 at 13:34
• Oh I see. You could use GroupBy instead of GatherBy then, so that you have control over which list satisfies the inequality. May 29 '19 at 16:21