Volúmenes Finitos y la ecuación de transporte

Métodos Numéricos Avanzados. 2023-2024

David Gómez-Castro

Universidad Complutense de Madrid

Volúmenes Finitos

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

integrando en \([x_{i-\frac 1 2}, x_{i+\frac 1 2}]\) tomamos la semi-discretización espacial \[\begin{equation} \frac{\diff}{\diff t} \rho_i (t) = - \frac{ \flux_{i+\frac 1 2} (t) - \flux_{i-\frac 1 2} (t) }{x_{i+\frac 1 2} - x_{i-\frac 1 2}}, \qquad i = 1, \cdots, N. \end{equation}\]

En nuestro tipo de leyes de conservación, será frecuente que estudiemos EDPs cuyas soluciones exactas sean positivas: \(\rho \ge 0\).

Muchos de los razonamientos se pueden generalizar al caso de cambio de signo.

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*}\]

Conservación de masa es equivalente a que \[ F_{N+\frac 1 2} (t) = F_{\frac 1 2} (t) \]

En algunas ocasiones se toman ambos flujos igual a 0 para que el problema esté bien determinado.

En otras ocasiones tomaremos condiciones periódicas en \(\rho\), es decir \[ \rho_0 = \rho_N, \qquad \rho_1 = \rho_{N+1}. \]

Ecuación de Transporte

Ecuación de transporte

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

\[\frac{\partial \rho}{\partial t} + a\frac{\partial \rho}{\partial x} = 0\]

Soluciones exactas

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

En este caso, buscamos curvas \(X(t)\) tales que \(\rho(t, X(t)) = \rho_0 (y)\), donde \(a(0) = 1\). \[ \partial_t \rho + \frac{\partial \rho}{\partial x} \frac{\diff X}{\diff t} = 0. \]

Luego \(\frac{\diff X}{\diff t}(t) = a\). Y \(X(t,y) = y + at\).

Finalmente resolviendo \(x = y + at\) tenemos \[ \rho(t,x) = \rho_0 (x - a t). \] Es decir, solución en el punto \((t,x)\) depende sólo del valor inicial en \((0,x - at)\).

Restricción a un intervalo cerrado

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

Puede, sin embargo, admitir condiciones de periódicas. Por ejemplo, para \(t \in (0,1)\) y \(a > 0\) \[ \rho(t,0) = \rho(t,1) = \rho_0(1 - at). \]

Imponer condiciones de flujo periódico consiste en decir \[ a \rho(t,1) = F(t,1) = F(t,0) = a \rho(t,0). \]

Desde el punto de vista exacto, es de nuevo imponer condiciones de periódicas en \(\rho\).

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

Para que sea explícito \(F_{i+\frac 1 2}^{n+1} = G_{i+\frac 1 2} ( \rho_i^n, \rho_{i+1}^n )\).

Es decir \(F_{i+\frac 1 2}^{n+1} = a\rho_i^n , a\rho_{i+1}^n\) u otro valor.

Experimentación numérica

Hagamos unos experimentos para ver qué elección de \(\rho_i\) funciona mejor.

using Plots

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

\(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. \]

Tomemos \(\rho^0 = (0,1,0,0,\cdots)\). Entonces \(\rho^1 = ( 0,1- a\frac{\Delta t}{\Delta x} , a\frac{\Delta t}{\Delta x} , 0, \cdots )\).

Así \(\rho^1 \ge 0\) si y sólo si \[ 0 \le a\frac{\Delta t}{\Delta x} \le 1. \]

La primera desigualdad dice que este método sólo funciona si \(a \ge 0\). Esto se llama up-winding

La segunda desigualdad dice que \(\Delta t\) debe ser “pequeño” \[ \Delta t \le \frac 1 a \Delta x. \]

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

Ejercicio. Comprobar que para \(\rho^{n+1}_i = \rho^n_i - a\frac{\Delta t}{\Delta x} ( \rho_{i}^n - \rho_{i-1}^n )\)

la condición es \(0 \le -a\frac{\Delta t}{\Delta x} \le 1\).

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.

Condición Courant–Friedrichs–Lewy

Tomamos \(a>0\) y el método con up-winding. Se cumple el siguiente mantra:

Un método numérico sólo puede ser convergente si el dominio de dependencia numérico contiene al dominio de dependencia la EDP, al menos cuando \(\Delta t\) y \(\Delta x\) tienden a \(0\).

Code

using Plots, Kronecker, LaTeXStrings

a = 0.9; T = 1; 
@gif for log2Δt=-3:1:0 , log2Δx=-3:1:1    
    Δt=2.0^log2Δt; Δx=2.0^log2Δx

    if Δt <= Δx/a
        inequality = "\\leq"
        color = :green
    else
        inequality = ">"
        color = :red
    end

    plot(ylim=(-.1,1.1), xlim=(-2.0,2.0),
        aspect_ratio=:equal,
        xlabel=L"x", ylabel=L"t", legend=:outerbottom, 
        title=latexstring("\\Delta t=$Δt " * inequality * " \\Delta x / a, \\Delta x = $(Δx), a = $a"),
        titlefontcolor=color)

    plot!([0.0,-a],[1.0,0.0],label="Característica",color=:blue,
        linewidth=3.0)

    for i in -(2/Δx):(2/Δx), j in 0:(T/Δt)
      scatter!([i*Δx], [j*Δt],  color=:white , label="")
    end
    
    for i in 0:(T/Δt), j in 0:(T/Δt-i)
      scatter!([-j*Δx], [i*Δt],  color=:red , label="")
    end
    scatter!([0], [1.0],  color=:red , label="Dependencia del método numérico")

end fps=0.5

Queremos que \(0 \le a \frac{\Delta t}{\Delta x} \le 1\) para que todo funcione bien.

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

Simulación que no cumple CFL

Δt = 1.1*Δ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

Simulación buena

N  = 128; T  = .5
x  = range(-1.0,1.0,N+1); x = vec(x[1:end-1])
Δx = x[2] - x[1]
Δt = Δx/5

ρ0 = min.(1,max.(-4*abs.(x).+2,0)); 
ρ = 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 every 2

Pregunta

Calcule la solución explícita de \(\frac{\partial \rho}{\partial t} + a \frac{\partial \rho}{\partial x} = 0\).

Cuando \(a \frac{\Delta t}{\Delta x} = 1\) ¿qué relación tiene el método con la solución explícita?

Up-winding

Puesto que \(F = \rho v\), elegida un aproximación \(v_{i+\frac 12}\), en la frontera de volúmenes \[ F_{i+\frac 1 2} = \Big(A_{i+\frac 1 2}\rho_i + B_{i+\frac 1 2} \rho_{i+1} \Big )v_{i+\frac 1 2}. \]

con \(A_{i+\frac 1 2}, B_{i+\frac 1 2} \ge 0\).

Al introducirlo en el método explícito obtenemos

\[ \begin{aligned} \frac{d \rho_i}{d t} &= \underbrace{ \frac{1}{\Delta x} \left( - A_{i+\frac 1 2} v_{i+\frac 1 2} + B_{i-\frac 1 2} v_{i-\frac 1 2} \right) }_{C_i} \rho_i \\ &\qquad + \frac{1}{\Delta x} \left( - B_{i+\frac 1 2} v_{i+\frac 1 2} \right) \rho_{i+1} + \frac{1}{\Delta x} \left( + A_{i-\frac 1 2} v_{i-\frac 1 2} \right) \rho_{i-1}. \end{aligned} \]

\[\begin{align*} \rho_i (t) &= e^{C_i t}\rho_i (0) \\ &+ \frac{\Delta t}{\Delta x} \int_0^t e^{C_i (t-s)} \left[ \left( - B _ {i+\frac 1 2} v _ {i+\frac 1 2} \right) \rho_{i+1} (s) + \left( + A _ {i-\frac 1 2} v _ {i-\frac 1 2} \right) \rho _ {i-1}(s) \right] d s \end{align*}\]

Tenemos

Así, cuando \(v_{i+\frac 1 2} > 0\) tenemos que elegir \(B_{i+\frac 1 2} = 0\),

y cuando \(v_{i-\frac 1 2} < 0\) tenemos que elegir \(A_{i-\frac 1 2} = 0\).

Por otro lado, el otro parámetro debe fijarse a 1, para que el esquema sea consistente.

Con la notación \(v_+ = \max\{ v , 0\}\) y \(v_- = \min\{v,0\}\) (es decir \(v = v_+ + v_-\)) se tiene la elección de up-winding \[ \tag{Up-wind} \label{eq:up-winding} F_{i+\frac 1 2} = \rho_i \left( v_{i+\frac 1 2} \right)_+ + \rho_{i+1} \left (v_{i+\frac 1 2} \right)_-. \]

Ejercicio. Método de Euler implícito

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

Despejando \[ \left( 1 + a \frac{\Delta t}{\Delta x} \right) \rho_i^{n+1} - a \frac{\Delta t}{\Delta x} \rho_{i-1}^{n+1} = \rho_i^n. \]

El sistema queda \[ A \rho^{n+1} = \rho^n, \qquad A = \begin{pmatrix} \left( 1 + a \frac{\Delta t}{\Delta x} \right) & & \cdots & - a \frac{\Delta t}{\Delta x} \\ - a \frac{\Delta t}{\Delta x} & \left( 1 + a \frac{\Delta t}{\Delta x} \right) \\ && \ddots \\ && - a \frac{\Delta t}{\Delta x} & \left( 1 + a \frac{\Delta t}{\Delta x} \right) \end{pmatrix} \]

Code

ρ = copy(ρ0); Δt = Δx/2

N = length(x)
A = zeros(N,N)
for i=1:N A[i,i]   = 1.0 + a*Δt/Δx end
for i=2:N A[i,i-1] = - a*Δt/Δx     end
A[1,N]= - a*Δt/Δx # Periodicidad

@gif for n=1:ceil(T/Δt)
    global ρ
    
    ρ = A \ ρ
    
    plot(x,ρ,
        ylim=(0.0,1.2),
        title="t=$(round(n*Δt,digits=3))",

        linewidth=2,thickness_scaling = 1.5,label="") 
end

Convergencia

Orden de convergencia

Para la ecuación de transporte en \((0,1)\) con condiciones periódicas, el método de Euler explícito con up-wind, y \(\Delta t = \frac 1 2 \Delta x\), calcular \(\alpha\) tal que \[ error \sim (\Delta x)^\alpha \]

Aproximamos el error por la norma \(L^1\) a tiempo \(T\), donde \(x_1, \cdots, x_M\) y \(t_1, \cdots, t_N = T\). \[ error = \Delta t \Delta x \sum_{n=1}^N \sum_{i=1}^M |\rho_i^n - \rho_{exacta} (t_n,x_i) | \]

Sugerencia: calcular la gráfica \((\log \Delta x, \log error)\).

Código

using Plots, LaTeXStrings, LinearAlgebra, CurveFit

T = .5; 
R = 0.5 # Coeficiente Δt = R Δx

ρ0f(x) = min(1,max(2.0-8*abs(x - 0.5),0.0)); 
ρexacta(t,x) = ρ0f.( mod.( x .- t , 1) );

Code

function num_error(N,R)
    x  = range(0.0,1.0,N+1); x = vec(x[1:end-1])
    Δx = x[2] - x[1]
    Δt = R*Δx
    ρ  = ρ0f.(x) 

    err = 0
    t = 0
    while t < T
        err = err + Δt * Δx * sum( abs.( ρ - ρexacta(t,x) ) )

        ρn = copy(ρ)
        ρ[1]     = ρn[1]     - Δt/Δx * ( ρn[1]     - ρn[N] ) 
        ρ[2:end] = ρn[2:end] - Δt/Δx * ( ρn[2:end] - ρn[1:end-1] ) 
        t += Δt;

        plot(x,ρ)    
    end;

    return (Δx, err)
end;

Code

Ns = 2 .^(4:15)
Δxs = zeros(size(Ns));  num_errors = zeros(size(Ns))
for n = 1:length(Ns)
    Δxs[n], num_errors[n] = num_error(Ns[n],R)
end

scatter(log.(Δxs),log.(num_errors), 
    label="Numerical error", 
    xlabel=L"\log (\Delta x)", 
    ylabel=L"\log(error)",
    legend=:topleft,
    thickness_scaling=1.5,
    title= latexstring("\\Delta t = $(R) \\Delta x")
);

b,a = linear_fit(log.(Δxs), log.(num_errors))

plot!( log.(Δxs) , a*log.(Δxs) .+ b, 
    label=latexstring("C (\\Delta x)^{ $(a) }")  )

Comentarios finales

CFL

Por extensión, algunos autores llaman condición CFL a cualquier relación \(\Delta t \le C (\Delta x)^\alpha\) que se necesaria para la convergencia.

Ejemplo. (ecuación de calor) Nótese que la solución exacta de ecuación de calor se escribe, es bien sabido que aunque \(\rho_0\ge 0\) sea de soporte compacto, \(\rho(t) > 0\) en todo el dominio. Al haber propagación de velocidad infinita, el dominio de dependencia de cualquier punto es total. Sin embargo, los métodos no toman en consideración todos los puntos de dominio, y por tanto no puede haber una condición de CFL. Sin embargo, es habitual en la literatura hablar de “condiciones CFL” para métodos explícitos.

CFL en la ecuación de ondas

La ecuación de ondas \[\begin{equation} \frac{\partial^2 \rho}{\partial t^2} = D \frac{\partial ^2 \rho}{\partial x^2} \end{equation}\] El dominio de dependencia del punto \((t,x)\) es el intervalo \((x-\sqrt D t,x+\sqrt D t)\) (véase la fórmula de D’Alambert). Esta ecuación puede discretizarse usando el método de volúmenes finitos para obtener \[\begin{equation*} \frac{\partial^2 \rho_i}{\partial t^2} = D \frac{\rho_{i+1} -2 \rho_i + \rho_{i-1}}{(\Delta x)^2} \end{equation*}\] Si empleamos la fórmula de diferencias centradas para aproximar la segunda derivada \[\begin{equation} \label{eq:wave punto medio} \frac{\rho_i^{n+1} - 2 \rho_i^n + \rho_i^{n-1} }{(\Delta t)^2}= D \frac{\rho_{i+1}^n -2 \rho_i^n + \rho_{i-1}^n}{(\Delta x)^2} \end{equation}\] Dependiendo del método que se emplee para la discretización temporal, se espera que el dominio de dependencia sea \((x - n \Delta x , x+ n \Delta x)\). De modo que la condición CFL es \(\sqrt{D} \frac{\Delta t}{\Delta x} \le 1\).

Recordatorio de estabilidad lineal

Al aplicar métodos monopaso en tiempo obtenemos \[ A_1 \mathbf \rho^{n+1} = A_0 \mathbf \rho^n. \] De esta forma, la estabilidad depende de los autovalores de \(A_1^{-1} A_0\). Recordamos el teorema fundamental > Un esquema lineal consistente es convergente si y sólo si es estable.

Ejemplo. Volviendo sobre \[ \rho^{n+1}_i = \left( 1 - a \frac{\Delta t}{\Delta x} \right) \rho_i^n + a \frac{\Delta t}{\Delta x}\rho_{i-1}^n. \] tenemos, para \(a > 0\) \[ \rho_i^{n+1} = \rho_i^n - a\frac{\Delta t}{\Delta x} ( \rho_{i}^n - \rho_{i-1}^{n} ) = \left(1 -a R \right) \rho_i^n + aR \rho_{i-1}^n, \] utilizando el coeficiente de Courant \(R = \frac{\Delta t}{\Delta x}\). Matricialmente, con la condición de flujo nulo \[ \mathbf \rho^{n+1} = A_0 \mathbf \rho^n, \qquad A_0 = \begin{pmatrix} 1-aR && & aR \\ aR & 1-aR \\ && \ddots \\ &&aR & 1-aR \end{pmatrix} \] Así, \(\|\rho^n\| \le |||A |||^n \| \rho^0 \|\). Si elegimos como norma la norma 1, entonces se puede calcular directamente que \[ ||| A |||_1 = |1-aR| + |aR|. \] Como hemos elegido \(a>0\), si \(1-aR < 0\) entonces tenemos \(||| A |||_1 = 2aR - 1 > 1\). Así, la condición de estabilidad es \(|a|\frac{\Delta t}{\Delta x} \le 1\), precisamente.