# How to make this distribution symmetric?

**URL:** <https://discourse.julialang.org/t/how-to-make-this-distribution-symmetric/76728>\
**Category:** New to Julia\
**Tags:** distributions\
**Created:** [February 19, 2022, 4:57am UTC](https://discourse.julialang.org/t/how-to-make-this-distribution-symmetric/76728 "2022-02-19T04:57:46Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![CentraCep](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/centracep/32/32454_2.png) [@CentraCep](https://discourse.julialang.org/u/CentraCep)\
**Post date:** [February 19, 2022, 4:57am UTC](https://discourse.julialang.org/t/how-to-make-this-distribution-symmetric/76728/1 "2022-02-19T04:57:46Z")

</div>

Let’s say I have the following data file:

https://ufile.io/b2v6t45e

Using this file, I can plot a Maxwell-Boltzmann functon and see how it changes as a function of `s` with this code:

````julia
# using necessary packages
using LaTeXStrings, Markdown, DataFrames, Queryverse, Plots, Printf
 
df = load("data.csv", spacedelim=true, header_exists=true) |> DataFrame

m_e = 9.11E-28 # electron mass [g]
k = 1.38E-16 # boltzmann constant [erg*K^-1]

function Maxwellian(v_Max, T_Max)
    ```Maxwellian distribution function for particles moving in only one direction```
    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

nu = 1 # collisional frequency [s^-1]
v_pos = range(0, 1e10, step=1e8) # positive velocity [cm/s]

f_s_pos = [Maxwellian.(v_pos, df.T[1])] # because boundary density is high at left footpoint, assume f_N = f_M

for is in 1:length(df.s) # for positive velocities, integrate from 0 to s_max
    local_F_s_pos = Maxwellian.(v_pos, df.T[is])
    push!(f_s_pos, f_s_pos[is] - (nu./v_pos).*(f_s_pos[is] - local_F_s_pos)*(df.ds[is]))
end

anim = @animate for i in 1:length(df.s)
           plot(v_pos, f_s_pos[i]; label=@sprintf("f_s at s = %.3e cm", df.s[i]))
           xlabel!("velocity (cm/s)")
           ylabel!("probability density (s/cm)")
           title!("Space-Dependent Maxwellian")
           ylims!(-1e-28, 2e-27)
       end;

gif(anim, "Simple_BGK_Spatial_Discretization_Pos.gif"; fps=20)

````

![Simple_BGK_Spatial_Discretization_Pos](https://global.discourse-cdn.com/julialang/original/3X/b/a/bacfdffefab895a55f7acb108ad8013ff882900f.gif)

Now, I want to try to create this animation again, but instead, we have negative values for velocity (`v_neg` instead of `v_pos`) and we have to step through s values from the max value of `s` to the min value of `s`. Here is my attempt:

```julia
nu = 1 # collisional frequency [s^-1]
v_neg = range(0, -1e10, step=-1e8) # negative velocity [cm/s]

f_s_neg = [Maxwellian.(v_neg, df.T[end])] # because boundary density is high at right footpoint, assume f_N = f_M

for is in reverse(eachindex(df.s)) # for negative velocities, integrate from s_max to 0 
    local_F_s_neg = Maxwellian.(v_neg, df.T[is])
    push!(f_s_neg, f_s_neg[length(df.s)+1-is] - (nu./v_neg).*(f_s_neg[length(df.s)+1-is] - local_F_s_neg)*(df.ds[is]))
end

anim = @animate for i in 1:length(df.s)
           plot(v_neg, f_s_neg[i]; label=@sprintf("f_s_neg at s = %.3e cm", df.s[i]))
           xlabel!("velocity (cm/s)")
           ylabel!("probability density (s/cm)")
           title!("Space-Dependent Maxwellian")
           ylims!(-1e-30, 2e-27)
       end;

gif(anim, "Simple_BGK_Spatial_Discretization_Neg.gif"; fps=20)

```

![Simple_BGK_Spatial_Discretization_Neg](https://global.discourse-cdn.com/julialang/original/3X/b/f/bfdff4997bdc503b82217caa03fd5b7a2784f123.gif)

I was expecting the function to look symmetric to the previous plot along the y axis but this is not the case. This might be a problem specific in my field (physics), but before investigating further, from a Julia programming perspective, are there any glaring issues with the way I set up the code for the negative velocity case that may be the problem?

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [February 20, 2022, 11:00am UTC](https://discourse.julialang.org/t/how-to-make-this-distribution-symmetric/76728/2 "2022-02-20T11:00:51Z")

</div>

> [@CentraCep](#):
>
> ```julia
> ... - (nu./v_neg).*(f_s_neg[length(df.s)+1-is] - local_F_s_neg)*(df.ds[is])
> 
> ```

Since you are going backwards, you probably need to change the `-` sign to `+`:

```julia
... + (nu./v_neg).*(f_s_neg[length(df.s)+1-is] - local_F_s_neg)*(df.ds[is]))
```

---

<div class="post-metadata">

**Author:** ![CentraCep](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/centracep/32/32454_2.png) [@CentraCep](https://discourse.julialang.org/u/CentraCep)\
**Post date:** [February 20, 2022, 11:34pm UTC](https://discourse.julialang.org/t/how-to-make-this-distribution-symmetric/76728/3 "2022-02-20T23:34:36Z")

</div>

@rafael.guerra This worked! Thank you so much!

Do you know how I could rewrite my code to combine the positive anfd negative velocity curves onto one graph/animation?

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [February 20, 2022, 11:40pm UTC](https://discourse.julialang.org/t/how-to-make-this-distribution-symmetric/76728/4 "2022-02-20T23:40:12Z")

</div>

The easiest might be to concatenate the vectors to plot, say:

```julia
plot([v_neg; v_pos], [f_s_neg[i]; f_s_pos[i]]; ...)

```

---

<div class="post-metadata">

**Author:** ![CentraCep](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/centracep/32/32454_2.png) [@CentraCep](https://discourse.julialang.org/u/CentraCep)\
**Post date:** [February 21, 2022, 2:58am UTC](https://discourse.julialang.org/t/how-to-make-this-distribution-symmetric/76728/5 "2022-02-21T02:58:39Z")

</div>

Hello again. So now, I am trying to recreate the plot but now using an implicit numerical method rather than an explicit numerical method.

```julia
# Spatial discretization 

v_pos = range(0, 1e10, step=1e8) # positive velocity [cm/s]
v_neg = range(0, -1e10, step=-1e8) # negative velocity [cm/s]

nu = 1 # constant collisional frequency [s^-1]

f_s_pos_imp = [Maxwellian.(v_pos, df.T[1])] # because boundary density is high at left footpoint, assume f_N = f_M
f_s_neg_imp = [Maxwellian.(v_neg, df.T[end])] # because boundary density is high at right footpoint, assume f_N = f_M

for is in 1:length(df.s) # for positive velocities, integrate from 0 to s_max
    local_f_s_pos = Maxwellian.(v_pos, df.T[is])
    #push!(f_s_pos_exp, f_s_pos_exp[is] - (nu ./v_pos).*(f_s_pos_exp[is] - local_f_s_pos_exp)*(df.ds[is]))   
    push!(f_s_pos_imp, (f_s_pos_imp[is] .+ (nu./v_pos).*local_f_s_pos.*df.ds[is])/(1 .+ (nu./v_pos).*(df.ds[is])))
    
end

for is in reverse(eachindex(df.s)) # for negative velocities, integrate from s_max to 0 
    local_f_s_neg = Maxwellian.(v_neg, df.T[is])
    #push!(f_s_neg_exp, f_s_neg_exp[length(df.s)+1-is] + (nu ./v_neg).*(f_s_neg_exp[length(df.s)+1-is] - local_f_s_neg_exp)
     # *(df.ds[is]))
    push!(f_s_neg_imp, (f_s_neg_imp[is] .+ (nu./v_neg).*local_f_s_neg.*df.ds[is])/(1 .+ (nu./v_neg).*(df.ds[is])))
end

anim = @animate for i in 1:length(df.s)
           plot([v_neg; v_pos], [f_s_neg_imp[i]; f_s_pos_imp[i]]; label=@sprintf("f_s at s = %.3e cm", df.s[i]))
           xlabel!("velocity (cm/s)")
           ylabel!("probability density (s/cm)")
           title!("Space-Dependent Maxwellian")
           ylims!(-1e-28, 2e-27)
       end;

gif(anim, "TEST.gif"; fps=30)

```

I thought I had ensured that all my lines in the loops were performing vectorized operations by adding `.`, but I receive this error:

```julia
MethodError: no method matching Vector{Float64}(::Matrix{Float64})
Closest candidates are:
  Array{T, N}(::AbstractArray{S, N}) where {T, N, S} at C:\Users\Acer\AppData\Local\Programs\Julia-1.7.0\share\julia\base\array.jl:563
  Vector{T}() where T at C:\Users\Acer\AppData\Local\Programs\Julia-1.7.0\share\julia\base\boot.jl:476
  Vector{T}(::SuiteSparse.CHOLMOD.Dense{T}) where T at C:\Users\Acer\AppData\Local\Programs\Julia-1.7.0\share\julia\stdlib\v1.7\SuiteSparse\src\cholmod.jl:850

```

What may I be missing here?

---

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [February 21, 2022, 3:05am UTC](https://discourse.julialang.org/t/how-to-make-this-distribution-symmetric/76728/6 "2022-02-21T03:05:00Z")

</div>

This error is telling you that you can’t make a Vector{Float64} by passing the constructor a Matrix{Float64}

since you didn’t get the full error, I’m not sure what line is causing the issue.

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [February 21, 2022, 3:13pm UTC](https://discourse.julialang.org/t/how-to-make-this-distribution-symmetric/76728/7 "2022-02-21T15:13:50Z")

</div>

> [@CentraCep](#):
>
> `... df.ds[is])/(1 .+ ( ...`

You seem to have forgotten one dot for the division.  
Instead of sprinkling dots all over the place and whenever possible, just use the `@.` macro at the beginning of the expression and leave the rest alone.

**NB:**

- In `nu./v_pos` there is division by 0
- As the result is symmetrical, you just need to model the positive velocities, and then mirror the resulting arrays for display

---

<div class="post-metadata">

**Author:** ![CentraCep](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/centracep/32/32454_2.png) [@CentraCep](https://discourse.julialang.org/u/CentraCep)\
**Post date:** [February 21, 2022, 10:44pm UTC](https://discourse.julialang.org/t/how-to-make-this-distribution-symmetric/76728/8 "2022-02-21T22:44:28Z")

</div>

> [@CentraCep](#):
>
> ```julia
> anim = @animate for i in 1:length(df.s)
> plot(v_pos, f_s_pos[i]; label=@sprintf("f_s at s = %.3e cm", df.s[i]))
> xlabel!("velocity (cm/s)")
> ylabel!("probability density (s/cm)")
> title!("Space-Dependent Maxwellian")
> ylims!(-1e-28, 2e-27)
> end;
> 
> gif(anim, "Simple_BGK_Spatial_Discretization_Pos.gif"; fps=20)
> 
> ```

Ah, thank you very much. The `@. macro` makes things convenient.