计算欧拉法的相对误差答案

【问题标题】：Calculating the relative error of a Euler method计算欧拉法的相对误差
【发布时间】：2021-02-19 05:30:27
【问题描述】：

我正在使用正向 Euler 方法来求解膜演化的 ODE。我有两个向量，X 和 Y，分别存储 x 和 y 值，并在 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

【问题讨论】：

请提供预期的minimal, reproducible example (MRE)。我们应该能够复制和粘贴您的代码的连续块，执行该文件，并重现您的问题以及跟踪问题点的输出。这让我们可以根据您的测试数据和所需的输出来测试我们的建议。您发布的代码死于未定义的符号。您还没有向我们展示输出跟踪。
我已经更改了代码，所以它现在应该可以工作了。我还添加了我得到的错误

标签： python numpy

【解决方案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。您必须在重新计算之前对第一行进行编码。

【讨论】：