计算科学 - 在python中用后向差分近似计算一阶导数 - 吾爱随笔录

我正在尝试编写一个函数来计算具有反向差分近似的一阶导数。对于第一点，我使用与我刚刚找到的近似的前向差分近似。 $u'(x_i) = \frac{u(x_i) - u(x_i - \Delta_x)}{\Delta x} \equiv D_- u(x_i).$

我的问题是函数不准确，以至于我无法通过测试用例。我收到错误消息

AssertionError: 
Not equal to tolerance rtol=1e-07, atol=0.001

Mismatched elements: 150 / 150 (100%)
Max absolute difference: 0.0585835752066481
Max relative difference: 0.0184951666851795
 x: array([0.519974, 0.519974, 0.538915, 0.557941, 0.577058, 0.596266,
       0.61557 , 0.634974, 0.654479, 0.67409 , 0.69381 , 0.713641,
       0.733589, 0.753655, 0.773844, 0.794158, 0.814602, 0.835178,...
 y: array([0.510532, 0.52943 , 0.548413, 0.567484, 0.586646, 0.605902,
       0.625255, 0.644709, 0.664267, 0.683931, 0.703707, 0.723596,
       0.743602, 0.763729, 0.78398 , 0.804358, 0.824868, 0.845512,...

所以我的函数本身可以工作，只是由于某种原因超出了容差。有没有逻辑错误？我实在想不通我们自己如何才能取得进步。

请帮助解决这个问题。

太感谢了

def compute_prime(x, f):
    """Compute the first derivative"""
   
    N = 150
    x_hat = numpy.linspace(x[0], x[-1], N)
    delta_x = x_hat[1] - x_hat[0]
    
    f_prime_hat = numpy.empty(x_hat.shape)
    for i in range(1, N):
        f_prime_hat[i] = (f(x_hat[i]) - f(x_hat[i-1]))/delta_x
    #f_prime_hat[0] = (f(x_hat[i+1])-f(x_hat[i]))/delta_x
    
    f_prime_hat[0] = (f(x_hat[1]) - f(x_hat[0])) / delta_x
    
    return f_prime_hat

测试用例是

f = lambda x: x**3 / numpy.sin(x)
f_prime = lambda x: -x**3 * numpy.cos(x) / numpy.sin(x)**2 + 3.0 * x**2 / numpy.sin(x)
x = numpy.linspace(0.25, 0.5 * numpy.pi, 150)
numpy.testing.assert_allclose(compute_prime(x, f), f_prime(x), atol=1e-3)
print("Success!")

二阶法

我正在尝试编写一个函数来计算具有反向差分近似的二阶一阶导数。 $u'(x_n) = \frac{3u(x_n) - 4 u(x_{n-1}) + u(x_{n-2})}{2 \Delta x}$

对于第一点，我使用前向差分近似反映了我刚刚找到的那个。 $u'(x_n) = \frac{-2u(x_n+2) + 4 u(x_{n+1}) - 3u(x_{n})}{2 \Delta x}$

我的问题是函数不准确，以至于我无法通过测试用例。

我收到错误消息

AssertionError: 
Not equal to tolerance rtol=1e-07, atol=0.001
Mismatched elements: 150 / 150 (100%)

Max absolute difference: 7.401621709928578
Max relative difference: 1.004627870636921

 x: array([-2.362674e-03,  1.695963e-02,  4.309081e-05,  4.458928e-05,
        4.609489e-05,  4.760789e-05,  4.912854e-05,  5.065709e-05,
        5.219379e-05,  5.373891e-05,  5.529271e-05,  5.685546e-05,...

 y: array([0.510532, 0.52943 , 0.548413, 0.567484, 0.586646, 0.605902,
       0.625255, 0.644709, 0.664267, 0.683931, 0.703707, 0.723596,
       0.743602, 0.763729, 0.78398 , 0.804358, 0.824868, 0.845512,...

我认为我如何找到第一点是不正确的，但真的不确定。请指出我在哪里犯了错误。

请帮助解决这个问题。

太感谢了

def compute_prime(x, f):
    """Compute the first derivative"""
    N = 150
    x_hat = numpy.linspace(x[0], x[-1], N)
    delta_x = x_hat[1] - x_hat[0]

    f_prime_hat = numpy.empty(x_hat.shape)
    for i in range(1, N):
        f_prime_hat[i] = (3*f(x_hat[i]) - 4*f(x_hat[i-1]) + f(x_hat[i-2]))/2*delta_x

    #f_prime_hat[0] = (-f(x_hat[i+2])-4*f(x_hat[i+1])- 3*f(x_hat[i]))/2*delta_x
    f_prime_hat[0] = (-f(x_hat[2])-4*f(x_hat[1])- 3*f(x_hat[0]))/2*delta_x

    return f_prime_hat

测试用例是

f = lambda x: x**3 / numpy.sin(x)

f_prime = lambda x: -x**3 * numpy.cos(x) / numpy.sin(x)**2 + 3.0 * x**2 / numpy.sin(x)

x = numpy.linspace(0.25, 0.5 * numpy.pi, 150)
numpy.testing.assert_allclose(compute_prime(x, f), f_prime(x), atol=1e-3)

print("Success!")