代码之家 › 专栏 › 技术社区 › PyRsquared

如何在python中计算数值导数的边界点?

numerical-methods math arrays python

PyRsquared · 技术社区 · 7 年前

Central difference formula . 问题是,我无法计算导数的边界点,因为中心差分公式使用的指数超出了界限。我的代码在下面

import numpy as np
n = 20000 # number of points in array
xs = np.linspace(start=-2*np.pi, stop=2*np.pi, num=n) # x values
y = np.array([np.sin(i) for i in xs]) # our function, sine

def deriv(f, h):
    """
    Calauclate the numerical derivative of any function
    :param f: numpy.array(float), the array of numbers we differentiate
    :param h: step size
    :rtype d: numpy.array(float)
    """
    d = np.zeros_like(f)
    # this loop misses the first and last points in f
    for i in range(1, f.shape[0]-1):
        # 2-point formula
        d[i] = (f[i+1] - f[i-1])/(2*h)

    return d

h = abs(xs[0] - xs[1]) # step size
y1 = deriv(y, h) # first derivative
y2 = deriv(y1, h) # second derivative
y3 = deriv(y2, h) # third derivative

y,y1,y2,y3 你可以看到它在终点爆炸

我想做的是把终点定在他们最近的邻居那里 deriv

...
d = np.zeros_like(f)
for i in range(1, f.shape[0]-1):
    d[i] = (f[i+1] - f[i-1])/(2*h)

d[0] = d[1]
d[-1] = d[-2]
...

中间的导数,远离边界,计算精细。问题在于边界。

我应该如何处理这里的边界条件?不同的数值微分格式比中心差分格式更有效吗?

编辑: 我正在寻找一个通用的方法来解决这个问题,而不仅仅是一个方法,可以应用于正弦函数或任何其他周期函数,因为我已经用这里来说明这个问题。

2 回复 | 直到 7 年前

b-fg 7 年前

无论如何,如果你的函数有周期性边界条件(它看起来是一个正弦波,所以在这种情况下你有周期性)只需创建一个新的数组和2个额外的元素:新的数组开始元素将是你的最后一个元素的原始数组和结束元素的新数组将是原始数组的开始元素。这里有一个方法

f_periodic = np.zeros(f.size+2)
f_periodic[1:-1], f_periodic[0], f_periodic[-1] = f, f[-1], f[0]

你现在可以在 f_periodic 为了什么 d[1] 和 d[-2] d[0] 和 d[-1] ) .

extrapolate 新边界的值。根据数字分化的顺序,需要更多的鬼细胞。对于当前方案,可以进行简单的线性外推(每个边界只需要1个幻影值):

f_new = np.zeros(f.size+2)
f_new[1:-1] = f
f_new[-1] = f[-2] + (f[-2]-f[-3])/(x[-2]-x[-3])*(x[-1]-x[-2])
f_new[0] = f[1] + (f[1]-f[2])/(x[1]-x[2])*(x[0]-x[1])

请注意,您还必须扩展 x . 但是,因为你有一个固定的间距,所以使用 h 而不是空间差异,例如。 x[-2]-x[-3] f_new 你将得到边界上导数的一阶近似值(因为你使用了线性外推来寻找幻影值)。

在边界上使用向前和向后方案

我不会在这里显示代码,但基本上您需要分别使用边界值和左右边界的right(forwards)或left(backwards)值进行区分。这是一阶近似。

Lutz Lehmann 7 年前

对于边界点,可以使用2阶正向和反向微分方案。基本上,我们知道

(f(x+h)-f(x))/h = f'(x) + h/2*f''(x) + O(hÂ²)                       (I)

(f(x+2h) - 2f(x+h) + f(x))/hÂ² = f''(x+h) + O(hÂ²) = f''(x) + O(h)   (II)

最后用二阶导数代替一阶项,即计算 (I)-h/2*(II) 得到

(-1/2*f(x+2h) + 2*f(x+h) -3/2*f(x))/h = f'(x) + O(hÂ²)

注意,一阶导数中的O(h)误差通常会导致除差的第二次迭代中的O(h)误差和第三次迭代中的O(1)。有人可能会说,误差项会适当地抵消,但这只会发生在内部点,单边导数会随着与边界距离的增加而“破坏”这种模式。