计算欧拉法的相对误差

Question

我使用正演欧拉法求解膜演化的ODE。我有两个向量，X和Y，存储x和y值resp，并在for循环运行10^6中使用该方法。现在，我想计算相对误差，并在误差小于选定值时停止循环。相对误差给出为

|r^t - r^(t-1)|/|r^t|

其中r^t是由所有x和y值组成的向量，在某时刻t。这是我的解

import numpy as np  
x0 = 2*np.pi*5/360/2 #resting membrane lengths
phi = np.linspace(0,2*np.pi, num=360, dtype=float)
X = (5 + 0.5*np.sin(20*phi))*np.cos(phi)
Y = (5 + 0.5*np.sin(20*phi))*np.sin(phi)
L = np.linspace(-1,358, num=360, dtype=int)
R = np.linspace(1,360, num=360,dtype=int) #right and left indexing vectors
R[359] = 0 #closed contour
ds = 1/360
r = np.hstack((X,Y)) #vector of X and Y
rt = np.zeros((10**5,720)) #r at time t
i = 0 
RE = 1 #relative error
e = 10**-10
while RE > e:

lengths = np.sqrt( (X[R]-X)**2 + (Y[R]-Y)**2 )

Ex = (1/10)/ds**2*(X[R] - 2*X + X[L] - x0*( (X[R]-X)/lengths - (X-X[L])/lengths[L]) )
Ey = (1/10)/ds**2*(Y[R] - 2*Y + Y[L] - x0*( (Y[R]-Y)/lengths - (Y-Y[L])/lengths[L]) )

X = X + e*Ex
Y = Y + e*Ey

r = np.hstack((X,Y))
i = i + 1

if i%10000000 ==0: #chooses time values to be multiples of 10^6
rt[i] = r
RE = np.sqrt(np.sum(np.power(rt[i]-rt[i-1],2)))/np.sqrt(np.sum(np.power(rt[i],2)))
print("RE=", RE)
print(i)

然而，循环永远不会结束。我得到这个错误

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

Answer 1

永远运行!问题是你正在运行一个非常耗时的循环:在我的笔记本电脑上，10^7次迭代将运行超过两个小时。

这个错误是因为你没有正确处理rt数组的索引。您只分配10^5行，但是当i是10^7的倍数(而不是10^6，如您的评论所述)时，您只填充行。这就是下标越界的原因:第一次尝试使用它时，它太大了。

也许您只需要在重新计算RE时才需要增加i，而不是每次迭代:

if i%10000000 ==0: #chooses time values to be multiples of 10^6
i = i + 1
rt[i] = r
RE = np.sqrt(np.sum(np.power(rt[i]-rt[i-1],2))) /     
np.sqrt(np.sum(np.power(rt[i],2)))
print("RE=", RE)
print(i)

现在，仍然必须处理初始条件。您最初设置它的方式是在后面增加i一行零。这个保证你的除法会失败，因为分母是0。在重新计算之前，必须对第一行进行编码。