³⁰² views

Kernel: Julia

In [9]:

using Plots, FileIO, Interact

In [0]:

Dinamica Browniana

Un coloide es un sistema formado por 2 o más fases. Típicamente una fase sólida y una líquida. La fase líquida se encuentra formada por partículas de unas pocas micras. En este caso, cada partícula sigue un movimiento browniano, pues la partícula es sufcientemente pequeña como para ver afectado su movimiento por la colisión de las moléculas del líquido, pero no tan pequeña que se pueda considerar parte del fluido.

En esta práctica desarrollarás un método para poder estudiar sistemas coloidales. La idea fundamental es resolver la ecuación de Langevin con el método de Euler y después aplicar esto a un conjunto de partículas que puede interactuar entre ellas.

Paso 1: Movimiento Browniano

En 1827 Robert Brown, un biólogo inglés, observó en el microscopio que las partículas de polen en agua se movían de forma "extraña", siguiendo trayectorias asarosas. Sus observaciones sirvieron como prueba de que la materia está constituida por moléculas que se mueven ergódicamente en todas direcciones. En 1905, Einstein fue capaz de hacer un modelo estadístico el cual describía adecuadamente el movimiento Browniano. Con él pudo mostrar que, las partículas Brownianas siguen un comportamiento difusivo.

Hoy en día esto se usa en muchas áreas de la física, como la termodinámica, la dinámica de fluidos, biofísica, y muchos otros.

Como un primer paso, tendrás que ser capaz de modelar una caminata aleatoria (o movimiento browniano).

[1] Dibuja una caminata de borracho en 2D. Para esto, genera pasos de tamaño fijo, en una dirección $\theta$ aleatoria.

In [10]:

function borracho2D(x0,h,n) #donde x0 es el opunto de inicio, h el largo del paso y n el número de pasos.
        
    X = []
    Y = []
        
    push!(X, x0[1])
    push!(Y, x0[2])
    
    x = x0
    
    for i in 1:n
    
        θ =  pi*rand()
        x += [h*cos(θ),h*sin(θ)]
        push!(X, x[1])
        push!(Y, x[2])
        
    end
    
    return hcat(X,Y)
end

Out[10]:

borracho2D (generic function with 1 method)

In [11]:

M1 = borracho2D([0,0],1,100)

Out[11]:

101×2 Array{Real,2}:
  0          0       
  0.261248   0.965272
 -0.41534    1.70163 
 -1.19376    2.32938 
 -1.94284    2.99186 
 -2.48718    3.83072 
 -2.73627    4.7992  
 -3.54187    5.39166 
 -4.41518    5.87883 
 -5.41493    5.901   
 -4.77998    6.67356 
 -3.79066    6.81929 
 -2.79141    6.85812 
  ⋮                  
 -8.23835   53.7196  
 -8.82442   54.5298  
 -7.99878   55.094   
 -8.09639   56.0892  
 -8.81213   56.7876  
 -9.17644   57.7189  
 -8.36468   58.3029  
 -7.42058   58.6325  
 -6.4286    58.7589  
 -5.66853   59.4088  
 -5.83554   60.3947  
 -6.57911   61.0634  

In [12]:

plot(M1[:,1],M1[:,2])

Out[12]:

INFO: Initializing package repository /home/user/.julia/v0.6
INFO: Cloning METADATA from https://github.com/JuliaLang/METADATA.jl

GitError(Code:ERROR, Class:Net, curl error: Failed to connect to github.com port 443: Connection timed out
)

Stacktrace:
 [1] macro expansion at ./libgit2/error.jl:99 [inlined]
 [2] clone(::String, ::String, ::Base.LibGit2.CloneOptions) at ./libgit2/repository.jl:276
 [3] #clone#100(::String, ::Bool, ::Ptr{Void}, ::Nullable{Base.LibGit2.AbstractCredentials}, ::Function, ::String, ::String) at ./libgit2/libgit2.jl:562
 [4] (::Base.LibGit2.#kw##clone)(::Array{Any,1}, ::Base.LibGit2.#clone, ::String, ::String) at ./<missing>:0
 [5] (::Base.Pkg.Dir.##8#10{String,String})() at ./pkg/dir.jl:55
 [6] cd(::Base.Pkg.Dir.##8#10{String,String}, ::String) at ./file.jl:70
 [7] init(::String, ::String) at ./pkg/dir.jl:53
 [8] #cd#1(::Array{Any,1}, ::Function, ::Function, ::String, ::Vararg{String,N} where N) at ./pkg/dir.jl:28
 [9] pickDefaultBackend() at /ext/julia/julia-0.6.2/share/julia/site/v0.6/Plots/src/backends.jl:155
 [10] backend() at /ext/julia/julia-0.6.2/share/julia/site/v0.6/Plots/src/backends.jl:174
 [11] Plots.Plot() at /ext/julia/julia-0.6.2/share/julia/site/v0.6/Plots/src/types.jl:81
 [12] #plot#176(::Array{Any,1}, ::Function, ::Array{Real,1}, ::Vararg{Array{Real,1},N} where N) at /ext/julia/julia-0.6.2/share/julia/site/v0.6/Plots/src/plot.jl:56
 [13] plot(::Array{Real,1}, ::Array{Real,1}, ::Vararg{Array{Real,1},N} where N) at /ext/julia/julia-0.6.2/share/julia/site/v0.6/Plots/src/plot.jl:52
 [14] include_string(::String, ::String) at ./loading.jl:522

[2] Calcula el desplazamiento cuadrático medio de estas caminatas aleatorias. Es decir, el promedio de $\langle (x(t)-x(0))^2 \rangle$ como función del tiempo ( $\langle \cdot \rangle$ significa el promedio sobre muchas caminatas aleatorias).

In [0]:

function promedio_caminata(M)
    p = 0
    n = length(M[:,1])
    
    for i in 1:n
    
        p+=(norm(M[i,:])-norm(M[1,:]))^2
    end
    return p/n
end

In [0]:

promedio_caminata(M1)

[3] Para tiempos grandes, el desplazamiento cuadrático medio va como $D t$ , donde $D$ es el coeficiente de difusión. Calcula este coeficiente.

In [0]:

function cD(M)
    p = promedio_caminata(M)
    t = length(M[:,1])
    return(p/t)
end

In [0]:

cD(M1)

[4] Generaliza esto para $N$ dimensiones.

In [0]:

function borrachoND(x0,h,n,N) # Donde las variables son las mismas que en borracho2D y N es el número de dimensiones.
    
    x = x0
    V = vcat(x)
    
    for i in 1:n
        y = rand(1,N)
        x += y*h/norm(y)
        V = vcat(V,x)
    end
    return V
end

In [0]:

M2 = borrachoND([0 0 0],2,5,3)

In [0]:

plot(M2[:,1],M2[:,2],M2[:,3])

Paso 2: Movimiento Browniano y la ecuación de Langevin

Hacer una caminata aleatroia es relativamente sencillo, pero ¿cómo se vería esa misma caminata si estuviramos bajo la acción de algún potencial? Tendría que haber una probabilidad mayor de moverse en determinada dirección. Calcular esa probabilidad como función del potencial no parece tan fácil.

La forma de resolver este problema, es escribiendo la segunda ley de Newton para nuestra partícula browniana de masa $m$ . Una vez escrita la segunda ley de Newton, se utiliza el teorema de equipartición que dice que por cada grado de libertad de cada partícula, la energía cinética está dada por $\frac{1}{2} k_B T$ para reducir términos ( $T$ es la temperatura y $k_b$ la constante de Boltzmann). Entonces, las fuerza $m \ddot{X}$ que actúan sobre una partícula grande de masa $m$ serán:

la fuerza debido al potencial, es decir $- \nabla U(X)$ .
la fuerza debido a la fricción de la partícula en el fluido, es decir $- \gamma \dot{X}$ y finalmente,
una contribución estocástica que depende de la fricción y la energía cinética de las partículas, es decir, $\sqrt{2 \gamma k_B T} R(t)$ , donde $R(t)$ es una variable aleatoria de tipo gaussiana.

Es decir:

$m\ddot{X} = - \nabla U(X) - \gamma \dot{X} + \sqrt{2 \gamma k_B T} R(t)$

Esta ecuación es complicada y en general, lo que nos interesa simular es el caso donde con respecto a nuestro sistema completo, $m$ es muy pequeña. En el caso donde $m$ tiende a 0, esta equación se reduce a $0 = - \nabla U(X) - \gamma \dot{X} + \sqrt{2 \gamma k_B T} R(t)$ , que es la equación de Lagevin, base para dinámica Browniana.

La dificultad de resolver esta clase de ecuaciones consiste en la $R(t)$ , pues no es una variable aleatoria homogenea, sino gaussiana.

[1] haz una función R(t) que te genere un valor aleatorio entre -Inf e Inf, con una distribución de probabilidad de tipo gaussiano.

In [0]:

function R(t,a,b)
 
    y = a*randn() + b
   
    return y
end

In [0]:

X4 = []
for i in 1:100000
    push!(X4, R(i,1,0))
end

histogram(X4,nbins= 100)

[2] Para el caso donde U(x) = 0, resuelve la ecuación de dinámica Browniana usando el método de Euler (recuerda que R depende del tiempo. Deberías obtener simplmente una partícula Browniana libre.

In [0]:

function Metodo_Euler(f, x0; t0 = 0, tf = 10, h = 0.0001) #los argumentos de la función son la función f, la condición
    #inicial x0, el tiempo inicial y el tiempo final t0 y tf y el tamaño de los pasos h.
    
    n = floor(Int,abs(tf - t0)/h) #definimos el tamaño de los arreglos como un número entero (para depués poder usarlo 
    #para definir índices) dividiendo el intervalo temporal entre el tamaño de los pasos
    
    ts = hcat(linspace(t0, tf, n)) #el arreglo de los valores del tiempo
    
    rango = zeros(n) #el arreglo donde se guardarán los valores de la función
    
    x = x0 + h*f(x0, ts[1]) #el primer valor de x se calcula a partir de la condición inicial x0
    rango[1] = x #y se guarda como el primer elemento del arreglo rango
    
    for i in 2:n #para todos los demás elementos
        x = x + h*f(x, ts[i]) #redefinimos el valor de x[n+1] a partir del valor de x[n]
        rango[i] = x #y los guardamos en el arreglo
    end
    
    return hcat(ts, rango) #la función devuelve un arreglo 2xn cuya primera columna son los valores de t y la segunda
    #los valores del rango
end

In [0]:

k = 1
T = 1
γ = 2
f(x,t) = sqrt(2*k*T/γ)*R(t,1,0)
M3 = Metodo_Euler(f,1)

In [0]:

plot(M3[:,1],M3[:,2])

[3] Calcula el coeficiente de difusión $D$ del paso anterior. Ahora tu coeficiente de difusión depende de la temperatura y de $\gamma$ . Para $\gamma$ fijo, ¿Cómo depende de la temperatura? ¿Qué relación hay con la $D$ del paso anterior?

In [0]:

cD(M3)

In [0]:

k = 1
T = 1
γ = 1
f(x,t) = sqrt(2*k*T/γ)*R(t,1,0)
M4 = Metodo_Euler(f,1)
@show cD(M4)
plot(M4[:,1],M4[:,2])

Paso 3: Dinámica Browniana.

En este punto ya eres casi un experto en movimiento Browniano. Ahora vamos a pasar a implementar potenciales. Para esto tendrás que usar una versión n-dimensional del método de Euler.

[1] Haz una función que calcule $\nabla U(X)$ , dada cualquier función $U$

In [0]:

function derivada_parcial(f, X, i; h = 1e-10)
    
    H = zeros(length(X)) #definimos un vector H cuyas entradas son todas iguales a 0
    
    H[i] = h*(H[i] + 1) #a la i-ésima entrada de H la convertimos de 0 a h

    der = (f(X + H) - f(X))/h #aplicamos la definición de derivada direccional en la dirección canónica j
    
    return der
end

In [0]:

function fuerza(U, X)
    n = length(X)
    F = []
    
    for i in 1:n
        push!(F, derivada_parcial(U, X, i))
    end
    
    return F
    
end

In [0]:

U(X) = 1/2*(X[1]^2 + X[2]^2)

In [0]:

fuerza(U, [0,1])

[2] Haz una función que aplique el método de Euler n-dimensional para obtener las soluciones de la ecuación de la dinámica Browniana, dado un potencial $U$ , una temperatura $T$ y un coeficiente de fricción $\gamma$ .

0 = - \nabla U(X) - \gamma \dot{X} + \sqrt{2 \gamma k_B T} R(t)

\dot{X} = {1 \over \gamma}\left(- \nabla U(X) + \sqrt{2 \gamma k_B T} R(t) \right)

Constante de Boltzmann: $k_B = 1.38 \times10^{-23}{J \over K}$

Coeficiente de fricción: $\gamma = 6\pi\eta a$

Donde $\eta$ es la viscosidad dinámica del medio y $a$ es el radio de la partícula browniana

In [0]:

function Euler_Vectorial(f, X0; h= 0.001, t0 = 0, tf = 10)
    
    n = floor(Int,abs(tf - t0)/h) #el número de veces que se repetirá la iteración
    
    dim = length(X0) #la dimensión del espacio sobre el cual está definido la función vectorial, la cual nos indica
    #cuántos rangos (uno por entrada de la función) necesitaremos 
    
    ts = hcat(linspace(t0, tf, n)) #hacemos la lista de los tiempos en los cuales aplicaremos el método de Euler
    
    rango = zeros(n, dim) #hacemos un arreglo donde guardaremos los resultados de los rangos, debe de tener dimension nxdim
    
    X = zeros(dim) #en el arreglo X se guardan las coordenadas modificadas del vector de posiciones
    
    for i in 1:dim
        X[i] = X0[i] #al tiempo t0, la función coordenada f[i] vale la condición inicial X0[i]
        rango[1,i] = X[i] 
    end
    
    for j in 2:n
        for i in 1:dim
            X[i] = X[i] + h*f(X, ts[j])[i] #esta sería la función Paso_Euler pero es lo suficientemente simple como
            #para que no valga la pena hacer una función separada para ella.
            rango[j,i] = X[i]
        end
    end
    return hcat(ts,rango) #guarda los resultados en un arreglo cuya primera columna sean los valores de ts y en las 
    #siguientes los resultados correspondientes de los rangos de las funciones coordenadas
end

In [0]:

function RN(t,a,b, N) #la misma función de valores aleatorios pero que ahora regrese un vector
    
    Y = []
    for i in 1:N
        push!(Y, a*randn() + b)
    end
   
    return Y
end

In [0]:

function dinamica_browniana(U, T, γ, X0; kb = 1.38e-23)
    
    dim = length(X0)
    
    f(X, t) = (1/γ)*(-fuerza(U, X) + sqrt(2*kb*γ*T)*RN(t, 1, 0, dim))
    
    Dim_Brown = Euler_Vectorial(f, X0)
    
    return Dim_Brown
end

[3] Haz una simulación (una animación), de una partícula de polen en agua bajo el cámpo gravitatorio de la tierra. En este punto suele haber problemas. A la hora de graficar, observa la escala y asegúrate que está realmente cayendo la partícula de Polen para campos suficientemente grandes (que la caída no sea de $10^{-n}$ con $n \sim 15$ .

$m_{polen} = 247 \times 10^{-12} kg$

$T = 300 K$

$U(x,y) = mgy$

In [0]:

ci = 1 #condición inicial
m = 247e-12
T = 300
η = 0.8701 #viscosidad dinámica del agua a 300 K
g = 9.81
a = 7e-6 #radio aproximado de una partícula de polen
γ = 6π*η*a
U(X) = m*g*X[2]

DB = dinamica_browniana(U, T, γ, [0,ci])

plt1 = plot(DB[:,2],DB[:,3], label = "Trayectoria de la particula de polen en agua")
display(plt1)

n = 100
@gif for i in 1:n
    
   scatter([0],[ci],label="origen")
    scatter!([DB[i,2]],[DB[i,3]], xlims = (minimum(DB[1:n,2]),maximum(DB[1:n,2])), ylims = (minimum(DB[1:n,3]),maximum(DB[1:n,3])), label = "particula de polen")
end

[4] Haz una simulación de una partícula browniana cargada negativamente, en un potencial de Coulomb de otra partícula (masiva) cargada positivamente. Tendrás problemas si las partículas se acercan demasiado, piensa como resolver ese caso de forma realista.

El potencial de una partícula de carga con carga $e$ en presencia de un potencial generado por una partícula de carga $q$ ubicada en el origen es:

U(X) = \frac{1}{4 \pi \epsilon_0} \frac{eq}{||X||}

$\epsilon_0 = 8.85\times 10^{-12} F/m$

$e = 1.602\times 10^{-19} C$ (carga de un electrón)

$a = 4.74\times 10^{-4}m$ (radio de la partícula cargada)

In [0]:

ci = 1e-4
cix = 1e-4
ϵ0 = 8.83e-12
e = -1.602e-19
q = 1.602e-10 #suponemos que la partícula "masiva" que genera el potencial tiene una carga 9 órdenes de
                #magnitud mayor que la otra partícula
R1 = 1e-6
T = 300
a = 4.74e-4
#η = 1.846e-3 viscosidad dinámica del aire a 300 K
γ = 5.5

U_out(X) = (e*q)/(4π*ϵ0*norm(X))#*exp(-norm(X)/9e-6)
U_in(X) = (e*q)/(4π*ϵ0*R1) 


DB_out = dinamica_browniana(U_out, T, γ, [cix,ci])#, kb = 1) potencial para norm(X)>R

DB_in = dinamica_browniana(U_in, T, γ, [cix,ci])#potencial para norm(x)<R

DB = []

for i in 1:length(DB_out[:,2])
   
    if norm(DB_out[i,:]) > R1 #para el listado se verifica si la norm(X)<R y se sustituyen los valores la de ecuación diferencial
        DB_in[i,2] = DB_out[i,2]
        DB_in[i,3] = DB_out[i,3]
    end
    DB = DB_in
end
plt1 = plot(DB[:,2],DB[:,3], label = "Trayectoria de la particula cargada")
display(plt1)

rang = 1:500

@gif for i in rang
    
    scatter([0],[0], label = "particula cargada positivamente")
    
    scatter!([DB[i,2]],[DB[i,3]], xlims = (minimum(DB[rang,2]),maximum(DB[rang,2])), ylims = (minimum(DB[rang,3]),maximum(DB[rang,3])), label = "particula cargada negativamente")
    #scatter!([cix],[ci], label = "posicion inicial")
end

Paso 4: Dinámica Browniana con muchas partículas.

Hasta ahora sólo hemos visto qué pasa cuando tenemos 1 partícula. Pero puede ser que las partículas interactuen entre ellas. ¿Cómo podemos hacer esto? El potenciál U dependerá de las posiciones de las demás partículas.

[1] Haz una simulación considerando 10 partículas (posiciones iniciales aleatorias). Las partículas deben tener algún potencial de esfera suave (el que gustes). ¿Qué sucede?

Suponiendo un conjunto de 10 particulas con carga negativa en presencia de una particula mas masiva de carga positiva tendremos que incluir en el potencial las interacciones entre las particulas negativas y las interacciones entre las particulas positivas, e,i $U = U_{Ne} + U_{ee}$ donde $U_{Ne} = \frac{1}{4πϵ_0}\sum_{i}\frac{q e_i}{r_i}$ y $U_{ee} = \frac{1}{4πϵ_0}\sum_{i \neq j}\frac{e_i e_j}{r_{ij}}$

In [0]:

1e-4 *rand(10,2)[1,:]

In [0]:


X0 = 1e-4 *rand(10,2)
ϵ0 = 8.83e-12
e = -1.602e-19
q = 1.602e-10 #suponemos que la partícula "masiva" que genera el potencial tiene una carga 9 órdenes de
                #magnitud mayor que la otra partícula

T = 300
a = 4.74e-4
#η = 1.846e-3 viscosidad dinámica del aire a 300 K
γ = 5.5
U1([x1 ,x2,x3,x4,x5,x6,x7,x8,x9,x10]) = (e*q)/(4π*ϵ0*abs(x1))+ (e^2 /4π*ϵ0)*(1/abs(x1 - x2) + 1/abs(x1 - x3) +1/abs(x1- x4) +1/norm(x1 - x5) +1/abs(x1 - x6)+ 1/abs(x1 - x7)+ 1/abs(x1 - x8)+ 1/abs(x1 - x9) + 1/abs(x1 - x10))
U2(X2) = (e*q)/(4π*ϵ0*norm(X))+ (e^2 /4π*ϵ0)*(1/norm(X - X0[1,:]) + 1/norm(X - X0[3,:]) +1/norm(X - X0[4,:]) +1/norm(X - X0[5,:]) +1/norm(X - X0[6,:])+ 1/norm(X - X0[7,:])+ 1/norm(X - X0[8,:])+ 1/norm(X - X0[9,:]))
U3(X3) = (e*q)/(4π*ϵ0*norm(X))+ (e^2 /4π*ϵ0)*(1/norm(X - X0[1,:]) + 1/norm(X - X0[2,:]) +1/norm(X - X0[4,:]) +1/norm(X - X0[5,:]) +1/norm(X - X0[6,:])+ 1/norm(X - X0[7,:])+ 1/norm(X - X0[8,:])+ 1/norm(X - X0[9,:]))
U4(X4) = (e*q)/(4π*ϵ0*norm(X))+ (e^2 /4π*ϵ0)*(1/norm(X - X0[1,:]) + 1/norm(X - X0[2,:]) +1/norm(X - X0[3,:]) +1/norm(X - X0[5,:]) +1/norm(X - X0[6,:])+ 1/norm(X - X0[7,:])+ 1/norm(X - X0[8,:])+ 1/norm(X - X0[9,:]))
U5(X5) = (e*q)/(4π*ϵ0*norm(X))+ (e^2 /4π*ϵ0)*(1/norm(X - X0[1,:]) + 1/norm(X - X0[2,:]) +1/norm(X - X0[3,:]) +1/norm(X - X0[4,:]) +1/norm(X - X0[6,:])+ 1/norm(X - X0[7,:])+ 1/norm(X - X0[8,:])+ 1/norm(X - X0[9,:]))
U6(X6) = (e*q)/(4π*ϵ0*norm(X))+ (e^2 /4π*ϵ0)*(1/norm(X - X0[1,:]) + 1/norm(X - X0[2,:]) +1/norm(X - X0[3,:]) +1/norm(X - X0[4,:]) +1/norm(X - X0[5,:])+ 1/norm(X - X0[7,:])+ 1/norm(X - X0[8,:])+ 1/norm(X - X0[9,:]))
U7(X7) = (e*q)/(4π*ϵ0*norm(X))+ (e^2 /4π*ϵ0)*(1/norm(X - X0[1,:]) + 1/norm(X - X0[2,:]) +1/norm(X - X0[3,:]) +1/norm(X - X0[4,:]) +1/norm(X - X0[5,:])+ 1/norm(X - X0[6,:])+ 1/norm(X - X0[8,:])+ 1/norm(X - X0[9,:]))
U8(X8) = (e*q)/(4π*ϵ0*norm(X))+ (e^2 /4π*ϵ0)*(1/norm(X - X0[1,:]) + 1/norm(X - X0[2,:]) +1/norm(X - X0[3,:]) +1/norm(X - X0[4,:]) +1/norm(X - X0[5,:])+ 1/norm(X - X0[6,:])+ 1/norm(X - X0[7,:])+ 1/norm(X - X0[8,:]))

#V = []
#p(X) =  0
#for i in 1:10
 #   for j in 1:10
  #      if !(i == j) 
   #         p += (e^2)/(4π*ϵ0*norm(X - X0[j,:]))
    #    end
   # end
   # U(X) = U(X) + p(X) 
   # DB = dinamica_browniana(U, T, γ, X0[i,:])
   # V = vcat(V,DB)
   # U(X) = (e*q)/(4π*ϵ0*norm(X))
#end
    
DB1 = dinamica_browniana(U1, T, γ, X0[1,:])
DB2 = dinamica_browniana(U2, T, γ, X0[2,:])
DB3 = dinamica_browniana(U3, T, γ, X0[3,:])
DB4 = dinamica_browniana(U4, T, γ, X0[4,:])
DB5 = dinamica_browniana(U5, T, γ, X0[5,:])
DB6 = dinamica_browniana(U6, T, γ, X0[6,:])
DB7 = dinamica_browniana(U7, T, γ, X0[7,:])
DB8 = dinamica_browniana(U8, T, γ, X0[8,:])


plt1 = plot(DB1[:,2],DB1[:,3], label = "Trayectoria de la particula cargada")
display(plt1)

rang = 1:500

plt2 = plot(DB2[:,2],DB2[:,3], label = "Trayectoria de la particula cargada")
display(plt2)

rang = 1:500

plt3 = plot(DB3[:,2],DB3[:,3], label = "Trayectoria de la particula cargada")
display(plt3)

rang = 1:500

plt4 = plot(DB4[:,2],DB4[:,3], label = "Trayectoria de la particula cargada")
display(plt4)

rang = 1:500

plt5 = plot(DB5[:,2],DB5[:,3], label = "Trayectoria de la particula cargada")
display(plt5)

rang = 1:500

plt6 = plot(DB6[:,2],DB6[:,3], label = "Trayectoria de la particula cargada")
display(plt6)

rang = 1:500

plt7 = plot(DB7[:,2],DB7[:,3], label = "Trayectoria de la particula cargada")
display(plt7)

rang = 1:500

plt8 = plot(DB8[:,2],DB8[:,3], label = "Trayectoria de la particula cargada")
display(plt8)

rang = 1:500


@gif for i in rang
    
    scatter([0],[0], label = "particula cargada positivamente")
    
    scatter!([DB[i,2]],[DB[i,3]], xlims = (minimum(DB[rang,2]),maximum(DB[rang,2])), ylims = (minimum(DB[rang,3]),maximum(DB[rang,3])), label = "particula cargada negativamente")
    #scatter!([cix],[ci], label = "posicion inicial")
end

In [0]:

U(X,Y) = (e^2 /4π*ϵ0)*(1/norm(X - Y)
function U(i,j)
        U = []
    for i in 1:10
            if i != j
                U += U(X, Xj)
                
end

[2] Ahora cambia de esfera suave a partículas cargadas, 5 positivas, 5 negativas.

In [0]:

X0 = 1e-4 *rand(10,2)
ϵ0 = 8.83e-12
e = -1.602e-19
q = 1.602e-10 #suponemos que la partícula "masiva" que genera el potencial tiene una carga 9 órdenes de
                #magnitud mayor que la otra partícula

T = 300
a = 4.74e-4
#η = 1.846e-3 viscosidad dinámica del aire a 300 K
γ = 5.5
U1(X) = (e*q)/(4π*ϵ0*norm(X))+ (e^2 /4π*ϵ0)*(1/norm(X - X0[2,:]) + 1/norm(X - X0[3,:]) +1/norm(X - X0[4,:]) +1/norm(X - X0[5,:]) -(1/norm(X - X0[6,:])+ 1/norm(X - X0[7,:])+ 1/norm(X - X0[8,:])+ 1/norm(X - X0[9,:])))
U2(X) = (e*q)/(4π*ϵ0*norm(X))+ (e^2 /4π*ϵ0)*(1/norm(X - X0[1,:]) + 1/norm(X - X0[3,:]) +1/norm(X - X0[4,:]) +1/norm(X - X0[5,:]) -(1/norm(X - X0[6,:])+ 1/norm(X - X0[7,:])+ 1/norm(X - X0[8,:])+ 1/norm(X - X0[9,:])))
U3(X) = (e*q)/(4π*ϵ0*norm(X))+ (e^2 /4π*ϵ0)*(1/norm(X - X0[1,:]) + 1/norm(X - X0[2,:]) +1/norm(X - X0[4,:]) +1/norm(X - X0[5,:]) -(1/norm(X - X0[6,:])+ 1/norm(X - X0[7,:])+ 1/norm(X - X0[8,:])+ 1/norm(X - X0[9,:])))
U4(X) = (e*q)/(4π*ϵ0*norm(X))+ (e^2 /4π*ϵ0)*(1/norm(X - X0[1,:]) + 1/norm(X - X0[2,:]) +1/norm(X - X0[3,:]) +1/norm(X - X0[5,:]) -(1/norm(X - X0[6,:])+ 1/norm(X - X0[7,:])+ 1/norm(X - X0[8,:])+ 1/norm(X - X0[9,:])))
U5(X) = (e*q)/(4π*ϵ0*norm(X))+ (e^2 /4π*ϵ0)*(1/norm(X - X0[1,:]) + 1/norm(X - X0[2,:]) +1/norm(X - X0[3,:]) +1/norm(X - X0[4,:]) -(1/norm(X - X0[6,:])+ 1/norm(X - X0[7,:])+ 1/norm(X - X0[8,:])+ 1/norm(X - X0[9,:])))
U6(X) = (e*q)/(4π*ϵ0*norm(X))+ (e^2 /4π*ϵ0)*(1/norm(X - X0[1,:]) + 1/norm(X - X0[2,:]) +1/norm(X - X0[3,:]) +1/norm(X - X0[4,:]) -(1/norm(X - X0[5,:])+ 1/norm(X - X0[7,:])+ 1/norm(X - X0[8,:])+ 1/norm(X - X0[9,:])))
U7(X) = (e*q)/(4π*ϵ0*norm(X))+ (e^2 /4π*ϵ0)*(1/norm(X - X0[1,:]) + 1/norm(X - X0[2,:]) +1/norm(X - X0[3,:]) +1/norm(X - X0[4,:]) -(1/norm(X - X0[5,:])+ 1/norm(X - X0[6,:])+ 1/norm(X - X0[8,:])+ 1/norm(X - X0[9,:])))
U8(X) = (e*q)/(4π*ϵ0*norm(X))+ (e^2 /4π*ϵ0)*(1/norm(X - X0[1,:]) + 1/norm(X - X0[2,:]) +1/norm(X - X0[3,:]) +1/norm(X - X0[4,:]) -(1/norm(X - X0[5,:])+ 1/norm(X - X0[6,:])+ 1/norm(X - X0[7,:])+ 1/norm(X - X0[8,:])))

#V = []
#p(X) =  0
#for i in 1:10
 #   for j in 1:10
  #      if !(i == j) 
   #         p += (e^2)/(4π*ϵ0*norm(X - X0[j,:]))
    #    end
   # end
   # U(X) = U(X) + p(X) 
   # DB = dinamica_browniana(U, T, γ, X0[i,:])
   # V = vcat(V,DB)
   # U(X) = (e*q)/(4π*ϵ0*norm(X))
#end
    
DB1 = dinamica_browniana(U1, T, γ, X0[1,:])
DB2 = dinamica_browniana(U2, T, γ, X0[2,:])
DB3 = dinamica_browniana(U3, T, γ, X0[3,:])
DB4 = dinamica_browniana(U4, T, γ, X0[4,:])
DB5 = dinamica_browniana(U5, T, γ, X0[5,:])
DB6 = dinamica_browniana(U6, T, γ, X0[6,:])
DB7 = dinamica_browniana(U7, T, γ, X0[7,:])
DB8 = dinamica_browniana(U8, T, γ, X0[8,:])


plt1 = plot(DB1[:,2],DB1[:,3], label = "Trayectoria de la particula cargada")
display(plt1)

rang = 1:500

plt2 = plot(DB2[:,2],DB2[:,3], label = "Trayectoria de la particula cargada")
display(plt2)

rang = 1:500

plt3 = plot(DB3[:,2],DB3[:,3], label = "Trayectoria de la particula cargada")
display(plt3)

rang = 1:500

plt4 = plot(DB4[:,2],DB4[:,3], label = "Trayectoria de la particula cargada")
display(plt4)

rang = 1:500

plt5 = plot(DB5[:,2],DB5[:,3], label = "Trayectoria de la particula cargada")
display(plt5)

rang = 1:500

plt6 = plot(DB6[:,2],DB6[:,3], label = "Trayectoria de la particula cargada")
display(plt6)

rang = 1:500

plt7 = plot(DB7[:,2],DB7[:,3], label = "Trayectoria de la particula cargada")
display(plt7)

rang = 1:500

plt8 = plot(DB8[:,2],DB8[:,3], label = "Trayectoria de la particula cargada")
display(plt8)

rang = 1:500


@gif for i in rang
    
    scatter([0],[0], label = "particula cargada positivamente")
    
    scatter!([DB[i,2]],[DB[i,3]], xlims = (minimum(DB[rang,2]),maximum(DB[rang,2])), ylims = (minimum(DB[rang,3]),maximum(DB[rang,3])), label = "particula cargada negativamente")
    #scatter!([cix],[ci], label = "posicion inicial")
end

Paso 5: Condiciones a la frontera.

En algunos casos, si las partículas no están confunadas, simplemente se escaparán del sistema. Una forma de confinar las partículas es poniendo un potencial de amarre para todas, pero eso es poco realista. Una mejor forma de confinar el sistema, es poniendo condiciones a la frontera. Un posible caso, es poner una barrera de potencial "suave" (lo pongo entre comillas, porque típicamente se busca que ese potencial se comporte como un potencial duro). Otra opción es considerar al sistema infinito, pero periódico, o lo que es lo mismo, poner condiciones periódicas (tipo pac-man) a la frontera.

[1] Haz las mismas simulaciones que en el paso 4, pero con condiciones periódicas a la frontera. ¿Qué sucede ahora? El arreglo final es igual?

[2] Repite las simulaciones del paso 4, pero esta vez poniendo barreras de potencial "suave" en la frontera, de tal forma quelas partículas estén confinadas. ¿Qué sucede esta vez?

[3] Confinamos las partículas coloidales entre 2 planos. Esto nos deja 2 direcciones donde las partículas se pueden mover libremente. En esas direcciones aplicamos condiciones periódicas a la frontera. pon los planos perpendiculares al eje $z$ . Utiliza un potencial de esfera suave y potencial gravitatorio. ¿Qué sucede?

Paso 5: Mediciones

Ya tenemos nuestras simulaciones, ahora, nos interesa hacer mediciones. Hay muchas posibilidades de medición aquí, porque hay muchas variables que podríamos considerar y muchos modelos complejos que podríamos hacer. Aquí te propondré uno relativamente sencillo. Utilizaremos la geometría del punto 3 del paso 5 y el potencial de esfera suabe de la forma: $U(x_1, x_2) = -C(r- |x_1 -x_2|)/(x_1 -x_2)^2$ con $C$ una constante dada, si $x_1 -x_2 < r$ $r$ el radio de las esferas y $U(x_1,x_2)=0$ si $x_1 -x_2>r$ .

[1] Mide el $\langle X^2 (t) \rangle$ para cada una de las simulaciones. ¿Cómo cambia la difusión (desplazamiento cuadrático medio) como función de $C$ y como función de $r$ ? ¿Como cambia si se "apaga" la gravedad?

[2] Incrementa la densidad del sistema y vuelve a medir $\langle X^2 (t) \rangle$ . ¿cómo cambia la difusión, como función de la densidad? (Para cambiar la densidad puedes incrementar las partículas o reducir el volumen).

In [0]:

ci = 1e-12
cix = 1e-12
ϵ0 = 8.83e-12
e = -1.602e-19
q = 1.602e-19 #suponemos que la partícula "masiva" que genera el potencial tiene una carga 9 órdenes de
                #magnitud mayor que la otra partícula
R1 = 1e-12
T = 300
a = 4.74e-4
#η = 1.846e-3 viscosidad dinámica del aire a 300 K
γ = 5.5

U_out(X) = (e*q)/(4π*ϵ0*norm(X))#*exp(-norm(X)/9e-6)
U_in(X) = (e*q)/(4π*ϵ0*R1) 


DB_out = dinamica_browniana(U_out, T, γ, [cix,ci])#, kb = 1) potencial para norm(X)>R

DB_in = dinamica_browniana(U_in, T, γ, [cix,ci])#potencial para norm(x)<R

DB = []

for i in 1:length(DB_out[:,2])
   
    if norm(DB_out[i,:]) < R1 #para el listado se verifica si la norm(X)<R y se sustituyen los valores la de ecuación diferencial
        DB_in[i,2] = DB_out[i,2]
        DB_in[i,3] = DB_out[i,3]
    end
    DB = DB_in*1e12
end
plt1 = plot(DB[:,2],DB[:,3], label = "Trayectoria de la particula cargada")
        scatter!([0],[0], label = "particula cargada positivamente")
display(plt1)

savefig("plt1.png") 

rang = 1:100

pot_elec = @animate for i in rang
    
    scatter([0],[0], label = "particula cargada positivamente")
    
    scatter!([DB[i,2]],[DB[i,3]], xlims = (minimum(DB[rang,2]),maximum(DB[rang,2])), ylims = (minimum(DB[rang,3]),maximum(DB[rang,3])), label = "particula cargada negativamente")
    #scatter!([cix],[ci], label = "posicion inicial")
end
gif(pot_elec, "pot_elec.gif", fps = 15)

In [0]:

#potencial de oscilador armónico 
U(X,Y) = 1/2*norm(X - Y)^2

U1(X1, X2, X3) = U(X1, X2) + U(X1, X3)
U2(X1, X2, X3) = U(X2, X1) + U(X2, X3)
U3(X1, X2, X3) = U(X3, X1) + U(X3, X2)

P = [U1, U2, U3]
P[1]([1,2],[7,2],[3,2])

In [0]:

fuerza(P,[2,4,5])

In [0]:

In [0]: