对于大小为 400 的轴 0，索引 400 超出范围 - Python 中的有限差分方法

问题描述

我明白为什么会出现这个错误。我正在使用有限差分来找到大小为 n = 400 的 numpy 数组的解决方案，并且有问题的术语产生了一个术语 i + 2，其中 for i in range(1,n - 1) 变成了 n + 1。但是，我正在努力寻找一个解决方案。

代码：

import numpy as np

import matplotlib.pyplot as plt

import copy

import sys



# Parameter Values

n = 400

r_max = 10

dr = r_max / n

r_min = dr

r = np.linspace(r_min,r_max,num=n)

#print(dr,r[1]-r[0])



# Imaginary Time Method



# Initial Variable Arrays

g0 = np.zeros(r.size)

fDerR = np.zeros(r.size)

fDerRR = np.zeros(r.size)

gDerR = np.zeros(r.size)

gDerRR = np.zeros(r.size)



k1 = np.zeros(r.size)

k2 = np.zeros(r.size)

k3 = np.zeros(r.size)

k4 = np.zeros(r.size)



j1 = np.zeros(r.size)

j2 = np.zeros(r.size)

j3 = np.zeros(r.size)

j4 = np.zeros(r.size)



timesteps = 5000

dt = 0.0005 #0.0005#0.004



deltaPlot = 1000 #1000#1600#200





# Winding Number

s = 3

# Interaction Strength

U_0 = 1

V_0 = 1

# Impurity Parameters

E = 1

N_I=1.





# Initial Conditions
# Pade' Approximation 

a1=11./32.

a2=11./384.

b1=1./3.

b2=a2

f = np.sqrt((np.square(r)*(a1+a2*np.square(r)))/(1.+b1*np.square(r)+b2*np.power(r,4)))


# Renormalised Initial Condition for Impurity so the mass is M

g = np.exp(-np.square(r)/(2.*np.square(0.5)))

g=np.sqrt(N_I/(2.*np.pi*np.trapz(np.multiply(np.square(g),r))))*g





# Preprare the plots

fig,(ax1,ax2) = plt.subplots(1,2,figsize=(14,5))



label='tau=%f' % (0*dt)

ax1.plot(r,f,label='$\u03C4$=0')
ax1.set_xlabel('r',fontsize=16)

ax1.set_ylabel('f',fontsize=16)

ax1.legend()

#title = 's=%d,V_0=%f' %(s,V_0)

#ax1.set_title(title)

ax2.plot(r,g,label='$\u03C4$=0')

#masstitle='N_I=%f' % (N_I)

#ax2.set_title(masstitle)

ax2.set_xlabel('r',fontsize=16)

ax2.set_ylabel('g',fontsize=16)

#title = 'U_0=%f,E=%f' %(U_0,E)

#ax2.set_title(title)

#ax2.legend()

#plt.show()





# Method Algorithm



for l in range(1,timesteps+1):



    # Compute k1

    tempf = np.copy(f)

    tempg = np.copy(g)



    fDerR[0] = -0.5 * (tempf[2] - 0.) / (2 * dr) + 2 * (tempf[1] - 0.) / (2 * dr) # as it must be f(r=0)=0

    gDerR[0] = -0.5 * (tempg[2] - 0. ) / (dr) + 2 * (tempg[1] - tempg[0]) / dr

    for i in range(1,n-1):

        fDerR[i] = 2 * (tempf[i+1] - tempf[i-1]) / (3 * dr) + tempf[i - 2] / 12 * dr - \
                  tempf[i + 2] / 12 * dr

        gDerR[i] = 2 * (tempg[i+1] - tempg[i-1]) / (3 * dr) + tempg[i - 2] / 12 * dr - \
                   tempg[i + 2] / 12 * dr


    fDerR[n-1] = 1.5 * tempf[n-1] / dr - 2 * tempf[n-2] / dr + 0.5 * tempf[n-3] / dr

    gDerR[n-1] = 1.5 * tempg[n-1] / dr - 2 * tempg[n-2] / dr + 0.5 * tempg[n-3] / dr



    fDerRR[0] = - fDerR[3] / dr + 4 * fDerR[2] / dr + (- 5 * fDerR[1] - 2 * fDerR[0]) / dr

    gDerRR[0] = - gDerR[3] / dr + 4 * gDerR[2] / dr + (- 5 * gDerR[1] - 2 * gDerR[0]) / dr

    for i in range(1,n-1):

        fDerRR[i] = 4 * (fDerR[i+1] + fDerR[i-1]) / (3 * dr) - 2.5 * fDerR[i] / dr - \
                    + fDerR[i-2] / 12 * dr - fDerR[i+2] / 12 * dr

        gDerRR[i] = 4 * (gDerR[i+1] + gDerR[i-1]) / (3 * dr) - 2.5 * gDerR[i] / dr - \
                gDerR[i - 2] / 12 * dr - gDerR[i+2] / 12 * dr

    fDerRR[n-1] = (-5 * fDerR[n-1] + 4 * fDerR[n-2] - fDerR[n-3]) / dr

    gDerRR[n-1] = (-5 * gDerR[n-1] + 4 * gDerR[n-2] - gDerR[n-3]) / dr



    k1=0.5*(fDerRR+fDerR/r-np.square(s/r)*tempf)+(E-V_0*np.square(tempf)-U_0*np.square(tempg))*tempf

    j1=0.5*(gDerRR+gDerR/r)-(U_0*np.square(tempf)-E)*tempg

我选择使用更高的精度（前向和后向为 2，中央为 4）有限差分系数 (https://en.wikipedia.org/wiki/Finite_difference_coefficient)，因为当我为 s = 3 运行此代码时，标准有限差分存在溢出问题由于（很可能）1/r 项的算法。

正如预期的那样，错误发生在

fDerR[i] = 2 * (tempf[i+1] - tempf[i-1]) / (3 * dr) + tempf[i - 2] / 12 * dr - \
                  tempf[i + 2] / 12 * dr.

我的第一个解决方案是对 tempf[i + 2] 项使用 np.append() 但是这在同一行产生了错误 ValueError: setting an array element with a sequence。自从我使用以来，我可能没有正确地实现它：

np.append(tempf,tempf[i + 2] / 12 * dr).

我目前的想法是定义一个大小为 401 的新数组，添加前一个数组中剩余的 400 项，然后添加第 (i + 2) 项，尽管我不确定如何实现它。

解决方法

暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！

如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。

小编邮箱:dio#foxmail.com (将#修改为@）