Code
using Plots, LaTeXStrings, Printf
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),ylim=[0.0,1.5],
title=L"\rho_0", linewidth=2, thickness_scaling = 1.5,label="")
Ecuación de calor en 1D
Métodos Numéricos Avanzados. 2023-2024
David Gómez-Castro 
Ecuación de calor
Modelización
\[ \newcommand{\diff}{\mathrm{d}} \newcommand{\bm}[1]{\boldsymbol{#1}} \]
La ecuación de calor corresponde a \(F = -\nabla \rho\)
\[ \frac{\partial \rho}{\partial t} = \Delta \rho \]
También se puede escribir como \(F = - \rho \nabla \ln \rho\), y esto es interesante en algunos contextos. No lo haremos así.
Pongamos como dato inicial
\[ \rho_0 = \min\Big(1,\max(2-4|x|,0)\Big) \]
Elección de \(F_{i+\frac 1 2}\)
Puesto que \(F = -\nabla \rho\), una elección natural es \[ F_{i+\frac 1 2} = -\frac{\rho_{i+\frac 1 2} - \rho_{i-\frac 1 2}}{\Delta x} \] Resulta el método \[ \frac{d \rho_i}{d t} = - \frac{2 \rho_i (t) - \rho_{i+1} (t) - \rho_{i-1} (t) }{(\Delta x)^2} \]
Condiciones de contorno
Tenemos tres opciones sencillas:
Periódicas, igual que antes \(\rho_0 = \rho_N\) y \(\rho_{N+1} = \rho_1\).
Neumann: fijamos el flujo igual a cero \(F_{\frac 1 2} = F_{N+\frac 1 2} = 0\). Es decir, formalmente \(\rho_0 = \rho_{N+1} = 0\) y trabajar con todas las incógnitas interiores.
Dirichlet: fijamos \(\rho(t,-1) = \rho(t,1) = 0\). Podemos elegir \(\rho_1 = \rho_{N+1} = 0\) y trabajar con un punto menos, o ajustar la malla.
Método de Euler explícito
Método de Euler explícito
\[ \rho_i^{n+1} = \rho_i^n - \frac{\Delta t}{(\Delta x)^2} \left( 2 \rho_i^n - {\rho_{i+1}^n - \rho_{i-1}^n } \right) \]
No hay una condición CFL estrictamente hablando.
Si queremos que sea estable en el sentido \(|\rho_i^{n+1}| \le \max_j |\rho_j^n|\) basta que
\[ 2 \frac{\Delta t}{(\Delta x)^2} \le 1 \]
Esto obliga a que \(\Delta t\) se muy pequeño.
Euler explícito con condición periódica
Code
T = 2
N = 16;
x = range(-1.0,1.0,N+1); x = vec(x[1:end-1])
Δx = x[2] - x[1]
C = 0.5;
Δt = C*(Δx)^2;
ρ = ρ0(x);
@gif for n=1:ceil(T/Δt)
plot(x, ρ, ylim=[0.0,1.5],
title=@sprintf("t = %0.3f", float((n-1)*Δt)),
thickness_scaling = 1.5,linewidth=3,label="")
ρn = copy(ρ)
ρ[1] = ρn[1] - C*( 2*ρn[1] - ρn[N] - ρn[2] )
ρ[2:end-1] = ρn[2:end-1] - C*( 2*ρn[2:end-1] - ρn[1:end-2] - ρn[3:end] )
ρ[N] = ρn[N] - C*( 2*ρn[N] - ρn[1] - ρn[N-1] )
end
Pero tenemos
@show Δx, Δt;(Δx, Δt) = (0.125, 0.0078125)Método de Euler implícito
\[ \rho_i^{n+1} + \frac{\Delta t}{(\Delta x)^2} \left( 2 \rho_i^{n+1} - {\rho_{i+1}^{n+1} - \rho_{i-1}^{n+1} } \right) = \rho_i^n \]
Es incondicionalmente estable. Funciona incluso con un \(\Delta t\) muy grande.
Construímos \(A\) tal que \(A \bm \rho^{n+1} = \rho^n\).
Δt = 0.1; C = Δt / (Δx)^2
A = zeros(N,N)
for i=1:N A[i,i] = 1 + 2*C end
for i=1:N-1 A[i,i+1] = A[i+1,i] = -C end
A[N,1] = A[1,N] = -C;Para no almacenar matrices llenas de 0s podemos el paquete SparseArrays
using SparseArrays
A = sparse(A)16×16 SparseMatrixCSC{Float64, Int64} with 48 stored entries:
⎡⠻⣦⡀⠀⠀⠀⠀⠈⎤
⎢⠀⠈⠻⣦⡀⠀⠀⠀⎥
⎢⠀⠀⠀⠈⠻⣦⡀⠀⎥
⎣⡀⠀⠀⠀⠀⠈⠻⣦⎦Code
T = 2
ρ = ρ0(x);
anim = @animate for n=1:ceil(T/Δt)
global ρ
plot(x, ρ, ylim=[0.0,1.5], title=@sprintf("t = %0.3f", float((n-1)*Δt)), thickness_scaling = 1.5,linewidth=3,label="")
ρ = A \ ρ
end
gif(anim,fps=5)
Problema con fuente
La ecuación de calor con una fuente de calor (a veces llamada forzamiento) \[\begin{equation} \tag{F} \label{eq:forced} \frac{\partial \rho}{\partial t} - \frac{\partial^2 u}{\partial x^2} = f(t,x) \end{equation}\] se plantea en un dominio acotado con las mismas condiciones de contorno.
Si contamos con unos flujos \(F_{i+\frac 1 2} \sim - \nabla u(x_{i+\frac 1 2})\) (que incluyen la condición de contorno) podemos plantear la ecuación \[ \frac{\diff \rho_i} {\diff t} + \frac{F_{i+\frac 1 2} - F_{i-\frac 1 2}}{\Delta x} = f(t,x_i). \]
Los métodos anteriores quedan, en realidad, como métodos de diferencias finitas \[ \frac{\diff \bm \rho} {\diff t} (t) + A \bm \rho(t) = \begin{pmatrix} f(t, x_1) \\ \vdots \\ f(t, x_N) \end{pmatrix}. \]
Balance de masa
Con condiciones periódicas, la solución continua satisface \[ \frac{\diff}{\diff t} \int_{-1}^1 \rho(t,x) \diff x = \int_{-1}^1 f(t,x) \diff x. \]
En el método también es cierto que \[ \Delta x \sum_i \rho^{n+1}_i = \Delta \sum_i \rho^n_i + \Delta t \Delta x \sum_i f(t_n, x_i) . \]
Método de Euler implícito
Como ya hemos calculado las matrices necesarias, extendemos el código anterior
Code
T = 5
f(t,x) = sin(π*x) # Tiene \int_{-1}^1 f(t,x) dx = 0.
ρ = ρ0(x);
anim = @animate for n=1:ceil(T/Δt)
global ρ
plot(x, ρ, ylim=[0.0,1.5],
title=@sprintf("t = %0.3f", float((n-1)*Δt)),
thickness_scaling = 1.5,linewidth=3,label="")
ρ = A \ ( ρ + Δt * f.( n*Δt, x ))
end
gif(anim,fps=5)
Problema estacionario
El comportamiento habitual es que cuando \(t \to \infty\) las soluciones tienda a uno de los estados estacionarios que resuelven \[\begin{equation} \tag{E} \label{eq:stationary} - \frac{\partial^2 u}{\partial x^2} = f \end{equation}\]
Esto es cierto con condiciones periódicas si \(f = f(x)\) y \(\int_{-1}^1 f(x) \diff x = 0\).
Dada \(f\), \(\eqref{eq:stationary}\) tiene solución única si:
Con condición de contorno Dirichlet (\(u(0) = u(1) = 0\)) hay una única solución siempre.
Con condiciones periódicas o Neumann (\(- \frac{\partial u}{\partial x} (-1) = - \frac{\partial u}{\partial x} (1) = 0\) en \(\partial \Omega\)) hay solución si y sólo si \(\int_{-1}^1 f = 0\).
Además, si \(u\) es solución y \(C \in \mathbb R\), entonces \(u + C\) es solución.
Si preescribimos \(\int_{-1}^1 u = M\) ya hay única solución.
Si contamos con unos flujos \(F_{i+\frac 1 2} \sim - \nabla u(x_{i+\frac 1 2})\) (que incluyen la condición de contorno) podemos plantear la ecuación \[ \frac{F_{i+\frac 1 2} - F_{i-\frac 1 2}}{\Delta x} = f(x_i). \]
La condición de masa, caso de ser necesaria, viene dada por \[ \Delta x \sum_{i} \rho_i = M. \]
Problema estacionario numérico
Para resolver \(\eqref{eq:forced}\) teníamos que estudiar \[ \left(I + \frac{\Delta t}{\Delta x^2} D_2\right) \bm \rho^{n+1} = \bm \rho^n + \Delta t f(t_n). \]
Pero para resolver \(\eqref{eq:stationary}\) debemos resolver \[ \frac{1}{\Delta x^2} D_2 \bm \rho = \bm f. \]
El problema con condiciones periódicas no está bien determinado
D2 = zeros(N,N)
for i=1:N D2[i,i] = 2/(Δx^2) end
for i=1:N-1 D2[i+1,i] = D2[i,i+1] = -1/(Δx^2) end
D2[1,N] = D2[N,1] = -1/(Δx^2);
using LinearAlgebra
@show N-rank(D2);N - rank(D2) = 1El núcleo de \(D_2\) es precisamente el espacio asociado a \((1, \cdots, 1)\).
Para que haya solución los datos deben cumplir \(\displaystyle \sum_i f_i = 0\)
Hacemos el sistema compatible determinado cambiando la última ecuación
D2[N,1:N] = Δx * ones(1,N)
@show N-rank(D2);N - rank(D2) = 0(esto tiene sentido si \(\sum_i f_i \sim 0\)).
f(x) = sin(π*x)
b = zeros(N)
b[1:N-1] = f.(x[1:end-1])
b[N] = Δx*sum(ρ0(x)) # Masa total
ρ = D2 \ b
using LaTeXStrings
plot(x, ρ , title=latexstring("\$-\\frac{\\partial^2 \\rho}{\\partial x^2} = \\sin(\\pi x)\$ \n con condición periódica"), thickness_scaling=1.5,linewidth=2,label="",xlabel=L"x",ylabel=L"\rho", ylim=(0,1))