理解Matlab语言中的DEL2函数，以便在C++中编写

Question

为了在c++中编写DEL2 matlab函数，我需要理解算法。我已经成功地为矩阵中不在边界或边缘上的元素编写了函数。我已经看过关于它的几个主题，并通过输入“编辑del2”或“类型del2”来阅读MATLAB代码，但我不理解为获得边界和边缘而进行的计算。

任何帮助都会很感谢，谢谢。

Answer 1

你只需要知道一个点右边(或左边)的u的值，就可以近似u'‘。为了有一个二阶近似，你需要3个方程(基本泰勒展开)：

u(i+1) = u(i) +h u‘+ (1/2) h^2 u'’+ (1/6) h^3 u‘+ O(h^4)

u(i+2) = u(i) +2 h u‘+ (4/2) h^2 u'’+ (8/6) h^3 u‘’+ O(h^4)

u(i+3) = u(i) +3 h u‘+ (9/2) h^2 u'’+ (27/6) h^3 u‘’+ O(h^4)

求解u'‘可得到(1)

h^2 u'‘= -5 u(i+1) +4 u(i+2) - u(i+3) +2 u(i) +O(h^4)

为了得到拉普拉斯公式，你需要在边界上用这个公式代替传统的公式。

例如，在"i = 0“的情况下，你会得到：

del2(u) (i=0，j) = -5 u(i+1，j) +4 u(i+2，j) - u(i+3，j) +2 u(i，j) + u(i，j+1) + u(i，j-1) - 2u(i，j) /h^2

编辑 clarifications：

拉普拉斯量是x和y方向上的二阶导数之和。可以用公式(2)计算二阶导数

u'‘= (u(i+1) + u(i-1) - 2u(i))/h^2

如果你同时有u(i+1)和u(i-1)。如果是i=0或i=imax，你可以使用我写的第一个公式来计算导数(请注意，由于二阶导数的对称性，如果i= imax，你可以将"i+k“替换为"i-k")。Y (j)方向也是如此：

在边缘，您可以混合使用公式(1)和(2)

del2(u) (i=imax，j) = -5 u(i-1，j) +4 u(i-2，j) - u(i-3，j) +2 u(i，j) + u(i，j+1) + u(i，j-1) - 2u(i，j) /h^2

del2(u) (i，j=0) = -5 u(i，j+1) +4 u(i，j+2) - u(i，j+3) +2 u(i，j) + u(i+1，j) + u(i-1，j) - 2u(i，j) /h^2

del2(u) (i，j=jmax) = -5 u(i，j-1) +4 u(i，j-2) - u(i，j-3) +2 u(i，j) + u(i+1，j) + u(i-1，j) - 2u(i，j) /h^2

在角落里，你只需要在两个方向上使用(1)两次。

del2(u) (i=0，j=0) = -5 u(i，j+1) +4 u(i，j+2) - u(i，j+3) +2 u(i，j) + -5 u(i，j+1) +4 u(i+2，j) - u(i+3，j) +2 u(i，j)/h^2

Del2是二阶离散拉普拉斯算子，即它允许近似实连续函数的拉普拉斯算子，给定它在正方形笛卡尔网格NxN上的值，其中两个相邻节点之间的距离是h。

h^2只是一个恒定的维度因子，你可以通过设置h^2 = 4从这些公式中获得matlab实现。

例如，如果你想在(0，L) x (0，L)正方形上计算u(x，y)的实拉普拉斯矩阵，你要做的就是在NxN笛卡尔网格上写下这个函数的值，即你计算u(0,0)，u(L/(N-1),0)，u(2L/(N-1),0) ...u( (N-1)L/(N-1) =L,0) ...u(0，L/(N-1))，u(L/(N-1)，L/(N-1))等等，你把这N^2个值记在一个矩阵A中。

然后你会得到ans = 4*del2(A)/h^2，其中h= L/(N-1)。

如果你的起始函数是线性的或二次的(x^2+y^2很好，x^3 + y^3不好)，del2将返回连续拉普拉斯的精确值。如果函数既不是线性的，也不是二次的，那么使用的点越多，结果就越准确(即在limit h -> 0中)

我希望这一点更清楚，请注意，我使用基于0的索引来访问矩阵(C/C++数组样式)，而matlab使用基于1的索引。

Answer 2

MatLab中的DEL2代表离散拉普拉斯算子，你可以在here中找到关于它的一些信息。

关于边的主要事情是，矩阵内部的元素有四个邻居，而边和角上的元素分别有三到两个邻居。因此，您以相同的方式计算角和边，但使用的元素更少。

Answer 3

这是我用Fortran90编写的一个模块，它复制了MATLAB中的"del2()“运算符，实现了上述思想。它只适用于至少4x4或更大的数组。当我运行它时，它工作得很成功，所以我想我应该把它贴出来，这样其他人就不用浪费时间制作自己的东西了。

module del2_mod
implicit none
real, private                       :: pi
integer, private                    :: nr, nc, i, j, k
contains
! nr is number of rows in array, while nc is the number of columns in the array.
!!---------------------------------------------------------- 

subroutine del2(in, out)
real, dimension(:,:)            :: in, out
real, dimension(nr,nc)          :: interior, left, right, top, bottom, ul_corner, br_corner, disp
integer                         :: i, j
real                            :: h, ul, ur, bl, br
! Zero out internal arrays
out = 0.0; interior=0.0; left = 0.0;  right = 0.0;  top = 0.0;  bottom = 0.0;  ul_corner = 0.0; br_corner = 0.0;
h=2.0

! Interior Points
do j=1,nc
    do i=1,nr
    ! Interior Point Calculations
    if( j>1 .and. j<nc .and. i>1 .and. i<nr )then
        interior(i,j) = ((in(i-1,j) + in(i+1,j) + in(i,j-1) + in(i,j+1)) - 4*in(i,j) )/(h**2)
    end if
    ! Boundary Conditions for Left and Right edges
    left(i,1) = (-5.0*in(i,2) + 4.0*in(i,3) - in(i,4) + 2.0*in(i,1) + in(i+1,1) + in(i-1,1) - 2.0*in(i,1) )/(h**2)
    right(i,nc) = (-5.0*in(i,nc-1) + 4.0*in(i,nc-2) - in(i,nc-3) + 2.0*in(i,nc) + in(i+1,nc) + in(i-1,nc) - 2.0*in(i,nc) )/(h**2)
    end do
! Boundary Conditions for Top and Bottom edges
top(1,j) = (-5.0*in(2,j) + 4.0*in(3,j) - in(4,j) + 2.0*in(1,j) + in(1,j+1) + in(1,j-1) - 2.0*in(1,j) )/(h**2)
bottom(nr,j) = (-5.0*in(nr-1,j) + 4.0*in(nr-2,j) - in(nr-3,j) + 2.0*in(nr,j) + in(nr,j+1) + in(nr,j-1) - 2.0*in(nr,j) )/(h**2)
end do
out = interior + left + right + top + bottom 
! Calculate BC for the corners
ul = (-5.0*in(1,2) + 4.0*in(1,3) - in(1,4) + 2.0*in(1,1) - 5.0*in(2,1) + 4.0*in(3,1) - in(4,1) + 2.0*in(1,1))/(h**2)
br = (-5.0*in(nr,nc-1) + 4.0*in(nr,nc-2) - in(nr,nc-3) + 2.0*in(nr,nc) - 5.0*in(nr-1,nc) + 4.0*in(nr-2,nc) - in(nr-3,nc) + 2.0*in(nr,nc))/(h**2)
bl = (-5.0*in(nr,2) + 4.0*in(nr,3) - in(nr,4) + 2.0*in(nr,1) - 5.0*in(nr-1,1) + 4.0*in(nr-2,1) - in(nr-3,1) + 2.0*in(nr,1))/(h**2)
ur = (-5.0*in(1,nc-1) + 4.0*in(1,nc-2) - in(1,nc-3) + 2.0*in(1,nc) - 5.0*in(2,nc) + 4.0*in(3,nc) - in(4,nc) + 2.0*in(1,nc))/(h**2)
! Apply BC for the corners
out(1,1)=ul
out(1,nc)=ur
out(nr,1)=bl
out(nr,nc)=br
end subroutine

end module