Ecuación de Burgers

Métodos Numéricos Avanzados. 2023-2024

Author

Affiliation

David Gómez-Castro

Universidad Complutense de Madrid

Ecuación de Burgers

La ecuación de Burgers corresponde a un fluido que se mueve con velocidad \(v = \rho/2\), es decir \(F = v \rho\) y

\[ \frac{\partial \rho}{\partial t} = - \frac{\partial}{\partial x} \frac {\rho^2} 2 \]

El \(2\) procede de \(\frac{\partial \rho}{\partial t} =- \rho \frac{\partial \rho }{\partial x}\).

Dato inicial

\[ \rho_0 = \min\Big(1,\max(2-4|x|,0)\Big) \]

using Plots, LaTeXStrings, Printf

plot_ρ(x,ρ,t) = plot(x, ρ, ylim=[0.0,1.5],
    title=@sprintf("t = %0.3f", t),
    linewidth=2,thickness_scaling = 1.5,label="")

plot_ρ (generic function with 1 method)

N  = 64; 
x  = range(-1.0,1.0,N+1); x = vec(x[1:end-1])
Δx = x[2] - x[1]
ρ0(x) = min.(1,max.(-4*abs.(x).+2,0))

plot_ρ(x,ρ0(x),0)

Soluciones exactas

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) \]

Luego \(\rho_t (y + \rho_0(y) t) = \rho_0(y)\)

y = copy(x)

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

Las características ya no son paralelas. ¡Se cortan y aparecen discontinuidades!

Si a un punto \((x,t)\) llegan dos características, es decir \(x = y_1 + \rho_0(y_1) t = y_2 + \rho_0(y_2) t,\)

para calcular \(\rho_t (x)\) tengamos que elegir entre \(\rho_0(y_1)\) y \(\rho_0(y_2)\).

Método de Euler explícito

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.

Resulta el método \[ \rho_{i}^{n+1} = \rho_i^n - \frac{\Delta t}{2 \Delta x} \left( {\rho_i^n} \left( \alpha \rho_i^n + (1-\alpha)\rho_{i+1}^n \right) - {\rho_{i-1}^n} \left( \alpha \rho_{i-1}^n + (1-\alpha)\rho_{i}^n \right) \right) \]

Condición CFL

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

Puesto que \(0 \le v_{i+\frac 1 2} (t) \le \frac {\max_i  \rho_i (t)} 2\)

una condición CFL del método de velocidad constante sería

\[ \frac {\max_i  \rho_i (t)} 2\frac{\Delta t}{\Delta x} \le 1 \]

Dado que la solución exacta cumple \(0 \le \rho_t \le 1\),

cabe esperar que \(\frac 1 2 \frac{\Delta t}{\Delta x} \le 1\) sea suficiente.

Experimento numérico.

      Calcular la solución numérica con \(\alpha = \frac 1 2\), y responder a la siguientes preguntas:

• ¿Es suficiente esta condición CFL para que el método converja?

• ¿Es posible calcular una condición CFL para este problema?

---

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

Las oscilaciones producen que se pierda CFL.

Así \(\rho_i\) deja de ser positivo.

De forma que estamos perdiendo el up-winding, y violando la condición CFL.

Monotonía

Monotonía

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

Las soluciones exactas cumplen \(\rho(t) \le \overline \rho = 1\) (principio de comparación).

Aunque se tiene CFL y \(\rho_0 \le 1\), no se tiene \(\rho_i^n \le \overline \rho_i^n = 1\).

Esto se conoce como pérdida de la monotonía.

¿Cómo preservar la monotonía en un método explícito?

      Tomemos un método explícito genérico

\[ \rho_i^{n+1} = H(\rho_{i+1}^n,\rho_i^n, \rho_{i-1}^n). \]

      Supongamos que \(\rho_i^n \le 1\) para todo \(i\), entonces queremos que

\[ 1 \ge \rho_{i}^{n+1} = H(\rho_{i+1}^n, \rho_{i-1}^n, \rho_{i-1}^n). \]

Para preservar la monotonía, \(H\) debe ser monótona no-decreciente en cada componente.

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) \]

Luego \[ H(p,q,r) = \left( 1 - \frac{\Delta t}{\Delta x} \frac{ {(1-\alpha) p + \alpha q}}{2} \right)q + \frac{\Delta t }{\Delta x} \frac{{(1-\alpha) q + \alpha r}}{2} r \]

Monotonía

Volvemos sobre \[ H(p,q,r) = \left( 1 - \frac{\Delta t}{\Delta x} \frac{ {(1-\alpha) p + \alpha q}}{2} \right)q + \frac{\Delta t }{\Delta x} \frac{{(1-\alpha) q + \alpha r}}{2} r \]

Debemos pedir que \[ 0 \le \frac{\partial H}{\partial p} = - \frac{\Delta t}{\Delta x} \frac{{1-\alpha}}{2} q \]

Necesariamente \(\alpha = 1\). Reducimos la fórmula a este caso

\[ H(p,q,r) = \left( 1 - \frac{\Delta t}{\Delta x} \frac q 2 \right)q + \frac{\Delta t }{\Delta x} \frac{r^2}{2}. \]

Debemos pedir también que \[ 0 \le \frac{\partial H}{\partial q} = 1 - \frac{\Delta t}{\Delta x} q \]

Esta condición dice que para todo \(i\) \[ \frac{\Delta t}{\Delta x} \rho_i^n \le 1. \] ¡Más exigente que CFL!

Por último \[ 0 \le \frac{\partial H}{\partial r} = \frac{\Delta t}{\Delta x} r \] que no impone nada nuevo.

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 pueden construir esquemas monótonos con aproximaciones monótonas del flujo, sin tener necesariamente en cuenta el up-winding de la ecuación de transporte. En tal caso \[ F_{i+\frac 1 2} = g_{i+\frac 1 2}^{n+1} ( \rho_1, \cdots, \rho_M ) . \]

Estos métodos son consistentes si \[ g_{i+\frac 1 2} (s, \cdots, s) = \Phi (s) . \] Habitualmente, \(g\) no depende de todas las componentes de \(\rho\), si no de una cantidad finita.

---

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) \]

El método es monótono si y sólo si \[ \frac{\partial g}{\partial a} \le 0, \qquad \frac{\partial g}{\partial b} \ge 0, \]

y, además, se cumple la condición \[ \frac{\Delta t}{\Delta x} \left( \frac{\partial g}{\partial b} (\rho_{i+1}^n, \rho_i^n) - \frac{\partial g}{\partial a} (\rho_{i}^n, \rho_{i-1}^n) \right) \le 1 . \]

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

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) \]

Es decir, en cada paso queremos resolver un sistema de ecuaciones.

Planteado como ecuación no lineal \[ 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) \]

Resolución de sistema no lineales

Para resolver \(f(\xi) = 0\) el método de Newton propone resolver el problema linealizado \[ h(\xi^n) + Dh(\xi^n) (\xi^{n+1} - \xi^n ) = 0. \]

Es decir \[ \xi^{n+1} = \xi^n - [Dh(\xi^n)]^{-1} h(\xi^n). \]

Normalmente iteramos hasta que \(\max_i|\xi_i^n - \xi_i^{n+1}| < \mathrm{tol}.\)

function fixed_point(ρ,h,Dh)
    tol  = 1e-5; max_iter = 1e3
    ρn   = copy(ρ)
    iter = 1
    δ    = ones(size(ρ))
    while (maximum(abs.(δ)) > tol) & (iter < max_iter)       
        δ = Dh(ρ,ρn) \ h(ρ,ρn)
        ρ = ρ - δ
        iter += 1
    end
    if (iter == max_iter)
        warning("Punto fijo no convergió")
    end

    return ρ
end;

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) \]

Luego \[ \frac{\partial h_i}{\partial \xi_i} = 1 + \frac{\Delta t}{\Delta x} \xi_i \qquad \qquad \frac{\partial h_i}{\partial \xi_{i-1}} = - \frac{\Delta t}{\Delta x} \xi_{i-1} \]

Es decir \[ Dh(\xi) = I + \frac{\Delta t}{\Delta x} \begin{pmatrix} 1 & 0 &\cdots & -1\\ -1 & 1 \\ &&\ddots \\ &&& - 1& 1 \end{pmatrix} \begin{pmatrix} \rho_1 \\ & \rho_2 \\ &&& \ddots \\ &&&&\rho_N \end{pmatrix} \]

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

Escribimos las derivadas de h

Code

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

Comentarios finales

Ecuación de Burgers viscosa

La ecuación de Burgers puede formar discontinuidades. Una forma de evitarlo es añadir un término “viscoso” \[ \frac{\partial \rho}{\partial t} + \frac{\partial}{\partial x} \frac{\rho^2}{2} = \varepsilon \Delta \rho \] Al límite \(\varepsilon \to 0\) se lo llama de vanishing viscosity

Un método de Euler explícito en este caso \[ \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 - 2 \rho_i^n + \rho_{i-1}^n }{(\Delta x)^2} \]

Empleando las técnicas anteriores, el método es monótono si \[ 2 \varepsilon \frac{\Delta t}{(\Delta x)^2} + \rho_i^n \frac{\Delta t}{\Delta x} \le 1. \]

Code

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

Δt = 0.9 / ( 2 * ϵ / (Δx)^2 + 1 / (Δx) ) ; # Condicion de monotonía
ρ = ρ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)))

    ρn = copy(ρ)

    F = ρn.^2/2 # puesto que v > 0.
    
    per(n) = mod(n-1,N)+1
    ρ = ρ - Δt/Δx * ( F - F[per.(0:N-1)] ) + ϵ*Δt/((Δx)^2)*( ρn[per.(2:N+1)] - 2*ρ + ρ[per.(0:N-1)]   )
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

Bibliografía

References

Eymard, R., Gallouët, T. & Herbin, R. (2000) Finite volume methods. In: Handbook of numerical analysis. Elsevier. pp. 713–1018. doi:10.1016/S1570-8659(00)07005-8.