Python - calculation time derivative and laplacien by finite differences

I would like to determine a temporal derivative and a Laplacian by the finite differences method to solve the heat equation in a 1-dimensional case. My aim is to get the sources term that is why I need to get all left terms of the equation.

When I plot the sources term as a function of time, the result is surprising... very far from expected. Is this a problem with the np.reshape? Thanks for help ;)

#profil1D_theta.shape[0] : time
#profil1D_theta.shape[1] : space

#profil1D_theta_reshape.shape
#Out[986]: (328, 411)

T_tempo = 1./115. #time period

#Time derivative

dt_theta_list = []

for i in range(profil1D_theta.shape[0]-1):
for j in range(profil1D_theta.shape[1]-1):
dt_theta = (profil1D_theta[i+1,j] - profil1D_theta[i,j])/T_tempo
dt_theta_list.append(dt_theta)

dt_theta = np.asarray(dt_theta_list)
dt_theta_reshape = dt_theta.reshape(328,411)

#Laplacian

T_spatial = 0.000021661 #spatial period

dx2test_theta_list =[]
for i in range(0,profil1D_theta.shape[0]-1):
for j in range(0,profil1D_theta.shape[1]-1):
dxtest_theta = (profil1D_theta[i,j+1] - 2*profil1D_theta[i,j] + profil1D_theta[i,j-1])/T_spatial**2
dx2test_list.append(dxtest_theta)

dx2test_theta = np.asarray(dx2test_theta_list)

dx2_theta_reshape = dx2test_theta.reshape(328,411)

rho = 7850.
C = 500.
k = 42.
D = k/(rho*C)
e = 3.*0.001
b = 6.*0.001
Boltz = 5.679373*10**(-8)
Emissivite_echant = 0.97 #
T0 = 273.15 + 26.
h = 5.

tau_th = rho*C*e/(2*h+8*Boltz*Emissivite_echant*T0**3)

#Sources

profil1D_theta_reshape = profil1D_theta[0:328,0:411]

w1_chaleur_list=[]

for i in range(0,profil1D_theta_reshape.shape[1]):
for j in range(0,profil1D_theta_reshape.shape[0]):
w_chaleur1D = rho*C*(dt_theta_reshape[j,i]+1./(tau_th)*profil1D_theta_reshape[j,i] - D*dx_theta_reshape[j,i])
w1_chaleur_list.append(w_chaleur1D)

d,e = dt_theta_reshape.shape
w1_chaleur = np.asarray(w1_chaleur_list)
w1_chaleur_reshape = np.reshape(w1_chaleur,(e,d))

clim_min = 0
clim_max = 1*10**6

plt.imshow(w1_chaleur_reshape);plt.clim([clim_min,clim_max]);plt.colorbar();
plt.ylabel(r'Axe longitudinal')
plt.xlabel(r'Temps (images)')
plt.show()

• You're not likely to get help with such a vague question and only hints of what's wrong. Please add more detail about your results (such as plots and what is "expected") and why you suspect the problem might be with np.reshape. As the question is written, you're basically asking us to just read and debug your code. – Doug Lipinski Jul 17 '14 at 23:34
• You derivatives are edge-centered. This means that the size of the result array is going to be smaller than the input array. Also, you Laplacian code will likely do a bunch of things wrong given the lack of bounds checking. See numpy.diff() for a smarter way to approach this all. – meawoppl Jul 20 '14 at 23:58
• I can add a bit of example code of how to do this if you are interested. – meawoppl Jul 21 '14 at 0:02
• I'm coding my own version of Heat Equation Solver, in 2D. You might want to take a look at how I solved the Laplacian problem for the border of the grid. https://github.com/matpompili/HeatFlow – Matteo Pompili Aug 15 '14 at 22:04
• You should check out my answer to this question about the same problem solving the heat equation. It's written for matlab but it should be easy to port to python. – nluigi Mar 8 '16 at 12:48