CODEWITHSUNDEEP
arraysfortrandifferential-equations
Jeff Faraci
Jeff Faraci

Reputation: 413

Finite difference derivative of an array

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?

Upvotes: 0

Views: 1266

Answers (1)

Vladimir F Героям слава
Vladimir F Героям слава

Reputation: 60088

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

Upvotes: 3

Related Questions