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;Volúmenes Finitos y la ecuación de transporte
Métodos Numéricos Avanzados. 2023-2024
David Gómez-Castro 
Volúmenes Finitos
Ecuación de continuidad
\[ \newcommand{\diff}{\,\mathrm{d}} \DeclareMathOperator{\diver}{div} \newcommand{\flux}{{F}} \]
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.
\(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
\[\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*}\]
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.