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

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

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

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

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

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

Дискретный случай. Рассмотрим дискретную ленту, на которой клетки располагаются с шагом h. Тогда координаты ячеек можно ​обозначить x_i=x_0 + ih, где i∈ℤ, а значение функции в x_i обозначим f_i. Тогда производную в точке x_i можно аппроксимировать используя центральную разность:

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

Порядок погрешности при этом будет O(h^2). Для ограниченной ленты необходимо определить граничные условия. Если лента состоит из N ячеек с индексами 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

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

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

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

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

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

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]

Эта аппроксимация также имеет погрешность порядка O(h^2).

Градиент

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

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

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}– ортогональные единичные векторы.

Дискретный случай. Рассмотрим прямоугольную сетку с шагом h. Координаты клеток можно ​обозначить x_i=x_0 + ih и y_i=y_0 + jh, где i, j∈ℤ, а значение функции в клетке (x_i, y_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}

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

∇_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}

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

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

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

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

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

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 массива X (окрестности (i,j)) с ядром K:

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

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

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

Дивергенция

Пусть имеется векторное поле mathbf{F}(x,y). Дивергенция в точке (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}.

Дискретная дивергенция в клетке (i,j) тогда задаётся как

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

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

Если сетка ограничена индексами 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)} применяется к компоненте u, а ядро K_y^{(1)} применяется к v. Тогда дискретная дивергенция получается как сумма результатов:

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

Лапласиан

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

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

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}.

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

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 в точке (i,j) с учётом выбранных граничных условий:

Delta_h f=f star K_{Delta}

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

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

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

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

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

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

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

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

Delta_h u=u * K_{Delta}

где ядро

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

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

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 для каждой клетки (i,j)

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

Заключение

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

Telegram

Автор: Ordevoir

Источник

Версия для печати Версия для печати
Rambler's Top100