Как реализовать дельта-функцию при сеточном решении УЧП

[misternickel]

Здравствуйте. Вопрос мой глобально звучит так: как численно реализовать дельта-функцию Дирака при решении уравнения в частных производных методом конечных разностей?

Суть вопроса: как именно реализовать дельта-функцию как источник в численном решении? Нужно ли сосредоточить источник с интегралом единица в одной конкретной ячейке сетки пространства (с координатами х0 и z0), или размазать источник по нескольким соседним ячейкам? Если размазать, то какую функцию выбрать, и насколько сильно размазывать?

Суть задачи: решается следующее (сопряжённое, но это не важно) уравнение:

$\frac{\partial C}{\partial t}+U\cdot\frac{\partial C}{\partial x} = \frac{\partial}{\partial z}K\frac{\partial C}{\partial z} + \delta (x-x_0) \delta (z-z_0)$

где x и z - независимые переменные пространства, t - время, C - зависимая переменная (концентрация вещ-ва), U - скорость ветра, K - коэффициент диффузии, $\delta ()$ - дельта-функция Дирака, x_0 , z_0 - известные константы. Есть какие-то начальные и граничные условия, я их дам, если нужно.

Уравнение решается с помощью обычного конечно-разностного алгоритма - делается сетка по независимым переменным, происходит дискретизация уравнения. Само решение осуществляется методом расщепления по направлениям, схема Кранка-Николсон по z, схема предиктор-корректор по х, если это важно.

Вроде бы, такая постановка задачи не нова, но нигде в известных учебниках по численным методам не написано, как аппроксимировать дельта-функцию на практике в сеточных задачах. Буду рад и подсказкам, и отсылкам к книгам-статьям.

[arseniiv]

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

Если подумать наивно, то сосредоточенная в одной точке сетки функция намного больше похожа на дельта-функцию, чем сосредоточенная в нескольких. Если подумать чуть менее наивно и взять, например, обратное преобразование Фурье от константы на сетке или поискать функцию, уместную в теореме $\delta*f = f$ ($*$ — свёртка, здесь дискретная) или, наконец, «дискретизировав» само её определение — получится опять функция, сосредоточенная в одной точке. По крайней мере, если сетка прямоугольная, и если ничего не упускаю.

[amon]

misternickel в сообщении #1258977 писал(а):

как именно реализовать дельта-функцию как источник в численном решении?

IMHO, ни один из Ваших способов не годится. Надо либо преобразовать уравнение, как советовал arseniiv. Если задача решается на всей плоскости - преобразование Фурье, тогда вообще численно надо только интегралы считать. Если задача решается в хитрой области с затейливыми гран.условиями, то дело - дрянь, и придётся $\delta$-функцию прятать в граничное условие в узле с координатами $x_0,z_0$ (разрыв по $x$ и разрыв производной по $z$). Что там получится с устойчивостью вычислительной схемы - бог весть.

[misternickel]

amon в сообщении #1259126 писал(а):

misternickel в сообщении #1258977 писал(а):

как именно реализовать дельта-функцию как источник в численном решении?

IMHO, ни один из Ваших способов не годится. Надо либо преобразовать уравнение, как советовал arseniiv. Если задача решается на всей плоскости - преобразование Фурье, тогда вообще численно надо только интегралы считать. Если задача решается в хитрой области с затейливыми гран.условиями, то дело - дрянь, и придётся $\delta$-функцию прятать в граничное условие в узле с координатами $x_0,z_0$ (разрыв по $x$ и разрыв производной по $z$). Что там получится с устойчивостью вычислительной схемы - бог весть.

Задача решается на всей плоскости, сетка прямоугольная с узлами фиксированного размера, граничные условия - элементарные (поскольку эта задача сопряжённая): по всем границам $C=0$ (условие первого рода), на нижней границе поток равен 0 (условие второго рода).

На ресерчгейте посоветовали со ссылкой на опыт довольно простой и естественный приём - аппроксимировать полиномом 7 степени с 6 условиями (ноль функции на краевых точках, ноль производной функции на краевых точках, максимум в серединном узле и нулевая производная там), но на сколько ячеек размазывать функцию, не сказали. Там же дали ссылку вот на это (DOI: 10.1017/S0962492902000077), где предлагается для аппроксимации дельта функции в очень сложной задаче вот такая функция, которая всем требованиям из предложения выше соответствует (также в узловых точках функция и её производная непрерывны):

$\delta_h(x)=\frac{1}{h}\cdot f(\frac{x}{h}) $

$ f(r)=$$\left\{\begin{array}{rcl} \frac{1}{4}\cdot(1+\cos(\frac{\pi\cdot r}{2})) &if& \left\lvert r \right\rvert \leqslant2\\ 0 &if& \left\lvert r \right\rvert >2\\ \end{array} \right.$$ $

Здесь $ h $ - шаг сетки, $ x $ - декартова координата. Это пятиточечный в одномерном случае шаблон, охватывающий 4 соседних ячейки. В сеточном виде это выглядит как треугольник, сеточный интеграл по которому равен единице, как и должно быть.

[misternickel]

arseniiv в сообщении #1259020 писал(а):

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

Если подумать наивно, то сосредоточенная в одной точке сетки функция намного больше похожа на дельта-функцию, чем сосредоточенная в нескольких. Если подумать чуть менее наивно и взять, например, обратное преобразование Фурье от константы на сетке или поискать функцию, уместную в теореме $\delta*f = f$ ($*$ — свёртка, здесь дискретная) или, наконец, «дискретизировав» само её определение — получится опять функция, сосредоточенная в одной точке. По крайней мере, если сетка прямоугольная, и если ничего не упускаю.

Я тоже подождал некоторое время. Математиком по образованию не являюсь, поэтому что-то утверждать в разных её областях с уверенностью не могу, но пытаюсь рассуждать наиболее естественным образом. Вот в сообщении выше настоящего я привёл функцию, которая естественным, на мой взгляд, требованиям удовлетворяет: имеет сеточный и аналитический интеграл единица (что должно быть исходя из определения дельта-функции), содержит всего три ненулевых узла (не сильно тянет по сетке функцию), даёт непрерывную производную, нулевую на краях, симметрична (что наверняка хорошо для устойчивости схемы).

Просто не понимаю, зачем здесь может быть нужно что-то ещё, на основании каких соображений. Какую бы функцию не удалось получить предложенным Вами способом, для неё остаётся актуальным вопрос про ширину размытия. И непонятно, какими свойствами она будет превосходить предложенную выше функцию. В этой связи интересен вопрос о полном наборе критериев, которым должна удовлетворять сеточная аппроксимация дельта-функции. Как-то он нигде не рассматривается.

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

В конце концов, вряд ли люди в процитированной публикации - полные идиоты. Имеет смысл им довериться.

[amon]

misternickel в сообщении #1259424 писал(а):

И непонятно, какими свойствами она будет превосходить предложенную выше функцию.

Беда в том, что $\delta$-функция это не функция вовсе. Это функционал, определяемый как $\int \delta(x)f(x)dx=f(0)$ для любой $f(x)$ из некоторого достаточно широкого класса функций. И что-то сильно сомневаюсь я, что что-нибудь кроме граничного условия в узле, получающегося из равенства $\delta(x)=\frac{d\theta(x)}{dx}$ ($\theta(x)$ это единичная ступенька) её как-нибудь на сетке аппроксимирует. Я с таким раз в жизни сталкивался, но в моем случае повезло - интегральное уравнение, получившееся после преобразования Фурье исходного уравнения, решилось аналитически, и все свелось к взятию на компьютере обратного преобразования.

[misternickel]

amon в сообщении #1259429 писал(а):

misternickel в сообщении #1259424 писал(а):

И непонятно, какими свойствами она будет превосходить предложенную выше функцию.

Беда в том, что $\delta$-функция это не функция вовсе. Это функционал, определяемый как $\int \delta(x)f(x)dx=f(0)$ для любой $f(x)$ из некоторого достаточно широкого класса функций. И что-то сильно сомневаюсь я, что что-нибудь кроме граничного условия в узле, получающегося из равенства $\delta(x)=\frac{d\theta(x)}{dx}$ ($\theta(x)$ это единичная ступенька) её как-нибудь на сетке аппроксимирует. Я с таким раз в жизни сталкивался, но в моем случае повезло - интегральное уравнение, получившееся после преобразования Фурье исходного уравнения, решилось аналитически, и все свелось к взятию на компьютере обратного преобразования.

Чувствую примерно, о чём Вы хотите сказать. Спасибо, буду держать в уме наличие такой опции, хоть и вряд ли в моём случае это возможно - вид функций скорости ветра и коэффициента диффузии очень сложный.

А можно Вас попросить подтвердить правильность одного утверждения? Вот оно: пусть есть непрерывная функция $C=C(x,z,t), x\in [0;X], z\in [0;Z], t \in [0;T],$, все обозначения совпадают с указанными в стартовом посте. Тогда

$$\int\limits_{0}^{T}\int\limits_{0}^{X}\int\limits_{0}^{Z} \delta(x-x_0) \delta(z-z_0) C(x,z,t)dxdzdt = \int\limits_{0}^{T} C(x_0,z_0,t)dt$$

[amon]

Если точки $x_0$ и $z_0$ лежат внутри области интегрирования, то все правильно.

[misternickel]

amon в сообщении #1259439 писал(а):

Если точки $x_0$ и $z_0$ лежат внутри области интегрирования, то все правильно.

Да, внутри.

Большое спасибо! Было бы глупо заниматься всем этим и не проверить ключевое преобразование, за счёт которого я связываю модель и реальные измерения. Дельта-функция хоть и плоха своей непонятностью, но чудо как хороша для введения наших реальных данных (одноточечных в пространстве) в модель, особенно в случае сопряжённых уравнений.

[misternickel]

misternickel в сообщении #1258977 писал(а):

Вроде бы, такая постановка задачи не нова, но нигде в известных учебниках по численным методам не написано, как аппроксимировать дельта-функцию на практике в сеточных задачах.

Подведу итог своим поискам, вдруг кто-то заинтересуется или столкнётся с такой задачей. Никакой особенной специфики включения дельта-функции в разностные схемы нет. Если есть конвективные члены, то используйте только монотонные схемы, поскольку решения быстропеременные из-за влияния точечного источника. Я использовал явно-неявную схему.

Рабочие варианты аппроксимации дельта-функции даны во многих научных статьях. Для широкого интервала точек см., например, статью Waldén J. On the approximation of singular source terms in differential equations // Numerical Methods for Partial Differential Equations. – 1999. – Т. 15. – №. 4. – С. 503-520. Я взял пятиточечную именно оттуда.

Как я понял, что всё верно? Очень просто. Решил прямую задачу при некоем значении параметра, затем взял это решение и по нему с помощью решения сопряжённого уравнения восстановил исходное значение параметра. Значения почти совпали: 0.69 - исходное значение параметра, 0.74$\pm$0.11 получилось для серии из 29 численных экспериментов, в которых сильно варьировали другие параметры. От дробности сетки результат не зависел.

Ещё мне на почту написал очень крутой спец из ВШЭ Владимир Гордин, увидев мой вопрос на ресрчгейте. Он дал вот такой совет:

Цитата:

Разностные схемы хорошо работают на гладких решениях, к каковым дельта-функция не относится.
Тут нужна комбинация аналитики и вычислений.
Сначала найдем аналитическое решение $u_1$ с коэффициентом К, зафиксированным в точке особенности.
Это дело нехитрое - преобразование Фурье и Лапласа - дело сделано.
Теперь разность $u-u_1$ удовлетворяет уравнению с гладкой правой частью. Его и вычисляйте разностями.

Я его не использовал, мне это показалось слишком сложным для моего слабого ума. Дай бог, поможет кому-то.