Ecuación de continuidad

Trabajamos en dimensión \(d = 1\) y estudiamos leyes de conservación de la forma \[\begin{equation} \tag{C} \label{eq:conservacion} \frac{\partial}{\partial t} \rho = - \diver \flux \end{equation}\]

Conservación

Estos métodos son naturales para leyes de conservación porque \[\begin{align*} \frac{\diff}{\diff t} \sum_{i=1}^N \Big(x_{i+\frac 1 2} - x_{i-\frac 1 2}\Big)\rho_i (t) &= - \left( \flux_{\frac 3 2} (t) - \flux_{\frac 1 2} (t) \right) - \left( \flux_{\frac 5 2} (t) - \flux_{\frac 3 2} (t) \right) - \cdots - \left( \flux_{N-\frac 1 2} (t) - \flux_{N + \frac 1 2} (t) \right) \\ &= - F_{N+\frac 1 2} (t) + F_{\frac 1 2} (t) \end{align*}\]

Ecuación de transporte

El caso más sencillo es \(\flux(t,x) = a \rho (t,x) \ge 0\) con \(a\) constante

Soluciones exactas

Suelen llamarse características a curvas \((t,X(t))\) tales que la EDP se puede resolver aisladamente a lo largo de curva.

Restricción a un intervalo cerrado

Si queremos trabajar un intervalo acotado, digamos \((0,1)\), no puede admitir condiciones de flujo nulo.

Discretización temporal

Hemos de elegir un método de discretizar \(\frac{\diff \rho_i}{\diff t}\).

Probaremos el método de Euler \[ \frac{\rho_i^{n+1} - \rho_i^n}{\Delta t} = - \frac{ F_{i+\frac 1 2}^{n+1} - F_{i - \frac 1 2}^{n+1} }{\Delta x} \]

\(F_{i+\frac 1 2}^{n+1} = a\rho_i^n\)

Tenemos el método \[ \rho^{n+1}_i = \rho^n_i - a\frac{\Delta t}{\Delta x} ( \rho_{i}^n - \rho_{i-1}^n ) \]

ρ = copy(ρ0); Δt = Δx/2
@gif for n=1:ceil(T/Δt)
    ρn = copy(ρ)
    
    ρ[1] = ρn[1] - Δt/Δx * a * (ρn[1] - ρn[end])
    for i = 2:length(ρ)
        ρ[i] = ρn[i] - a * Δt/Δx * ( ρn[i] - ρn[i-1] )
    end
    
    plot(x,ρ,
        ylim=(0.0,1.2),
        title="t=$(round(n*Δt,digits=3))",
        linewidth=2,thickness_scaling = 1.5,label="") 
end

Ejercicio. \(F_{i+\frac 1 2}^{n+1} = a\rho_{i+1}^n\)

\[\rho^{n+1}_i = \rho^n_i - a\frac{\Delta t}{\Delta x} ( \rho_{i+1}^n - \rho_{i}^n )\]

Code

ρ = copy(ρ0); Δt = Δx/2
@gif for n=1:ceil(T/Δt)
    ρn = copy(ρ)
    
    ρ[1:end-1] = ρn[1:end-1] - Δt/Δx * a * ( ρn[2:end] - ρn[1:end-1] )
    ρ[end] = ρn[end] - a * Δt/Δx * ( ρn[1] - ρn[end] )
    
    plot(x,ρ,
        ylim=(0.0,1.2),
        title="t=$(round(n*Δt,digits=3))",
        linewidth=2,thickness_scaling = 1.5,label="") 
end

Up-winding. Euler explícito

Tomamos el método \[ \rho^{n+1}_i = \rho^n_i - a\frac{\Delta t}{\Delta x} ( \rho_{i}^n - \rho_{i-1}^n ) %= (1- a\frac{\Delta t}{\Delta x}) \rho_i^n + a\frac{\Delta t}{\Delta x} \rho_{i-1}^n. \]

Ejercicios

Ejercicio. Comprobar que esta condición es suficiente para si \(\rho_i^0 \ge 0\) para todo \(i\), entonces \(\rho_i^1 \ge 0\) para todo \(i\).

Up-winding. Semi-discretización espacial

Tomemos \(a > 0\). Observamos que \[ \begin{aligned} \frac{d \rho_i}{dt} (t) &= \frac{a}{\Delta x} \rho_{i} (t) \,\, \underbrace{- \frac{a}{\Delta x} \rho_{i+1}(t) }_{negativo} \end{aligned} \]

La solución es \[ \rho_i(t) = e^{\frac{a}{\Delta x}t} \rho_i(0) - \frac{a}{\Delta x} \int_0^t e^{\frac{a}{\Delta x}(t-s)} \rho_{i+1} (s) d s\]

Podría cambiar de signo.

Up-winding. Semi-discretización espacial

Sin embargo \[\begin{aligned} \frac{d \rho_i}{dt} (t) = - \frac{a}{\Delta x} \rho_{i} (t) \,\, \underbrace{+ \frac{a}{\Delta x} \rho_{i-1}(t) }_{positivo} \end{aligned} \]

La solución es \[ \rho_i(t) = e^{-\frac{a}{\Delta x}t} \rho_i(0) + \frac{a}{\Delta x} \int_0^t e^{-\frac{a}{\Delta x}(t-s)} \rho_{i-1} (s) d s \]

Si \(\rho_i (0) \ge 0\) para todo \(i\). Entonces \(\rho_i \ge 0\).

Esta elección de flujo se conoce como up-winding.

Simulación que cumple CFL

Δt = .75*Δx
ρ = copy(ρ0); 

@gif for n=1:ceil(T/Δt)
    ρn = copy(ρ)
    
    ρ[1]     = ρn[1]     - a * Δt/Δx * ( ρn[1]   - ρn[end] )
    ρ[2:end] = ρn[2:end] - a * Δt/Δx * ( ρn[2:end] - ρn[1:end-1] )
    
    plot(x,ρ,
        ylim=(0.0,1.2),
        title="t=$(round(n*Δt,digits=3))",
        linewidth=2,thickness_scaling = 1.5,label="") 
end