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 '14 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 '14 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 '14 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 '14 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 '16 at 12:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.