derivative of an array

Question

Hi 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$ and $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?

Answer 1

Purely in terms of terminology, it's probably better to talk about taking discrete partial derivatives of variable fields stored in an array, rather than differentiating an array itself.

Regarding the code itself, you appear to have dropped the element assignment inside your loop,

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)

although that may be a typo. Your version is also bugged and unsafe at the end points of the ranges (i.e i=-nx, i=nx, j=ny, j=-ny) since in these places you are accessing elements of the u array which do not exist. At these locations you either need to use a non centred derivative, something like v(-nx,j) = (u(-nx+1,j) - u(-nx,j)/h, or use any knowledge of the boundary conditions.

As a side note, if you are writing Fortran code frequently, then in terms of code optimisations, it is useful to remember that the index furthest to right should be outermost in loops for the code to have the best chance of running fast, which implies swapping your i and j loops. In many simple cases it's possible the compiler will fix this for you, and it's certainly less important than having code which works, but it's probably a good habit to get into.

Answer 2

It's called column major due to the Fortran array model. Brieftly saying about this: Let $a$ be a matrix of size $m \times n$. Let $a_{ij}$ be element of $a$ for $i=\overline{1,m}$ and $j=\overline{1,n}$, then $a_{ij} $ is stored at the $(i+(j-1)m)$-th position of one-dimenstional storage $A(:)$, i.e. $ a_{ij} = A(i+(j-1)m )$. Writting the loop as the orginal post, the storage $A(:)$ is processed with $m$-increment. However, swapping loop order, the storage $A(:)$ is processed continuously. Hence this approach is performed faster than the first one.