Soluciones exactas. Shocks

Si intentamos buscar soluciones exactas \(\rho_t (X_t ) = \rho_0 (y)\) deducimos que \[ \frac{dX_t}{dt} (t) = \rho_0(y) \]

y = copy(x)

@gif for t=0:0.01:2
    plot_ρ( y + ρ0(y)*t , ρ0(y), t )
end

Elección de \(v_{i+\frac 1 2}\)

Escribimos método de Euler explícito para el tiempo y up-winding tomando como aproximación de la velocidad \[v_{i+\frac 1 2} = \frac 1 2 \left( \alpha \rho_i + (1-\alpha)\rho_{i+1}\right)\] para algún \(\alpha \in (0,1)\).

Nótese que \(v_{i+\frac 1 2} \ge 0\) si \(\rho_i \ge 0\). De manera que el up-winding se simplifica.

Condición CFL

Ejercicio. Calcular una condición CFL heredada del caso de velocidad constante:

El método de Euler explícito con \(\alpha = \frac 1 2\) resulta \[ \rho_{i}^{n+1} = \rho_i^n - \frac{\Delta t}{4 \Delta x} \left( {\rho_i^n} \left( \rho_i^n + \rho_{i+1}^n \right) - {\rho_{i-1}^n} \left( \rho_{i-1}^n + \rho_{i}^n \right) \right) \]

Code

T = 2
Δt = Δx;
ρ = ρ0(x);

F = zeros(N,1); #F[i] = F_{i+\frac 1 2}
@gif for n=0:ceil(T/Δt)
    plot_ρ(x, ρ, n*Δt) 

    F[1:end-1]= ρ[1:end-1].* 
                (ρ[1:end-1] + ρ[2:end])/4
    F[end]    = ρ[end]*(ρ[end] + ρ[1])/4
    
    ρ[1]     -= Δt/Δx*( F[1] - F[end] ) 
    ρ[2:end] -= Δt/Δx*( F[2:end] - F[1:end-1] ) 
end

Monotonía

Nótese que \(\overline \rho_i^n = 1\) es solución del método.

Burgers

Si \(v_{i+\frac 1 2}\) es combinación de \(\rho_i\) y \(\rho_{i+1}\) teníamos que \[ \rho_{i}^{n+1} = \rho_i^n - \frac{\Delta t}{\Delta x} \left( {\rho_i^n}\dfrac{\alpha \rho_i^n + (1-\alpha)\rho_{i+1}^n}{2} - {\rho_{i-1}^n}\dfrac{\alpha \rho_{i-1}^n + (1-\alpha)\rho_{i}^n}{2} \right) \]

Experimento numérico

Concluimos que el método monótono es \[ \rho_{i}^{n+1} = \rho_i^n - \frac{\Delta t}{\Delta x} \left( \frac{(\rho_i^n)^2}{2} - \frac{(\rho_{i-1}^n)^2}{2} \right) \] Y añadimos la condición CFL \(\Delta t \le \frac{\Delta x}{\max_i \rho_i^0}\).

Ejercicio. Programar el método anterior añadiendo condiciones de contorno periódicas

Code

Δt = Δx; # Condición de monotonía
ρ = ρ0(x);
@gif for n=0:ceil(T/Δt)  
    global ρ
    plot_ρ(x, ρ, n*Δt) 
    
    F = ρ.^2/2 # puesto que v > 0.
    
    ρ[1]   -= Δt/Δx * ( F[1] - F[end] ) # Condición periódica
    ρ[2:end] -= Δt/Δx * ( F[2:end] - F[1:end-1] ) 
end

Monotonía: más allá del up-winding

En general, se pueden estudiar leyes de conservación de la forma \[ \frac{\partial \rho}{\partial t} + \frac{\partial}{\partial x} \Phi(\rho) = 0, \] es decir \(F = \Phi (\rho)\).

No siempre es necesario verlas como problemas de transporte donde \(F = \rho \frac{\Phi(\rho)}{\rho}\)

(es decir \(v = \frac{\Phi(\rho)}{\rho}\)).

Se habla de estas elecciones según el número de punto. Típicamente la fórmula de \(g\) no cambia. Por ejemplo, un método de 3 puntos resulta \[ \rho_i^{n+1} = \rho_i^n - \frac{\Delta t}{\Delta x} \Big( g(\rho_{i+1}^{n} ,\rho_i^n ) - g(\rho_{i}^{n} ,\rho_{i-1}^n ) \Big) \]

Algunos ejemplos clásicos

Ver [Eymard, Gallouët & Herbin, 2000]

• Esquema de flux-splitting. Si podemos escribir \(f = f_1 + f_2\) con \(f_1' \ge 0 \ge f_2'\) entonces \[ g(a,b) = f_1(a) + f_2(b). \] Entendiendo el problema como una ley de transporte, es decir \(F = \rho v\) y \(v = \frac{ \Phi(\rho) }{\rho}\), se observa que el up-winding corresponde a la descomposición \(f_1(s) = (\Phi(s))_+\) y \(f_2 (s) = (\Phi(s))_-\).

• Esquema de Godunov. Se toma \[ g(a,b) = \begin{dcases} \min_{[a,b]} f & \text{si } a \le b , \\ \max_{[a,b]} f & \text{si } a > b. \end{dcases} \]

• Esquema de Lax-Friedrichs modificado: se toma \[ g(a,b) = \frac{f(a) + f(b)}{2} + D(a-b) \] donde \(2 D \ge |f'(s)|\).

Método de Euler implícito

La versión implícita del método anterior es \[ \rho_{i}^{n+1} = \rho_i^n - \frac{\Delta t}{\Delta x} \left( \frac{(\rho_i^{n+1})^2}{2} - \frac{(\rho_{i-1}^{n+1})^2}{2} \right) \]

Derivadas de \(h\)

Nosotros queremos resolver \[ h(\xi) = 0 \qquad \text{ donde } h_i (\xi) = -\rho_i^n + \xi_{i} + \frac{\Delta t}{\Delta x} \left( \frac{(\xi_i)^2}{2} - \frac{(\xi_{i-1})^2}{2} \right) \]

Ejercicio. Programar el método de Euler implícito.

Escribimos las derivadas de h

N = length(x)
D = zeros(N,N)
D[1,N] = -Δt/Δx
for i=1:N D[i,i]   = +Δt/Δx end
for i=2:N D[i,i-1] = -Δt/Δx end

using LinearAlgebra
Id = Matrix(1.0I, N, N)

 h(ρ,ρn) = ρ  + D * ρ.^2/2 - ρn 
Dh(ρ,ρn) = Id + D * Diagonal(ρ);

Y lanzamos el método

Code

ρ = ρ0(x);
@gif for n=0:ceil(T/Δt) 
    global ρ
    plot_ρ(x, ρ, n*Δt) 
    ρ = fixed_point(ρ,h,Dh)
end

Método Implícito

Otra opción es hacer el método implícito en la difusión \[ \frac{\rho_{i}^{n+1} - \rho_i^n}{\Delta t} + \frac{ \frac{(\rho_i^n)^2}{2} - \frac{(\rho_{i-1}^n)^2}{2}}{\Delta x} = \varepsilon \frac{ \rho_{i+1}^{n+1} - 2 \rho_i^{n+1} + \rho_{i-1}^{n+1} }{(\Delta x)^2} \]

Code

ϵ  = 5e-2
N  = 128 
T = 1

x  = range(-1.0,1.0,N+1); x = vec(x[1:end-1])
Δx = x[2] - x[1]

Δt = Δx ; 

N = length(x)
D = zeros(N,N)
D[1,N] = -Δt/Δx
for i=1:N D[i,i]   = +Δt/Δx end
for i=2:N D[i,i-1] = -Δt/Δx end

D2 = zeros(N,N)
D2[1,N] = -1
for i=1:N D2[i,i]   = +2 end
for i=2:N D2[i,i-1] = -1 end
for i=1:N-1 D2[i,i+1] = -1 end
D2[N,1] = -1

using LinearAlgebra
Id = Matrix(1.0I, N, N)

A = Id + ϵ*Δt/((Δx)^2)*D2

ρ = ρ0(x);
@gif for n=0:ceil(T/Δt) 
    global ρ

    plot_ρ(x, ρ, n*Δt) 
    title!(latexstring(@sprintf("\\varepsilon =%0.3f, \\quad t = %0.3f", ϵ, (n-1)*Δt)))

    F = ρ.^2/2 # puesto que v > 0.
    
    ρ = A \ ( ρ - D*F )    
end