np.polyfit 不会绘制特征，但会给出值

Question

我遇到的问题是，当我使用下面的代码查找错误的 norm-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]))

Answer 1

我对下面代码的演变给出了一条斜率接近 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]))