Plotting Maxwellian

New to Julia
1

I have a list of positive and negative values and a single temperature. I am trying to plot the [Maxwell-Boltzmann Distribution][1] using the equation for particles moving in only one direction.

m_e = 9.11E-28 # electron mass [g]
k = 1.38E-16 # boltzmann constant [erg*K^-1]
v = range(1e10, -1e10, step=-1e8) # velocity [cm/s]

T_M = 1e6 # temperature of Maxwellian [K]

function Maxwellian(v_Max, T_Max)
    normal = (m_e/(2*pi*k*T_Max))^1.5
    exp_term = exp(-((m_e).*v_Max.*v_Max)/(3*k*T_Max))
    return normal*exp_term
end

# Initially comparing chosen distribution f_s to Maxwellian F_s
plot(v, Maxwellian.(v, T_M), label= L"F_s" * " (Maxwellian)")
xlabel!("velocity (cm/s)")
ylabel!("probability density")

However, when, plotting this, my whole function is 0:

[![enter image description here][2]][2]

I tested out if I wrote my function correctly by replacing return normal*exp_term with return exp_term (i.e. ignoring any normalization constants) and this seems to produce the distinct of the bell curve:

[![enter image description here][3]][3]

Yet, without the normalization constant, this will not preserve the area under the curve. I was wondering what may I be doing incorrectly with setting up my Maxwellian function and the constant in front of the exponential.
[1]: Maxwell–Boltzmann distribution - Wikipedia
[2]: https://i.stack.imgur.com/tYogq.png
[3]: https://i.stack.imgur.com/wyKeE.png

2

The data extrema are:

(2.9269998814999053e-123, 1.0769341115495682e-27)

Did you try manually setting the plot() keyword ylim accordingly?

3

Maybe there’s a choice of units that gives more reasonable size numbers? If so, Unitful and UnitfulRecipes could be of use.

4

Settings ylims does show the correct plot, as mentioned in my StackOverflow response to this.

Since we have the topic already, I might as well ask this here (I was planning to make a separate post): why are the default ylims of Plots.jl off by several orders of magnitude? Plotting this directly with julia> GR.plot(v, Maxwellian.(v, T_M)) automatically uses appropriate ylims, that lead to a visible plot.

Screenshot_20220214_115901

5

I think Plots has a lower threshold where it assumes whatever you’re plotting is “supposed” to be zeros. Note that eps(1.0) is about 1e-16, which could be the (somewhat flawed) reason for the specific choice. Either way, I would probably change the units to give a more intuitive scale.

6

Linking a related Plots.jl issue.

1 Like
7 
using Plots, Unitful, UnitfulRecipes

m_e = 9.11E-28u"g"
k = 1.38E-16u"erg/K" # boltzmann constant [erg*K^-1]
v = range(1e10, -1e10, step=-1e8) .* 1u"cm/s" # velocity [cm/s]

T_M = 1e6u"K" # temperature of Maxwellian [K]

function Maxwellian(v_Max, T_Max)
    normal = (m_e/(2*pi*k*T_Max))^1.5
    exp_term = exp(-((m_e).*v_Max.*v_Max)/(3*k*T_Max))
    return normal*exp_term |> upreferred
end

results in a maximum at 1.077e-21 s^3 m^-3 which is not much more reasonable size number than 1e-26s^3 cm^-3. But in any case Unitful is your friend in doing any physical calculations.

P.S. s3/m3 appears wrong units

for particles moving in only one direction"

8

Good thing about SI units: 1e-21 s^3/m^3 = 1000 ÎĽs^3/m^3

1 Like