Search code examples
arraysfortrannumerical-methodsautomatic-differentiation

How to treat boundary points in numerical differentiation?


I have a query on numerical differentiation that's beyond the language used. Suppose I have an array of n points x and f(x) and I want to take the first derivative of f(x). Every method will consume points making the derivative array shorter than the function, so how can a "lengthen" the array in a smart way. For example I want to take the derivate using the Five-point stencil, that is

    f'(0) = 1/12 h (-f(-2) + 8 f(-1)- 8 f(1) + f(2))

where f(n) is the function evaluated at the n-th point. So with this method the f' array is shorter by 4 points. How can I lengthen this array in a smart way and if it's possible in a way that produces an error similar as this 5 point stencil method?


Solution

  • SUBROUTINE Diff(X, N, XPrime)
    INTEGER           , INTENT(IN   ) :: N
    REAL, DIMENSION(N), INTENT(IN   ) :: X
    REAL, DIMENSION(N), INTENT(INOUT) :: XPrime
    enter code here
    REAL, DIMENSION(-1:N+2)           :: Temp
    INTEGER                           :: I
    
    !Use temp for X
    Temp(1:N) = X
    !... Temp(O)   = X(1) - (X(2) - X(1)  )
    !... Temp(N+1) = X(N) + (X(N) - X(N-1))
    
    !Your code here
    
    !output XPrime from 1:N
    
    END SUBROUTINE Diff
    

    The middle of the vector is easy, just the ends needs something special.

    For X' maybe Temp(0 :N+1).

    For X'' maybe Temp(-1:N+2).

    Of course it should not take long to realise that you could get rid of Temp entirely and do the ends manually. It depends on how long your vectors are, and whether you need some alignment. In some parallel world you may want the temp array as a function, For a simple serial implementation the temp may be conceptually easier to grasp. Also you mentioned lengthening the array, which is really accordioning the vector by grabbing each end and moving your paws further apart. Extending it means that you have chin-scratch to keep track of the index, where as in the above implimentation the X, X', and Temp are all aligned in index value. Since fortran can go from any#:AnyOther# this the perfect example of when you may want to do that.