Finite difference derivative of an array

Question

I am trying to take a derivative of an array but am having trouble. The array is two dimensional, x and y directions. I would like to take a derivative along x and along y using central difference discretization. The array has random values of numbers, no values are NaN. I will provide a basic portion of the code below to illustrate my point (assume the array u is defined and has some initial values already inputted into it)

integer :: i,j
integer, parameter :: nx=10, ny=10
real, dimension(-nx:nx, -ny:ny) :: u,v,w
real, parameter :: h

do i=-nx,nx
  do j=-ny,ny

    v = (u(i+1,j)-u(i-1,j))/(2*h)
    w = (u(i,j+1)-u(i,j-1))/(2*h)

  end do 
end do

Note, assume the array u is defined and filled up before I find v,w. v,w are supposed to be derivatives of the array u along x and along y,respectively. Is this the correct way to take a derivative of an array?

I am afraid you need to post your question on one of the stackexchange forums closer to numerical methods. Stackoverflow is more into programming. Just trying to help, not to push you out. Good luck! :-) — Tarik, Oct 21 '16 at 16:06
@HighPerformanceMark The trouble I am having is to implement this chosen method. Thanks for your help — Jeff Faraci, Oct 21 '16 at 16:35
@HighPerformanceMark Okay I will check the boundaries carefully, that makes sense. — Jeff Faraci, Oct 21 '16 at 16:40
Approach negative array indices with additional caution, see: http://stackoverflow.com/questions/38140951/fortran-subroutine-returning-wrong-values/38154225#38154225 — jlokimlin, Oct 21 '16 at 17:12
@agentp Probably yes. But also the loop bounds must be changed, there is an out of bound access on the right hand side. — Vladimir F Героям слава, Oct 21 '16 at 17:56
yes, I meant in addition to the already noted boundary issue. — agentp, Oct 21 '16 at 18:02
In future, please do not cross-post questions. For more information, see [here](http://meta.stackexchange.com/q/64068). — Matt, Oct 23 '16 at 14:56

Vladimir F Героям слава · Accepted Answer · 2016-10-21T22:14:06.973

I can see several problems in your code.

1.You must be careful what you have on the left hand side.

v = (u(i+1,j)-u(i-1,j))/(2*h)

means that the whole array v will be set to the same number everywhere. You don't want this in a loop. In a loop you want to set just one point at a time

v(i,j) = (u(i+1,j)-u(i-1,j)) / (2*h)

and 2) You are accessing the array out of bounds. You can keep the simple loop, but you must use the boundary points as "ghost points" which store the boundary values. If I assume that points -nx,nx,-nyandny` are lying on the boundary, then you can only compute the derivative using the central difference inside the domain:

do i=-nx+1,nx-1
  do j=-ny+1,ny-1

    v(i,j) = (u(i+1,j)-u(i-1,j)) / (2*h)
    w(i,j) = (u(i,j+1)-u(i,j-1)) / (2*h)

  end do 
end do

If you need the derivative on the boundary, you must use a on-sided difference like

  do j=-ny+1,ny-1

    v(nx,j) = (u(nx,j)-u(nx-1,j)) / h
    w(nx,j) = (u(nx,j+1)-u(nx,j-1)) / h

  end do

Thanks a lot. It's all clear, I just have one question. I have a second derivative in my equation that I use a 5 point discretization for, so my boundary conditions have an extra layer, thus the interior of my grid is nx-2,-nx+2, etc. So should I do the derivatives of v(i,j), w(i,j) in loops with i=-nx+2,nx-2, etc? Thanks. — Jeff Faraci, Oct 23 '16 at 03:35

Finite difference derivative of an array

1 Answers1