Уравнение Пуассона: наложение полного градиента как граничного условия с помощью множителей Лагранжа

У меня есть проблемы определяется уравнением Пуассона в двух измерениях

$$-\nabla^2 u = f(x,y), \; in \; \Omega$$ у меня есть измерения двух компонентов градиента$\partial{u}/\partial{x}$ и$\partial{u}/\partial{y}$ вдоль некоторой части границы,$\Gamma_m$ , поэтому хотелосьввести $$\frac{\partial{u}}{\partial{x_i}}_0 = g_m, \; on \; \Gamma_m$$ и размножаются в дальнем поле.

Тангенциальный градиентный компонент, $\frac{\partial{u}}{\partial{x}}_{(t,0)}$, я могу просто интегрировать, а затем применить через условие Дирихле, такое, что Чтобы одновременно наложить нормальную составляющую, ∂

$$\int_{\Gamma_m}\frac{\partial{u}}{\partial{x}}_{(t,0)} \, ds = u_0$$ , я понял, что мне нужно пройти через множители Лагранжа.$\frac{\partial{u}}{\partial{x}}_{(n,0)}$

Так что я думаю , что вариационная форма , то Я потратил много времени, пытаясь собрать воедино информацию о связанных проблемах, такую ​​как https://answers.launchpad.net/fenics/+question/212434https://answers.launchpad.net/fenics/+question / 216323

$$\int_\Omega \nabla{u} \cdot\nabla{v}\, dx - \lambda \int_{\Gamma_m} ( \frac{\partial{u}}{\partial{x}}_{(n,0)}-g_m ) v \,ds = \int_\Omega f\, v\, dx$$

но до сих пор не вижу, где я иду не так. Моя попытка решения до сих пор:

from dolfin import *

# Create mesh and define function space
mesh = UnitSquareMesh(64, 64)
V = FunctionSpace(mesh, "Lagrange", 1)
R = FunctionSpace(mesh, "R", 0)
W = V * R

# Create mesh function over cell facets
boundary_parts = MeshFunction("uint", mesh, mesh.topology().dim()-1)

# Mark left boundary facets as subdomain 0
class LeftBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and x[0] < DOLFIN_EPS

Gamma_Left = LeftBoundary()
Gamma_Left.mark(boundary_parts, 0)

class FarField(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and ( (x[0] > 1.0-DOLFIN_EPS) \
               or (x[1]<DOLFIN_EPS) or (x[1]> 1.0-DOLFIN_EPS) )

Gamma_FF = FarField()
Gamma_FF.mark(boundary_parts, 1)

# Define boundary condition
u0 = Expression("sin(x[1]*pi)")
bcs = [DirichletBC(V, u0, Gamma_Left)]

# Define variational problem
(u, lmbd) = TrialFunctions(W)
(v, d) = TestFunctions(W)

f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)")
g = Constant(0.0)
h = Constant(-4.0)
n = FacetNormal(mesh)

F = inner(grad(u), grad(v))*dx + d*dot(grad(u),n)*ds(0) + lmbd*dot(grad(v),n)*ds(0)-\
   (f*v*dx + g*v*ds(1) + h*d*ds(0) + lmbd*h*ds(0))

a = lhs(F)
L = rhs(F)

# Compute solution
A = assemble(a, exterior_facet_domains=boundary_parts)
b = assemble(L, exterior_facet_domains=boundary_parts)
for bc in bcs: bc.apply(A, b)

w = Function(W)
solve(A, w.vector(), b, 'lu')
(u,lmbd) = w.split()

# Plot solution
plot(u, interactive=True)

который работает, но дает шумный результат, совсем не похожий на решение уравнения Пуассона. Кажется, что-то связано с объединенными функциональными пространствами, но я не могу найти ошибку.
Я был бы признателен за любую помощь или указатели в правильном направлении - большое спасибо уже!
Ура
Маркус

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

В зависимости от оператора, часто существует целый ряд действительных условий, которые вы можете наложить (чтобы почувствовать вкус, см. Трехтомную монографию Lions and Magenes). Однако то, что вы пытаетесь сделать (указать полный градиент, который эквивалентен условиям Дирихле и Неймана на одной и той же (части) границе для эллиптического PDE второго порядка), не входит в их число - это называется боковая проблема Кошии является некорректным: нет никакой гарантии, что данная пара граничных данных допускает решение, и даже если они существуют, нет стабильности в отношении небольших возмущений в данных. (На самом деле, это первоначальная некорректная задача в смысле Адамара и классический пример того, почему вы не можете игнорировать граничные условия при обсуждении корректности. Вы можете найти явный пример в его Лекциях по проблеме Коши в линейном дифференциальном уравнении с частными производными уравнения с 1920-х годов.)

Теперь к вашей конкретной проблеме (которая может быть примером из XY в учебнике ). Если я правильно понимаю, вы хотите решить внешнюю задачу для уравнения Пуассона, и вы берете двумерную расчетную область . У вас есть полные данные (Дирихле - с использованием трюка с тангенциальной производной - и Неймана) по (скажем) x = r . Если R достаточно велико, вы можете обосновать условие излучения (которое предписывает асимптотику вашего решения при x → ∞ ); см., например, книгу$(r,R)\times(a,b)$$x=r$$R$$x\to \infty$Численные методы для внешних задач Лунъань Ин. Вопрос в том, что вы знаете об оставшихся краевых частях: и y = b .$y=a$$y=b$

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

2. Если вы не знаете поведения при и y = b , у вас действительно есть боковая задача Коши, и вы не можете вычислить решение, используя стандартные методы конечных элементов. Стандартным способом решения этой проблемы является метод квазиобратимости (введенный Латте и Лайонсом в 1960-е годы; см. Их книгу ), который заключается в аппроксимации задачи второго порядка задачей четвертого порядка (где вы можете - и надо - прописать два граничных условия). В вашем случае это будет означать замену - Δ u = f на - Δ u + ε Δ $y=a$$y=b$$-\Delta u=f$ для небольшого ε > 0 . (Это также можно интерпретировать как минимизацию остатка в подходящей норме и добавлениечлена регуляризации H 2 ; это связано с вашим комментарием по поводу трактовки проблемы как обратной задачи.) Вы можете показать это для совместимых данных (т. Е. Пары границ условия, которые фактически соответствуют решению u уравнения Пуассона), решения квазиобратимости u ε сходятся к u при ε → 0 .$-\Delta u + \varepsilon \Delta^2 u=f$$\varepsilon>0$$H^2$$u$$u_\varepsilon$$u$$\varepsilon\to0$

   Поскольку теперь это проблема четвертого порядка с решениями в , лучший способ ее численного решения состоит в использовании смешанной конечно-элементной формулировки, такой как описанная в статье Дарде, Ханнукайнена и Хювонена. (В Интернете также есть несколько слайдов .) Не должно быть слишком сложно реализовать этот подход с использованием FEniCS.$H^2$

Вы не можете ожидать, что решение вашей измененной проблемы будет решением проблемы Пуассона, потому что вам нужно как-то изменить проблему, чтобы сделать ее корректной.

$$F(u, \lambda) = \int_\Omega \frac{1}{2}|\nabla u|^2 \;\mathrm{d}x - \int_\Omega f u \;\mathrm{d}x - \int_{\Gamma_\mathrm{N}} g u \;\mathrm{d}S + \int_{\Gamma_\mathrm{N}} \lambda (u-u_\mathrm{D})\;\mathrm{d}S$$ $(u,\lambda) \in V\times L^2(\Gamma_\mathrm{N})$$V = \{v\in H^1;v|_{\Gamma_\mathrm{D}}=0\}$$\Gamma_\mathrm{D}$$u_\mathrm{D}$$\Gamma_\mathrm{N}$$F(u)$ $$0 = \mathrm{D}F(u)[v] = \int_\Omega \nabla u\cdot\nabla v\;\mathrm{d}x - \int_\Omega fv\;\mathrm{d}x - \int_{\Gamma_\mathrm{N}} gv\;\mathrm{d}S \quad \forall v\in V,$$ ${\Gamma_\mathrm{N}}$${\Gamma_\mathrm{D}}$ $$0 = \mathrm{D}F(u,\lambda)[v,\mu] = \int_\Omega \nabla u\cdot\nabla v\;\mathrm{d}x - \int_\Omega fv\;\mathrm{d}x - \int_{\Gamma_\mathrm{N}} gv\;\mathrm{d}S \quad + \int_{\Gamma_\mathrm{N}} \lambda v\;\mathrm{d}S + \int_{\Gamma_\mathrm{N}} \mu (u-u_\mathrm{D})\;\mathrm{d}S \quad \forall (v,\mu)\in V\times L^2(\Gamma_\mathrm{N}),$$ $-\Delta u=f$$\frac{\partial u}{\partial\mathbf{n}} = g-\lambda$$\Gamma_\mathrm{N}$$\Gamma_\mathrm{N}$

$\lambda\ll|g|$

$\Gamma_\mathrm{D}$$v\in V$$\Gamma_\mathrm{D}$

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

---

$F(u,\lambda)$$\mathrm{D}F(u,\lambda)[v,\mu]$derivative()

$F(u,\lambda)$$\lambda\in L^2(\Gamma_\mathrm{N})$$\lambda\in L^2(\Omega)$$\lambda\in\mathbb{R}$

Ваш подход не может работать, определенно из-за реализации и, вероятно, из-за вашей формулировки.

Наложение условий Дирихле на dolfin, в конечном итоге устанавливает соответствующие DOFs вашего тестового пространства на ноль.

Это выдержка из руководства по фенике :

Глава 3.3.9 (конец): Применение граничного условия Дирихле к линейной системе идентифицирует все степени свободы, которые должны быть установлены на заданное значение, и модифицирует линейную систему так, чтобы ее решение учитывало граничное условие. Это достигается путем обнуления и вставки 1 по диагонали строк матрицы, соответствующих значениям Дирихле, и вставкой значения Дирихле в соответствующую запись правого вектора вектора [...]

$v$$\Gamma_m$

Таким образом, используя стандартную процедуру в dolfin, вы не можете применять и Dirichlet, и другие условия на одной границе.

Однако, прежде чем пытаться исправить это в своей реализации, найдите правильные тестовые пространства для вашей математической формулировки (как только что упомянул @Jan Blechta).