0
$\begingroup$

I am trying to write a function to compute 1st derivative with backward difference approximation. $ u'(x_i) = \frac{u(x_i) - u(x_i - \Delta_x)}{\Delta x} \equiv D_- u(x_i).$ And for the first point, I use the forward difference approximation that mirrors the one I just found.

My problem is the function is not accurate, so that I cannot pass the test case. I got error message

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,...

So my function itself works, just for some reason the error beyond the tolerance. Is there any logic fault? I really cannot figure our by myself how can I make progress.

Please help with this.

Thank you so much

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

And the test case is

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!")

Second order method

I am trying to write a function to compute 2nd order 1st derivative with backward difference approximation. $u'(x_n) = \frac{3u(x_n) - 4 u(x_{n-1}) + u(x_{n-2})}{2 \Delta x}$

And for the first point, I use the forward difference approximation $u'(x_n) = \frac{-2u(x_n+2) + 4 u(x_{n+1}) - 3u(x_{n})}{2 \Delta x} $ that mirrors the one I just found.

My problem is the function is not accurate, so that I cannot pass the test case.

I got error message

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,...

I think how I find the first point is not correct, but really not sure. Please point it out where did I make mistakes.

Please help with this.

Thank you so much

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

And the test case is

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!")
$\endgroup$
7
  • $\begingroup$ First write down the mathematical functions you want to express, generally “fix my code” posts here does not get good responses. Then at least explain a little about the error you are getting from the code. $\endgroup$
    – boyfarrell
    Jan 16 at 13:20
  • $\begingroup$ Thank you for letting me know. I do not mean to just 'throw' a problem. It's just the math function is not complicate. But I am really new here, thank you so much for pointing it out. I will pay more attention on how to behave properly. $\endgroup$
    – Jingyi Liu
    Jan 16 at 14:56
  • $\begingroup$ First of all Jingyi, welcome to scicomp.stackexchange! This is a great place to ask a very knowledgable group of people all sorts of questions related to computational science. Given that you are a student in my class, and of course this is my homework question, I feel a particular obligation not only to respond but also to help to perhaps provide some guidelines to better form a question for this stack exchange. Rest assured you are perfect welcome to use this resource to ask questions. For the future here are some suggestions for asking homework questions on here: ### Suggestions on Asking H $\endgroup$ Jan 17 at 0:29
  • 1
    $\begingroup$ @JingyiLiu The key missing element here is second order accuracy, which you also left out in your description above. Your finite difference approximation above is only first order accurate. I would suggest editing your question at least to provide the second order accurate analytical version of the finite difference formulas and the code. $\endgroup$ Jan 17 at 0:31
  • $\begingroup$ Thank you for pointing it out. I now changed to 2nd order accurate. $\endgroup$
    – Jingyi Liu
    Jan 18 at 15:19
3
$\begingroup$

You did nothing wrong, the numbers just are what they are. Plotting the error for the backwards difference gives the plot

enter image description here

which is as expected for a method of order one with step size about $h=0.01$. For the second order methods one expects an error of the magnitude $h^2=10^{-4}$, the plot below confirms this

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.