Search code examples
pythonnumpynumerical-methods

Solving heat equation

I solve the heat equation for a metal rod as one end is kept at 100 °C and the other at 0 °C as

enter image description here

import numpy as np
import matplotlib.pyplot as plt

dt = 0.0005
dy = 0.0005
k = 10**(-4)
y_max = 0.04
t_max = 1
T0 = 100

def FTCS(dt,dy,t_max,y_max,k,T0):
    s = k*dt/dy**2
    y = np.arange(0,y_max+dy,dy) 
    t = np.arange(0,t_max+dt,dt)
    r = len(t)
    c = len(y)
    T = np.zeros([r,c])
    T[:,0] = T0
    for n in range(0,r-1):
        for j in range(1,c-1):
            T[n+1,j] = T[n,j] + s*(T[n,j-1] - 2*T[n,j] + T[n,j+1]) 
    return y,T,r,s

y,T,r,s = FTCS(dt,dy,t_max,y_max,k,T0)

plot_times = np.arange(0.01,1.0,0.01)
for t in plot_times:
    plt.plot(y,T[t/dt,:])

If changing the Neumann boundary condition as one end is insulated (not flux),

enter image description here

then, how the calculating term

T[n+1,j] = T[n,j] + s*(T[n,j-1] - 2*T[n,j] + T[n,j+1])

should be modified?

Solution

  • A typical approach to Neumann boundary condition is to imagine a "ghost point" one step beyond the domain, and calculate the value for it using the boundary condition; then proceed normally (using the PDE) for the points that are inside the grid, including the Neumann boundary.

    The ghost point allows us to use the symmetric finite difference approximation to the derivative at the boundary, that is (T[n, j+1] - T[n, j-1]) / (2*dy) if y is the space variable. Non-symmetric approximation (T[n, j] - T[n, j-1]) / dy, which does not involve a ghost point, is much less accurate: the error it introduces is an order of magnitude worse than the error involved in the discretization of the PDE itself.

    So, when j is the maximal possible index for T, the boundary condition says that "T[n, j+1]" should be understood as T[n, j-1] and this is what is done below. 

    for j in range(1, c-1):
    T[n+1,j] = T[n,j] + s*(T[n,j-1] - 2*T[n,j] + T[n,j+1])  # as before
j = c-1 
T[n+1, j] = T[n,j] + s*(T[n,j-1] - 2*T[n,j] + T[n,j-1])  # note the last term here