+关注
已关注

分类  

暂无分类

标签  

暂无标签

日期归档  

暂无数据

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

发布于2021-01-23 19:52     阅读(54)     评论(0)     点赞(18)     收藏(5)


0

1

2

3

4

5

6

我遇到的问题是,当我使用以下代码查找错误的规范-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]))

0

1

2

3

4

5

6

7



所属网站分类: 技术文章 > 问答

作者:黑洞官方问答小能手

链接: https://www.pythonheidong.com/blog/article/787088/10e618854b11c368130d/

来源: python黑洞网

任何形式的转载都请注明出处,如有侵权 一经发现 必将追究其法律责任

18 0
收藏该文
已收藏

评论内容:(最多支持255个字符)