我遇到的问题是，当我使用以下代码查找错误的规范-1时。首先，当我针对步长h绘制误差时，误差值非常小，范围在10 ^ -14到10 ^ -16之间。其次，在下面，您可以看到我将np.polyfit应用于图形的尝试，该图形在运行时将不符合特征，但将输出值。p [0]的值并不完美，因此我认为有些错误，但是它与3的期望输出“接近”。这仅仅是输入错误还是数据错误？

def rk3(A,bvector,y0,interval,N):
    x0=interval[0]
    x_end=interval[1]
    x=np.linspace(x0,x_end,N+1)
    h=(x_end-x0)/N
    y=np.zeros((N+1,len(y0)))
    y[0, :] = y0
    for n in range(N):
        y_1=y[n,:]+h*(np.dot(A,y[n,:])+bvector(x[n]))
        y_2=(3/4)*y[n,:]+(1/4)*y_1+(1/4)*h*(np.dot(A,y_1)+bvector(x[n]+h))
        y[n+1,:]=(1/3)*y[n,:]+(2/3)*y_2+(2/3)*h*(np.dot(A,y_2)+bvector(x[n]+(1/2)*h))
    return x,y

err_vals = []
h_vals = []
for k in range(2,11): #for the range of N=40k, where k=1,...,10
    N=40*k
    x, y = rk3(A,bvector,y0,[0,0.1],N)
    yc = y[-1,:]
    h = (x[-1]-x[0])/N 
    h_vals.append(h)
    yvals.append(yc)
    yn = y[:,1]
    abs_err = np.zeros(N) 
    print("The value of y at k=",k," is ",yc)
for j in range(1,N):
    y_exact=np.array([np.exp(-1000*x[j]), (1000/999)*(np.exp(-x[j])-np.exp(-1000*x[j]))])
    y_exact_2 = y_exact[1]
    abs_err[j] = np.abs((y[j, 1] - y_exact_2)/y_exact_2)
Error = h*np.sum(abs_err[j])
err_vals.append(Error)

p = np.polyfit(np.log(h_vals), np.log(err_vals), 1)
pyplot.loglog(h_vals,err_vals,"kx")
pyplot.xlabel("h")
pyplot.ylabel("Error")
pyplot.loglog(h,np.exp(p[1])*h**(p[0]), 'r--')
print("Best fit line slope ",format(p[0]))

解决方案

我在下面的代码演变过程中给出了一条完整的直线，斜率接近3，用于区间[0,0.01]的积分。

对于给定的间隔[0,0.1]，斜率值大约大1/3。误差分布图（即绝对误差除以步长的预期全局误差幂）给出了收敛模式，确认了该方法的3阶收敛性。

误差2e7*h^3范围相当大，这说明为什么对于较大的步长，问题和方法的组合会变得非常成问题。

误差是通过函数差的L1准则和精确解来计算的，

Error = sum(abs((y-y_exact(x))[:,1]))/sum(abs(y[:,1]))

给出数学上合理的数量。局部相对误差之和可能导致总误差的失真，而精确解的根或值很小。但是，即使使用积分局部相对误差的计算方法，也忽略了第一个数据点零，

Error = sum(abs((y[1:,1]/y_exact(x)[1:,1]-1)))*h

给出类似的线性图，范围向下移至1e-7..1e-9，斜率保持在3.0293

请注意，如果要h_vals在类似计算的列表中使用列表来绘制拟合线，则必须先将其转换为numpy数组。

h=np.asarray(h_vals)

完整的代码

def rk3(A,bvector,y0,interval,N):
    """Solves an IVP y'=f(x, y(x)) on x \in [0, x_end] with y(0) = y0 using N points, using Runge-Kutta method."""    
    x=np.linspace(*interval,N+1)
    h=x[1]-x[0]
    y=np.zeros((N+1,len(y0)))
    y[0, :] = y0
    for n in range(N):
        y_1=y[n]+h*(np.dot(A,y[n])+bvector(x[n]))
        y_2=(3/4)*y[n,:]+(1/4)*y_1+(1/4)*h*(np.dot(A,y_1)+bvector(x[n]+h))
        y[n+1]=(1/3)*y[n]+(2/3)*y_2+(2/3)*h*(np.dot(A,y_2)+bvector(x[n]+0.5*h))
    return x,y

A = np.array([[-1000.0,0.0],[1000.0,-1.0]]);
bvector = lambda x: 0
y_exact = lambda x: np.array([np.exp(-1000*x), (1000/999)*(np.exp(-x)-np.exp(-1000*x))]).T
y0 = y_exact(0)

plt.figure(figsize=(6,3));
h_vals, y_vals, err_vals = [],[],[] 
for k in range(2,11): #for the range of N=40k, where k=1,...,10
    N=40*k
    x, y = rk3(A,bvector,y0,[0,0.01],N)
    yc = y[-1,:]
    h = x[1]-x[0]; 
    plt.plot(x,(y-y_exact(x))[:,1]/h**3)
    h_vals.append(h)
    y_vals.append(yc)
    yn = y[:,1]
    print("The value of y at k=",k," is ",yc)
    Error = sum(abs((y-y_exact(x))[:,1]))/sum(abs(y[:,1]))
    err_vals.append(Error)
plt.grid(); plt.show()

p = np.polyfit(np.log(h_vals), np.log(err_vals), 1)

plt.figure(figsize=(6,4))
plt.loglog(h_vals,err_vals,"kx")
h=np.asarray(h_vals)
plt.plot(h,np.exp(p[1])*h**(p[0]), '--r', lw=0.5)
plt.xlabel("h")
plt.ylabel("Error")
plt.grid(); plt.show()
print("Best fit line slope ",format(p[0]))