2D heat equation

INTRODUCTION


Unknown function u(x,y,t), x in [0,L], y in [0,L], t in [0,T],

heat equation u_t = D * (u_xx + u_yy).

Initial condition u(x,y,0) = f(x,y).
Boundary condition u(boundary,t) = g(boundary,t).

x = i*dx, i = 0..Nx (Nx+1 points), Nx*dx = L,
y = j*dy, j = 0..Ny (Ny+1 points), Ny*dy = L, dx=dy, Nx=Ny,
t = k*dt, k = 0..Nt (Nt+1 points), Nt*dt = T,

u_t = [ u(i,j,k+1) - u(i,j,k) ] / dt   # forward difference O(dt)
u_xx = [ u(i-1,j,k) -2*u(i,j,k) + u(i+1,j,k) ] / (dx*dx)   # central difference O(dx^2)
u_yy = [ u(i,j-1,k) -2*u(i,j,k) + u(i,j+1,k) ] / (dy*dy)   # central difference O(dy^2)

r = D*dt / (dx*dx), formula is stable if r ≤ 0.25 (?)

u(i,j,k+1) - u(i,j,k) = r * [
u(i-1,j,k) -2*u(i,j,k) + u(i+1,j,k)
+ u(i,j-1,k) -2*u(i,j,k) + u(i,j+1,k) ]

Final formula (inside mesh)
u(i,j,k+1) = r*u(i-1,j,k) + r*u(i+1,j,k) + r*u(i,j-1,k) + r*u(i,j+1,k) + (1-4*r)*u(i,j,k)

Initial condition (t=0, k=0)
for i in range(Nx+1):
    for j in range(Ny+1):
        u(i,j,0) = f(i*dx,j*dy)

(a) Constant boundary conditions
for k in range(Nt+1):
    for i in range(Nx+1):
        u(i,0,k) = f(i*dx,0)
        u(i,Ny,k) = f(i*dx,L)
    for j in range(Ny+1):
        u(0,j,k) = f(0,j*dy)
        u(Nx,j,k) = f(L,j*dy)

(b) Periodic boundary conditions - final formula
for k in range(Nt+1):
    for i in range(Nx+1):
        for j in range(Ny+1):
            u(i,j,k+1) = [ r * u((i-1+Nx) % (Nx+1), j, k) 
                         + r * u((i+1) % (Nx+1), j, k) 
                         + r * u(i, (j-1+Ny) % (Ny+1), k) 
                         + r * u(i, (j+1) % (Ny+1), k) 
                         + (1-4*r) * u(i,j,k) ]

How to improve derivative estimates when using finite divided differences:
1. decrease the step size,
2. use a higher-order formula that employs more points,
3. use Richardson extrapolation.