2

为了在 C++ 中编写 DEL2 matlab 函数，我需要了解算法。我已经设法为不在边界或边缘上的矩阵元素编写函数。我已经看过几个关于它的主题，并通过键入“edit del2”或“type del2”来阅读 MATLAB 代码，但我不明白为获得边界和边缘而进行的计算。

任何帮助将不胜感激，谢谢。

4

4 回答 4

5

您想近似 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

编辑说明：

拉普拉斯算子是 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”替换为“ik”） . 这同样适用于 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 是二阶离散拉普拉斯算子，即它允许在两个相邻节点之间的距离为h的方形笛卡尔网格 NxN 上给定实际连续函数的拉普拉斯算子近似值。

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^矩阵 A 中的 2 个值。

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

如果您的起始函数是线性或二次函数（x^2+y^2 很好，x^3 + y^3 不好），del2 将返回连续拉普拉斯算子的确切值。如果函数不是线性的也不是二次的，使用的点越多，结果就越准确（即在极限 h -> 0 内）

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

于 2012-11-12T11:36:33.540 回答
1

MatLab 中的 DEL2 代表离散拉普拉斯算子，您可以在此处找到有关它的一些信息。

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

于 2012-11-12T11:35:47.940 回答
0

它是如此艰苦！我浪费了几个小时来理解和用 Java 实现它。

这里是：https ://gist.github.com/emersonmoretto/dec8f7125c032775da0d

测试并与原始函数 DEL2 (Matlab) 进行比较

我在 sbabbi 回复中发现了一个错字：

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(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+1,j) + 4 u(i+2,j) - u(i+3,j) + 2 u(i,j)]/h^2
于 2014-09-03T04:34:22.723 回答
0

这是我在 Fortran 90 中编写的一个模块，它复制了 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
于 2013-08-02T00:39:29.883 回答