阵列的有限差分导数

Question

我正在尝试一个数组的导数，但有麻烦。该阵列为二维、x和y方向。我想用中心差分离散化，沿着x和y取一个导数。数组有随机数值，没有值是NaN。我将提供下面代码的一个基本部分来说明我的观点(假设数组u已经定义，并且已经输入了一些初始值)

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

注意，假设数组u是在我找到v，w之前定义和填充的。v、w分别是沿x和沿y的阵列u的导数。这是获取数组的导数的正确方法吗？

Answer 1

我可以在您的代码中看到几个问题。

1.你必须小心左手边的东西。

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

这意味着整个数组v都将被设置为相同的数字。你不想把这个循环起来。在循环中，您希望一次只设置一个点。

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

2)您正在访问超出界限的数组。您可以保留简单的循环，但是必须使用边界点作为“鬼点”来存储边界值。如果假设点-nx、nx,-nyandny` `位于边界上，则只能使用域中的中心差分计算导数：

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

如果需要在边界上使用导数，则必须使用如下所示的边差

  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