为了在 C++ 中编写 DEL2 matlab 函数,我需要了解算法。我已经设法为不在边界或边缘上的矩阵元素编写函数。我已经看过几个关于它的主题,并通过键入“edit del2”或“type del2”来阅读 MATLAB 代码,但我不明白为获得边界和边缘而进行的计算。
任何帮助将不胜感激,谢谢。
您想近似 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 的索引。
MatLab 中的 DEL2 代表离散拉普拉斯算子,您可以在此处找到有关它的一些信息。
边的主要内容是矩阵内部的元素有四个邻居,而边和角上的元素分别有三个或两个邻居。因此,您以相同的方式计算角和边,但使用的元素更少。
它是如此艰苦!我浪费了几个小时来理解和用 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
这是我在 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