1
$\begingroup$

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?

$\endgroup$
5
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Hi the answer you gave is very clear, thanks. I just have a question about swapping my i,j loops. Why is it that the code will run faster if the j loop is outermost? I would like to learn more about this, could you perhaps point me in a direction where I could read about this? I have never heard this before, and its very helpful!!! Thanks again. $\endgroup$ – Jeff Faraci Oct 22 '16 at 20:07
  • $\begingroup$ @integrals tqviet's answer expresses the basic idea. You could also try looking at the wikipedia page for the C style alternative en.wikipedia.org/wiki/Row-major_order and the links it contains. But again, this isn't something to worry too much about until you're writing big, complicated programs. $\endgroup$ – origimbo Oct 23 '16 at 19:21
2
$\begingroup$

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.

$\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.