2017-06-09 29 views
1

这是为我正在处理的项目。我正在试图在地壳内的一些主体或乡村岩石中模拟岩浆冷却的侵入堤防。我对编码相当陌生。我尽我所能将这种代码格式从另一种编码语言转换为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 

回答

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

这将让您内部阵列NX元件(二维阵列)的一个阵列,并且您需要访问它像TNEW [0] [i]于

只是它更改为

Tnew = np.zeros((nx)) 

检查numpy.zeros DOC

NO TE:解决这个错误后,你要面对的是索引错误在“T [i + 1]”,这是因为当你的'i'到达最后一个元素时,你不会在' T”。

0

你给Tnew形状(1,nx),即它是一个1逐nx 2D阵列。

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

相关问题