Понимание функции DEL2 в Matlab для кодирования ее в переполнении стека
Для того, чтобы кодировать функцию MATLAB DEL2 в C ++, мне нужно понять алгоритм. Мне удалось закодировать функцию для элементов матрицы, которые не находятся на границах или краях.
Я видел несколько тем об этом и читал код MATLAB, набирая «edit del2» или «type del2», но я не понимаю, какие вычисления сделаны для получения границ и краев.
Любая помощь будет оценена, спасибо.
Решение
Вы хотите приблизиться к 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’ » + + (h ^ 4)
u (i + 3) = u (i) + 3 h u ‘+ (9/2) h ^ 2 u’ ‘+ (27/6) h ^ 3 u’ » + + (h ^ 4)
Решение для вас дает (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
РЕДАКТИРОВАТЬ уточнения:
Лапласиан — это сумма 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, вы можете использовать первую формулу, которую я написал для вычисления производных (обратите внимание, что из-за симметрии 2-й производной, если 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 — дискретный лапласиан 2-го порядка, т. Е. Он позволяет аппроксимировать лапласиан реальный непрерывная функция, учитывая ее значения на квадратной декартовой сетке NxN, где расстояние между двумя соседними узлами равно час.
h ^ 2 — это просто постоянный размерный фактор, вы можете получить реализацию matlab из этих формул, установив h ^ 2 = 4.
Например, если вы хотите вычислить реальный Лапласиан u (x, y) на квадрате (0, L) x (0, L), что вы делаете, записываете значения этой функции в декартовой сетке 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).
Del2 вернет точное значение непрерывный Laplacian, если ваша начальная функция является линейной или квадратичной (x ^ 2 + y ^ 2 хорошо, x ^ 3 + y ^ 3 не хорошо). Если функция не является ни линейной, ни квадратичной, результат будет тем точнее, чем больше точек вы используете (т. Е. В пределе h -> 0)
Я надеюсь, что это более понятно, обратите внимание, что я использовал индексы на основе 0 для доступа к матрице (стиль массива C / C ++), в то время как в Matlab используется основанный на 1.
Другие решения
DEL2 в MatLab представляет собой дискретный оператор Лапласа, вы можете найти некоторую информацию о нем Вот.
Главное в ребрах состоит в том, что элементы внутри матрицы имеют четырех соседей, в то время как элементы на ребрах и углах имеют трех или двух соседей с уважением. Таким образом, вы рассчитываете углы и ребра одинаково, но с использованием меньшего количества элементов.
Вот модуль, написанный мною в Fortran 90, который реплицирует оператор «del2 ()» в MATLAB, реализуя приведенные выше идеи. Это работает только для массивов, которые по крайней мере 4×4 или больше. Он работает успешно, когда я запускаю его, поэтому я решил опубликовать его, чтобы другие люди не тратили время на создание своего.
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
Это очень трудно! Я потратил несколько часов, чтобы понять и реализовать его на 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
detector