Ecuación de calor en 1D

Métodos Numéricos Avanzados. 2022-2023

Author

David Gómez-Castro

Ecuación de calor

Modelización

La ecuación de calor corresponde a \(F = -\Delta \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

  • 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