不确定为什么会出现此错误，我有一个想法，但不确定如何解决

Question

这是我正在从事的项目。我正在尝试在地壳内的某些主岩或乡村岩石中模拟岩浆冷却的侵入堤坝。我对编码很陌生。我尽力将此代码从另一种编码语言转换为 python。我对正在发生的事情有一个基本的了解。我知道我正在尝试索引超出范围的内容，但不确定在哪里以及如何修复它。我能得到的任何帮助将不胜感激！提前致谢。

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 的形状是 (1,nx)，即它是一个 1×nx 二维数组。

如果您只需要一维数组，请改为设置 Tnew = np.zeros(nx)。或者，如果您想保持二维，请访问 Tnew[0,i]。

Answer 2

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

This will gives you one array of nx element inside array (Two dimensional array) and you need to access it like Tnew[0][i]

just change it to

Tnew = np.zeros((nx))

检查numpy.zeros DOC

注意：解决此错误后，您将面临 "T[i+1]" 处的索引错误，这是因为当您的 'i' 到达最后一个元素时，您不会在 'T'.

中获得 'i+1' 元素

不确定为什么会出现此错误，我有一个想法，但不确定如何解决

Not sure why I am getting this error, I have an idea but am not sure how to fix it

python

numerical-methods