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

#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 = np.asarray(dt_theta_list)
dt_theta_reshape = dt_theta.reshape(328,411)


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


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


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

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.ylabel(r'Axe longitudinal')
plt.xlabel(r'Temps (images)')
  • 13
    $\begingroup$ 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. $\endgroup$ Jul 17, 2014 at 23:34
  • $\begingroup$ 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. $\endgroup$
    – meawoppl
    Jul 20, 2014 at 23:58
  • $\begingroup$ I can add a bit of example code of how to do this if you are interested. $\endgroup$
    – meawoppl
    Jul 21, 2014 at 0:02
  • $\begingroup$ 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 $\endgroup$ Aug 15, 2014 at 22:04
  • $\begingroup$ 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. $\endgroup$
    – nluigi
    Mar 8, 2016 at 12:48


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.