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)