Ecuación de calor en 1D

Métodos Numéricos Avanzados. 2023-2024

David Gómez-Castro

https://gomezcastro.xyz

Universidad Complutense de Madrid

Ecuación de calor

Modelización

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

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="")

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

\[ \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) = 1

El 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))