Дискретные дифференциальные операторы / forpes.ru

Главная
Дискретные дифференциальные операторы

Дискретные дифференциальные операторы +3

23.11.2025 07:16

Ordevoir 0 267 Источник

Каждый раздел содержит по три подраздела: непрерывный случай, дискретный случай и кросс-корреляция.

Производная первого порядка

Непрерывный случай. Производная функции в точке :

$f'(x) =\frac{df}{dx} = \lim_{h \to 0} \frac{f(x+h) - f(x)}{h}$

Дискретный случай. Рассмотрим дискретную ленту, на которой клетки располагаются с шагом . Тогда координаты ячеек можно обозначить , где , а значение функции в обозначим . Тогда производную в точке можно аппроксимировать используя центральную разность:

$D \, f_i \approx \frac{f_{i+1}-f_{i-1}}{2h}$

Порядок погрешности при этом будет . Для ограниченной ленты необходимо определить граничные условия. Если лента состоит из ячеек с индексами $i = 0, 1, 2, \ldots, N-1$ , то можно замкнуть ленту задав значения

$f_{-1} = f_{N-1}, \quad f_{N} = f_0$

Кросс-корреляция. Вычислять производные на ленте удобно при помощи кросс-корреляции

$D_h f = f \star K^{(1)}$

с ядром

$K^{(1)} = \left[ -\frac{1}{2h}, 0, \frac{1}{2h} \right]$

оператор $\star$ обозначает кросс-корреляцию.

Кросс-корреляция 1D

Результатом кросс-корреляции массива с ядром в позиции будет массив $Y = X \star K$ , элементы которого определяются скалярным произведением фрагмента $(x_{i-1}, x_i, x_{i+1})$ массива с массивом :

$y_i \; = \; k_0 x_{i-1} + k_1 x_{i} + k_2 x_{i+1}$

Значения на краях массиваопределяются в соответствии с заданными граничными условиями.

Производная второго порядка

Непрерывный случай. Вторая производная функции :

$f''(x) = \frac{d}{dx} \left( \frac{df}{dx} \right) = \frac{d^2f}{dx^2}$

Дискретный случай на ленте по аналогии с разностным оператором первого порядка:

$D^{(2)}f_i \approx \frac{f_{i-1}-2f_i+f_{i+1}}{h^2}$

Кросс-корреляция. Вычислять вторые производные на ленте удобно при помощи кросс-корреляции

$D^{(2)}_h f = f \star K^{(2)}$

с ядром

$K^{(2)} = \left[ \frac{1}{h^2}, -\frac{2}{h^2} , \frac{1}{h^2} \right]$

Эта аппроксимация также имеет погрешность порядка .

Градиент

Пусть имеется скалярное поле. Градиент в точке – это вектор, указывающий направление наискорейшего возрастания функциив заданной точке.

Непрерывный случай. Будем рассматривать градиент в двумерном пространстве:

$\mathrm{grad} \, f = ∇ ·f(x, y) = \left( \frac{∂f}{∂x} , \frac{∂f}{∂y} \right)$

где $∇=\mathbf{i}\frac{∂}{∂x} +\mathbf{j} \frac{∂ }{∂y}$ – двумерный оператор набла (оператор Гамильтона), $\mathbf{i}, \mathbf{j}$ – ортогональные единичные векторы.

Дискретный случай. Рассмотрим прямоугольную сетку с шагом . Координаты клеток можно обозначить и , где , а значение функции в клетке обозначим $f_{i,j}$ . Тогда, центральные разности можно записать так:

$D_x \, f_{i,j} \approx \frac{f_{i+1,j}-f_{i-1,j}}{2h}, \quad D_y \, f_{i,j} \approx \frac{f_{i,j+1}-f_{i,j-1}}{2h}$

Дискретный градиент в клетке будет представлять собой вектор:

$∇_h f_{i,j} = \left( D_x \, f_{i,j}, \; D_y \, f_{i,j} \right)$

Если сетка ограничена индексами $i = 0, 1, 2, \ldots, N-1$ и $j = 0, 1, 2, \ldots, M-1$ , то необходимо задать граничные условия. По аналогии с замкнутой лентой, можно замкнуть сетку в тор:

$f_{-1,j} = f_{N-1,j}, \quad f_{N,j} = f_{0,j}$ $f_{i,-1} = f_{i,M-1}, \quad f_{i,M} = f_{i,0}$

Кросс-корреляция. Вычислять градиент на сетке удобно при помощи кросс-корреляции с ядрами. Для производной по можно использовать ядро

$K_x^{(1)} = \frac{1}{2h} \begin{bmatrix} 0 & 0 & 0 \\ -1 & 0 & 1 \\ 0 & 0 & 0 \end{bmatrix}$

а для производной по

$K_y^{(1)} = \frac{1}{2h}\begin{bmatrix} 0 & -1 & 0 \\ 0 & 0 & 0 \\ 0 & 1 & 0 \end{bmatrix}$

Тогда компоненты дискретного градиента для всей сетки можно вычислять так:

$D_x \, f = f \star K_x^{(1)}, \quad D_y \, f = f \star K_y^{(1)}$

Кросс-корреляция 2D

Результатом кросс-корреляции матрицы с ядром

$K = \begin{bmatrix} k_{0,0} & k_{0,1} & k_{0,2} \\ k_{1,0} & k_{1,1} & k_{1,2} \\ k_{2,0} & k_{2,1} & k_{2,2} \end{bmatrix}.$

будет матрица $Y = X \star K$ , элементы $y_{i,j}$ которого определяются как скалярное произведение фрагмента $3 \times 3$ массива (окрестности ) с ядром :

$y_{i,j} \; = \;\sum_{p=0}^{2} \sum_{q=0}^{2} k_{p, \,q} \, · x_{i - 1 + p, \, j - 1 + q}$

Значения на краях матрицы определяются в соответствии с заданными граничными условиями.

Так же для оценки частных производных можно воспользоваться ядрами Собеля. Такой подход обеспечивает дополнительное сглаживание и часто используется при обработке изображений.

Дивергенция

Пусть имеется векторное поле $\mathbf{F}(x,y)$ . Дивергенция в точке – это скалярная величина, показывающая, является ли точка источником или стоком потока.

Непрерывный случай. Пусть задано векторное поле в двумерном евклидовом пространстве

$\mathbf{F}(x,y) = (u(x,y), \, v(x,y)) = \mathbf{i}u + \mathbf{j}v$

Дивергенция этого поля

$\mathrm{div} \,\mathbf{F}(x,y) = \nabla \cdot \mathbf{F}(x,y) =\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y}.$

Дискретный случай. Для прямоугольной сетки, частные производные можно аппроксимировать центральными разностями, аналогично дискретному градиенту:

$D_x \, u_{i,j} \, \approx \, \frac{u_{i+1,j} - u_{i-1,j}}{2h},\quad D_y \, v_{i,j} \, \approx \, \frac{v_{i,j+1} - v_{i,j-1}}{2h}.$

Дискретная дивергенция в клетке тогда задаётся как

$\nabla_h \cdot \mathbf{F}_{i,j} \, \approx \, D_x \, u_{i,j} + D_y \, v_{i,j}$

Порядок погрешности для такой схемы – , как и у центральных разностей для первой производной.

Если сетка ограничена индексами $i = 0, 1, 2, \ldots, N-1$ и $j = 0, 1, 2, \ldots, M-1$ , то необходимо задать граничные условия. Как и для скалярной функции, можно замкнуть её в тор, задав периодические условия для обеих компонент вектора:

$u_{-1,j} = u_{N-1,j}, \quad u_{N,j} = u_{0,j}, \qquad u_{i,-1} = u_{i,M-1}, \quad u_{i,M} = u_{i,0},$ $v_{-1,j} = v_{N-1,j}, \quad v_{N,j} = v_{0,j}, \qquad v_{i,-1} = v_{i,M-1}, \quad v_{i,M} = v_{i,0}.$

Кросс-корреляция. Аналогично градиенту, удобно выразить дискретную дивергенцию через кросс-корреляцию с маленькими ядрами. Можно использовать те же ядра, что и для градиента. Ядро $K_x^{(1)}$ применяется к компоненте , а ядро $K_y^{(1)}$ применяется к . Тогда дискретная дивергенция получается как сумма результатов:

$\nabla_h \cdot \mathbf{F} = u \star K_x^{(1)} + v \star K_y^{(1)},$

Лапласиан

Пусть имеется скалярное поле . Лапласиан в точке – это скаляр, показывающий, насколько значение функции в этой точке отличается от среднего в малой окрестности, т.е. степень ее локальной вогнутости или выпуклости. Если $\Delta f(x, y) > 0$ , то функция выпукла вниз в точке , а при $\Delta f(x, y) < 0$ – выпукла вверх.

Непрерывный случай. В двумерном пространстве лапласиан скалярной функции определяется как дивергенция градиента:

$\Delta f(x, y) = \nabla \cdot (\nabla f) = \frac{\partial^2 f}{\partial x^2} + \frac{\partial^2 f}{\partial y^2}.$

Дискретный случай. Для прямоугольной сетки, частные производные можно аппроксимировать центральными разностями второго порядка:

$D^{(2)}_x \, f_{i,j}\approx \frac{f_{i+1,j} - 2f_{i,j} + f_{i-1,j}}{h^2},$ $D^{(2)}_y \, f_{i,j}\approx \frac{f_{i,j+1} - 2f_{i,j} + f_{i,j-1}}{h^2}.$

Тогда дискретный лапласиан в точке можно будет представить в виде суммы:

$\Delta_h \, f_{i,j} = D^{(2)}_x \, f_{i,j} + D^{(2)}_y \, f_{i,j} \approx \frac{f_{i+1,j} + f_{i-1,j} + f_{i,j+1} + f_{i,j-1} - 4 f_{i,j}}{h^2}.$

В качестве граничных условий можно задать, к примеру, условия Дирихле, фон Неймана или замкнуть плоскость в тор.

Кросс-корреляция. Вычислять лапласиан на сетке удобно при помощи кросс-корреляции с ядром

$K_{\Delta} = \frac{1}{h^2} \begin{bmatrix} 0 & 1 & 0 \\ 1 & -4 & 1 \\ 0 & 1 & 0 \end{bmatrix}$

Тогда значение дискретного лапласиана $\Delta_h \, f_{i,j}$ получается как результат кросс-корреляции матрицы $f_{i,j}$ с ядром $K_\Delta$ в точке с учётом выбранных граничных условий:

$\Delta_h f = f \star K_{\Delta}$

Симуляция диффузии на клеточном автомате

Классическое уравнение диффузии:

$\frac{\partial u(\mathbf{r}, t)}{\partial t} = D \,\Delta u(\mathbf{r}, t)$

где $u(\mathbf{r}, t)$ – концентрация в точке $\mathbf{r}$ в момент , – коэффициент диффузии.

В двумерном случае $\mathbf{r} = (x, y)$ , поэтому

$\Delta u =\frac{\partial^2 u}{\partial x^2} +\frac{\partial^2 u}{\partial y^2}$

Пусть имеется прямоугольная сетка с шагом . Лапласиан функции может быть вычислен через кросс-корреляцию:

$\Delta_h u = u * K_{\Delta}$

где ядро

$K_{\Delta} =\begin{bmatrix} 0 & 1 & 0 \\ 1 & -4 & 1 \\ 0 & 1 & 0 \end{bmatrix}$

При этом, для каждой клетки на временно́м шаге получим значение

$\Delta_h u^{(n)}_{i,j} = u_{i+1,j}^n + u_{i-1,j}^n + u_{i,j+1}^n + u_{i,j-1}^n - 4u_{i,j}^n$

Производную по времени заменяем разностью

$\frac{∂ u}{∂ t} = \frac{u^{(n)} - u^{(n-1)}}{\delta t}$

При $\delta t = 1$ для каждой клетки

$\frac{∂ u}{∂ t} = u_{i,j}^{(n)} - u_{i,j}^{(n-1)}$

Получаем итерационную формулу

$u_{i,j}^{(n)} \; ← \; u_{i,j}^{(n-1)} + D(u_{i+1,j}^n + u_{i-1,j}^n + u_{i,j+1}^n + u_{i,j-1}^n - 4u_{i,j}^n )$

или через кросс-корреляцию

$u^{(n)} \; ← \; u^{(n-1)} + D· u^{(n)} \star K_{\Delta}$

Код реализации на Python диффузии двух веществ на клеточном автомате есть на GitHub. Здесь я приведу содержательный фрагмент – определение функции, выполняющей шаг диффузии. Концентрация отдельного вещества в соответствующей клетке хранится в двумерном массиве u. На каждом шаге вычисляется дискретный лапласиан массива путем вычисления двумерной кросс-корреляции. Для этого используется функция correlate2d из модуля signal библиотеки scipy.

import numpyas np
from scipy.signal import correlate2d

def diffusion_step_box(u, D, kernel=None):
    if kernel is None:
        laplacian_kernel = np.array([
            [0.0,  1.0, 0.0],
            [1.0, -4.0, 1.0],
            [0.0,  1.0, 0.0]
        ], dtype=float)

    laplacian = correlate2d(u, laplacian_kernel,
                            mode='same',      # размер выхода = размер входа
                            boundary='symm')  # зеркальное продолжение (Neumann)
    return u + D * laplacian

Заключение

Здесь были рассмотрены лишь базовые дифференциальные операторы и простые ядра. Но если уловить суть, то остальное можно оперативно разобрать при необходимости. Для каждого оператора существуют множество вариантов ядер, причем различных размеров. Чаще, это квадратные матрицы с нечетным числом рядов и колонок, но в общем случае, применяются и прямоугольные матрицы.