2
$\begingroup$

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

example of potential curve

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.

Thanks a lot for any help you can provide!

Improve this question
$\endgroup$
3
  • 3
    $\begingroup$ Matteo, StackExchange is not a free coding service. Do you have any concrete questions about what you've tried? $\endgroup$
    – Roman
    May 22, 2019 at 11:12
  • $\begingroup$ 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 $\endgroup$
    – Matteo S
    May 22, 2019 at 16:27
  • 1
    $\begingroup$ I suggest to learn some functional programming reference.wolfram.com/language/howto/WorkWithPureFunctions.html and use functions like Gather , Select , Map and so on $\endgroup$
    – Fortsaint
    May 29, 2019 at 9:45

1 Answer 1

Reset to default
2
$\begingroup$

This is a restyled code you can start with

samplot2D[n_, lam_]:=

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

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

    Graphics[Thread@{{Red, Blue}, Point/@#}
             , PlotRange -> {{-1.1, 1.1}, {-0, 1}}
             , AspectRatio -> 1
             , Frame -> True]&

There is no need for auxiliary variables and procedural commands like Table, If etc.

Improve this answer
$\endgroup$
4
  • $\begingroup$ Thanks a lot @Fortsaint! that is incredibly more efficient than my code. I am not used with the short notation (but I am studying now the link you posted above). But do I understand it correctly that the gathering step is done by the first point? I mean that if the first point satisfies the inequality than the gathering and subsequent plotting shows in red all points within the double well. Instead if the first random point does not satisfies, then the color will be inverted. Is this correct? $\endgroup$
    – Matteo S
    May 29, 2019 at 13:02
  • $\begingroup$ @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. $\endgroup$
    – Fortsaint
    May 29, 2019 at 13:34
  • $\begingroup$ It is just that if you run several time samplot2D, sometimes the colors get inverted, thus I assumed that GatherBy was indeed dividing the list into the 2 subsets but that the order between satisfying and not satisfying was determined by the first random point $\endgroup$
    – Matteo S
    May 29, 2019 at 14:27
  • $\begingroup$ Oh I see. You could use GroupBy instead of GatherBy then, so that you have control over which list satisfies the inequality. $\endgroup$
    – Fortsaint
    May 29, 2019 at 16:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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