Je ne sais pas pourquoi je reçois cette erreur, j'ai une idée, mais je ne suis pas sûr de savoir comment y remédier

Question

Ceci est pour un projet sur lequel je travaille. J'essaie de modéliser une digue intrusive de refroidissement du magma dans un rocher de l'hôte ou du pays dans la croûte terrestre. Je suis relativement nouveau au codage. J'ai fait de mon mieux pour convertir ce code en un autre langage de codage en python. J'ai une idée de base de ce qui se passe. Je sais que j'essaie d'indexer quelque chose hors de portée, mais je ne suis pas sûr de savoir où et comment y remédier. Toute aide que je pourrais obtenir serait grandement appréciée! Merci d'avance.

import numpy as np 
import matplotlib.pyplot as plt 

#Physical parameters 
L = 100 #Length of modeled domain [m] 
Tmagma = 1200 #Temp. magma [°C] 
Trock = 300 #Temp. of country rock [°C] 
kappa = 1e-6 #Thermal diffusivity of rock [m^2/s] 
W = 5 #Width of dike [m] 
day = 3600*24 #seconds per day 
dt = 1*day #Timestep 
print(kappa) 
print(day) 
print(dt) 



#Numerical parameters 
nx = 201 #Number of gridpoints in x-direction 
nt = 500 #Number of timesteps to compute 
dx = L/(nx-1) #Spacing of grid 
x = np.linspace(-50,50,100) #Grid 
print(dx) 
print(x) 

#Setup initial temp. profile 
T = np.ones(np.shape(x))*Trock 
T[x>=-W/2] = 1200 
T[x>=W/2] = 300 
print(T) 


time = 0 
for n in range(0,nt): #Timestep loop 
    #compute new temp. 
    Tnew = np.zeros((1,nx)) 
    print(Tnew) 
    for i in range(2,nx-1): 
     Tnew[i] = T[i]+kappa*dt*((T[i+1]-2*T[i]+T[i-1])/(dt**2))    # Set BC's 
    Tnew[1] = T[1] 
    Tnew[nx] = T[nx] 

    #update temp. and time 
    T = Tnew 
    time = time+dt 

    #Plot solution 
    plt.figure() 
    plt.plot(x,Tnew) 
    plt.xlabel('x [m]') 
    plt.ylabel('Temerpature [°C]') 
    plt.legend() 
    surf = ax.plot_surface(X, Y, cmap=cm.coolwarm, linewidth=0, antialiased=False) 
    fig.colorbar(surf, shrink=0.5, aspect=5) 
    plt.show() 

---

IndexError        Traceback (most recent call last) 
<ipython-input-51-e80d6234a5b4> in <module>() 
    37  print(Tnew) 
    38  for i in range(2,nx-1): 
---> 39   Tnew[i] = T[i]+kappa*dt*((T[i+1]-2*T[i]+T[i-1])/(dt**2)) 
    40 
    41  # Set BC's 

IndexError: index 2 is out of bounds for axis 0 with size 

Answer 1

Tnew = np.zeros((1,nx)) 

Cela vous donne un réseau d'élément nx intérieur tableau (Deux tableau dimensions) et vous devez y accéder comme Tnouv [0] [i]

changer juste à

Tnew = np.zeros((nx)) 

check numpy.zeros DOC

NO TE: Après avoir résolu cette erreur, vous allez faire face à l'erreur d'index "T [i + 1]", car lorsque votre "i" atteint le dernier élément, vous n'obtiendrez pas l'élément 'i + 1' dans ' T '.

Answer 2

La forme que vous avez donné Tnew est (1,nx), à savoir qu'il est un tableau 2D 1 sous-nx.

Si vous souhaitez uniquement une baie 1D, définissez Tnew = np.zeros(nx) à la place. Ou si vous voulez le garder 2D, accédez au Tnew[0,i].