I am new to solving PDEs with finite difference methods. I implemented the code below to solve the heat equation following the explicit scheme, but when I plot the result I am suprised that the decay in time is very small. Could anyone tell me what I am doing wrong ?
Here is my code in Python. Below I did plot for four steps in time (2, 50, 100, 150) what the heat looks as a function in space :
from __future__ import division import matplotlib.pyplot as plt import numpy as np import math L =1. N =200 J =200 k = 1./4. dx = float(L)/float(N) T = 0.5 * float(J)/k * dx * dx - 0.00001 #ensures that alpha <0.5 dt = float(T)/float(J) alpha = k * dt / (dx*dx) # stability condition: alpha <0.5 print ("stability condition (alpha<0.5) --> alpha = ", alpha) u = np.zeros((N,J)) def f(phi): return math.cos(np.pi*phi/float(L)) x = np.arange(0,L,dx) t = np.arange(0,T,dt) for n in range(0,N): u[n] = 6.*math.sin(np.pi*x[n]/float(L)) for j in range(0, J): u[j] = u[N-1][j] = 0. for j in range(0, J-1): for n in range(1,N-1): u[n][j+1] = (1.-2.*alpha)*u[n][j]+ alpha * ( u[n+1][j] +u[n-1][j]) u1 =  u2 =  u3 =  u4 =  x1 =  for n in range(0,N): u1.append(u[n]) u2.append(u[n]) u3.append(u[n]) u4.append(u[n]) x1.append(x[n]) fig = plt.figure(0) plt.plot(x1,u1) plt.plot(x1,u2) plt.plot(x1,u3) plt.plot(x1,u4) plt.show()