Métodos Numéricos Avanzados. 2023-2024
\[ \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.
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}. \]
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\]
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)\).
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\).
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.
Hagamos unos experimentos para ver qué elección de \(\rho_i\) funciona mejor.
Tenemos el método \[ \rho^{n+1}_i = \rho^n_i - a\frac{\Delta t}{\Delta x} ( \rho_{i}^n - \rho_{i-1}^n ) \]
\[\rho^{n+1}_i = \rho^n_i - a\frac{\Delta t}{\Delta x} ( \rho_{i+1}^n - \rho_{i}^n )\]
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. \]
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\).
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.
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.
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\).
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.
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
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?
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)_-. \]
\[ \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} \]
ρ = 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
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)\).
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;
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) }") )
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.
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\).
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.
Métodos Numéricos Avanzados. David Gómez-Castro