Solving heat equation with python (NumPy)

Question

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

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),

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?

This is a programming site, not a physics site, so you need to give the relevant equations, both for the code you show and for the code you want. Many of us do not understand what "one end is insulated (not flux)" means. — Rory Daulton, Mar 24 '18 at 11:15
@RoryDaulton you are absolutely right. I modified the question with the corresponding equations. — Googlebot, Mar 24 '18 at 11:34
Your heat equation is incorrect. The LHS is ∂u/∂t, not ∂u/∂x. Otherwise your heat doesn't change with time at all, and there is no flux. — dROOOze, Mar 24 '18 at 11:37
@droooze sorry for the typo, I was looking for a quick way to add math equations here. — Googlebot, Mar 24 '18 at 11:45
My point was that googling Neumann boundaries should give you the answer... — Andras Deak -- Слава Україні, Mar 24 '18 at 11:56
Since you're using a finite difference approximation, see [this](https://en.wikipedia.org/wiki/Finite_difference_method). All you have to do is to figure out what the boundary condition is in the finite difference approximation, then replace the expression with `0` when the finite difference approximation reaches these conditions. This is not really a `python` or an implementation question, since you haven't yet figured out the FD discretisation before you've attempted to code it (or at least haven't shown the discretisation of the changed BC in the question). — dROOOze, Mar 24 '18 at 12:14
Take a look at [the example in the `scipy.integrate` tutorial](https://docs.scipy.org/doc/scipy/reference/tutorial/integrate.html#solving-a-system-with-a-banded-jacobian-matrix). It shows the spatial discretization for a system of PDEs with Neumann ("no flux") boundary conditions. (The time evolution is solved using `scipy.integrate.odeint`; the tutorial is an example of the "method of lines".) — Warren Weckesser, Mar 24 '18 at 13:39

score 9 · Accepted Answer · answered Mar 24 '18 at 17:12

9

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

answered Mar 24 '18 at 17:12

Very subtle solution. It solved my problem, but have a question out of curiosity. This code applies to the boundary of `u_x=0` (no flux). How does it work for non-zero Neumann boundaries such as `u_x= constant` (constant flux)? In principle, these two cases should be similar, as the former is in the category of the latter, in general. – Googlebot Mar 24 '18 at 17:59
1

If u_x is A then `(T[n, j+1] - T[n, j-1]) / (2*dx) = A` implies `T[n, j+1]` is `T[n, j-1] + 2*A*dx`, so that is the value that replaces `T[n, j+1]` in the formula on the last line. – Mar 24 '18 at 18:18
Sorry for too many questions, but I am fascinated by the simplicity of this solution and my stupidity to comprehend the whole picture. Thus, I try to inspect different boundary conditions. How will be the case of a constant flux from one side and no flux on the other side? `u_x(0,t)=A` and `u_x(l,t)=0`? I struggle to comprehend how the code works when `T0` is not a contact value. – Googlebot Mar 25 '18 at 17:20
1

There is no T0 in this case. The array can be originally filled with anything, say zeros. Then the initial values are filled in. After that, the diffusion equation is used to fill the next row. On the left boundary, when j is 0, it refers to the ghost point with j=-1. Because of the boundary condition, T[n, j-1] gets replaced by T[n, j+1] - 2*A*dx when j is 0. – Mar 25 '18 at 17:38

Solving heat equation with python (NumPy)

1 Answers1