0

我明白为什么会发生这个错误。我正在使用有限差分来找到大小为 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.

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

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

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

4

0 回答 0