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.