Równanie Poissona: nałóż pełny gradient jako warunek brzegowy za pomocą mnożników Lagrange'a

Mam problem fizyczny rządzony równaniem Poissona w dwóch wymiarach Mam pomiary dwóch składników gradientu i wzdłuż pewnej części granicy, , więc chciałbym narzucić

- \nabla^{2} u = fa (x, y), ja n Ω

$-\nabla^2 u = f(x,y), \; in \; \Omega$

\partial u / \partial x

$\partial{u}/\partial{x}$

\partial u / \partial y

$\partial{u}/\partial{y}$

Γ_{m}

$\Gamma_m$

rozmnażają się w dalekim polu.

{\frac{\partial u}{\partial x_{ja}}}_{0} = {sol}_{m}, o n Γ_{m}

$\frac{\partial{u}}{\partial{x_i}}_0 = g_m, \; on \; \Gamma_m$

Styczny składnik gradientu , mogę po prostu zintegrować, a następnie wymusić przez warunek Dirichleta, taki że $\frac{\partial{u}}{\partial{x}}_{(t,0)}$ Aby jednocześnie narzucić normalny składnik,

\int_{Γ_{m}} {\frac{\partial u}{\partial x}}_{(t, 0)} re s = u_{0}

$\int_{\Gamma_m}\frac{\partial{u}}{\partial{x}}_{(t,0)} \, ds = u_0$

, zrozumiałem, że będę musiał przejść przez mnożniki Lagrange'a.

{\frac{\partial u}{\partial x}}_{(n, 0)}

$\frac{\partial{u}}{\partial{x}}_{(n,0)}$

Myślę więc, że forma wariacyjna to Spędziłem dużo czasu próbując poskładać je razem z informacjami na temat powiązanych problemów, takich jak https://answers.launchpad.net/fenics/+question/212434 https://answers.launchpad.net/fenics/+question / 216323

\int_{Ω} \nabla u \cdot \nabla v re x - λ \int_{Γ_{m}} ({\frac{\partial u}{\partial x}}_{(n, 0)} - {sol}_{m}) v re s = \int_{Ω} fa v re x

$\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$

ale wciąż nie widzę, gdzie się mylę. Moja dotychczasowa próba rozwiązania to:

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)

która działa, ale daje głośny wynik, który wcale nie przypomina rozwiązania równania Poissona. Wydaje się, że ma to coś wspólnego z połączonymi przestrzeniami funkcji, ale nie mogę znaleźć błędu.
Byłbym wdzięczny za każdą pomoc lub wskazówki we właściwym kierunku - wielkie dzięki już!
Pozdrawiam
Markus

finite-element python fenics

— Markus
źródło

Pozwól mi to wyjaśnić: masz zarówno dane Dirichleta, jak i Neumanna, ale tylko część granicy?

— Christian Clason

Jak zrozumiałem PO, na granicy podano gradient. Dane Dirichleta są używane do narzucenia stycznej pochodnej. Pomyślałem, że to dziwne narzucenie Dirichletowi i Neumannowi jednej części granicy, ale może w tej konkretnej sytuacji jest to spójne. Problemem jest raczej zastosowanie danych gradientu na granicy (przez mnożniki).

— stycznia

To prawda, że dałoby to spójne dane, ale nadal masz problem z brakiem stabilności i faktem, że masz warunki brzegowe tylko na części granicy.

— Christian Clason

Ok, dam więcej informacji na temat konkretnego problemu fizycznego, który próbuję rozwiązać. Mam statyczne pole magnetyczne, które mogę rozsądnie założyć, że jest obrotowo symetryczne, a więc 2D. Mierzę składowe promieniowe i osiowe wektora gęstości pola magnetycznego wzdłuż linii, dość blisko osi obrotu i chciałbym widzieć to pole magnetyczne w znacznej odległości od tej osi obrotu. Połączenie Dirichleta i Neumanna BC było moim pomysłem na podejście do problemu, jak wymownie opisał Jan - narzucone dane gradientu na granicy.

— Markus

OK, to znacznie zmienia rzeczy. Masz więc nieograniczoną domenę i informacje pochodne na całej „skończonej” części granicy?

— Christian Clason

Odpowiedzi:

Po pierwsze, ogólna uwaga: nie można określić arbitralnych warunków brzegowych dla operatora różniczkowania cząstkowego i oczekiwać, że równanie różniczkowe cząstkowe (które zawsze obejmuje zarówno warunki operatora, jak i warunki brzegowe) jest dobrze postawione, tzn. Dopuszcza unikalne rozwiązanie, które zależy od dane - z których wszystkie są niezbędnym warunkiem rzeczywistej próby obliczenia czegoś.

W zależności od operatora często można nałożyć wiele ważnych warunków (aby posmakować, zobacz trzytomową monografię Lions i Magenes). Jednak to, co próbujesz zrobić (określić pełny gradient, który jest równoważny zarówno warunkom Dirichleta, jak i Neumanna na tej samej (części) granicy eliptycznej PDE drugiego rzędu), nie jest wśród nich - jest to znane jako boczny problem Cauchy'egoi jest źle postawiony: nie ma gwarancji, że dana para danych granicznych dopuszcza rozwiązanie, a nawet jeśli takie istnieje, nie ma stabilności w odniesieniu do niewielkich zaburzeń w danych. (W rzeczywistości jest to pierwotnie źle postawiony problem w sensie Hadamarda i klasyczny przykład, dlaczego nie można ignorować warunków brzegowych przy omawianiu dobrej pozycji. Można znaleźć wyraźny przykład w jego Wykładach na temat problemu Cauchy'ego w częściowej różnicy liniowej równania z lat dwudziestych.)

Przejdźmy teraz do konkretnego problemu (który może być podręcznikowym przykładem problemu XY ). Jeśli dobrze rozumiem, chcesz rozwiązać problem zewnętrzny równania Poissona i weź dwuwymiarową domenę konformacyjną . Masz pełne dane (Dirichlet - używając sztuczki stycznej pochodnej - i Neumann) na (powiedzmy) . Jeśli jest wystarczająco duży, możesz uzasadnić stan promieniowania (który określa asymptotyczne zachowanie twojego rozwiązania przy ); patrz np. książka $(r,R)\times(a,b)$ $x=r$ $R$ $x\to \infty$ Metody numeryczne problemów zewnętrznych autorstwa Long'an Ying. Pytanie brzmi: co wiesz o pozostałych częściach granicznych, i . $y=a$ $y=b$

Jeśli potrafisz narzucić warunki brzegowe (Neumann, Robin, Dirichlet - które, nawiasem mówiąc, musielibyście ustalić stałą w całkowaniu pochodnej stycznej), wystarczy użyć albo normalnych składników gradientu jako warunku Neumanna (jeśli można naprawić tryb stały) lub zintegrować komponenty styczne jako warunek Dirichleta. Ponieważ oba warunki prawdopodobnie odpowiadają tej samej funkcji, w obu przypadkach otrzymujesz to samo rozwiązanie.
Jeśli nie znasz zachowania w i , naprawdę masz boczny problem Cauchy'ego i nie możesz obliczyć rozwiązania przy użyciu standardowych metod elementów skończonych. Standardowym sposobem radzenia sobie z tym jest metoda quasi-odwracalności (wprowadzona przez Lattésa i Lionsa w latach 60. XX wieku; patrz ich książka ), która polega na przybliżeniu problemu drugiego rzędu przez problem czwartego rzędu (gdzie można - i trzeba - przepisać dwa warunki brzegowe). W twoim przypadku oznaczałoby to zastąpienie przez $y=a$ $y=b$ $-\Delta u=f$ dla niektórych małych . (Można to również zinterpretować jako zminimalizowanie pozostałości w odpowiedniej normie i dodanieterminu regularyzacji ; związanego z twoim komentarzem na temat traktowania problemu jako problemu odwrotnego.) Możesz to pokazać dla kompatybilnych danych (tj. Pary granic warunki, które faktycznie odpowiadają rozwiązaniu równania Poissona), quasi-odwracalne rozwiązania zbieżne do jako . $-\Delta u + \varepsilon \Delta^2 u=f$ $\varepsilon>0$ $H^2$ $u$ $u_\varepsilon$ $u$ $\varepsilon\to0$

$H^2$

— Christian Clason
źródło

Wdrożenie przez mieszane elementy w FEniCS patrz biharmonicdemo. Prawdopodobnie jest to bez terminu Laplace'a, ale myślę, że można go łatwo dodać.

— Jan Blechta

Cześć Christian, dziękuję za twoją sugestię! Miałem wrażenie, że równanie Poissona jest łagodne, jeśli chodzi o stabilność numeryczną - dziękuję za zwrócenie na to uwagi. Przeczytam o tym, jak zasugerowałeś. Nie jestem pewien, czy to znacznie zmienia sytuację, ale jak wspomniano w komentarzu powyżej, kombinacja Dirichleta-Neumanna może być myląca. „Wszystko”, którego szukam, to sposób nałożenia danych gradientu na granicy.

— Markus

Równanie Poissona jest łagodne, ale nie jest to równanie, które próbujesz rozwiązać :) (Warunki brzegowe są integralną częścią równania; sam operator nie jest wystarczający, aby zdecydować o stabilności.)

— Christian Clason

Dobra, to daje mi coś do przeżuwania. Dziękuję wszystkim za poświęcony czas, radę i cierpliwość - i przepraszam za wpadnięcie w pułapkę XY ...

— Markus

Nie możesz oczekiwać, że rozwiązanie twojego zmienionego problemu będzie rozwiązaniem problemu Poissona, ponieważ musisz jakoś zmienić problem, aby był dobrze postawiony.

fa (u, λ) = \int_{Ω} \frac{1}{2)} | \nabla u |^{2)} re x - \int_{Ω} fa u re x - \int_{Γ_{N.}} sol u re S. + \int_{Γ_{N.}} λ (u - u_{re}) re S.

$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, λ) \in V \times L^{2} (Γ_{N})

$(u,\lambda) \in V\times L^2(\Gamma_\mathrm{N})$

V = {v \in H^{1}; v |_{Γ_{D}} = 0}

$V = \{v\in H^1;v|_{\Gamma_\mathrm{D}}=0\}$

Γ_{D}

$\Gamma_\mathrm{D}$

u_{D}

$u_\mathrm{D}$

Γ_{N}

$\Gamma_\mathrm{N}$

F (u)

$F(u)$

0 = re fa (u) [v] = \int_{Ω} \nabla u \cdot \nabla v re x - \int_{Ω} fa v re x - \int_{Γ_{N.}} sol v re S. \forall v \in V.,

$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,$

Γ_{N}

${\Gamma_\mathrm{N}}$

Γ_{D}

${\Gamma_\mathrm{D}}$

0 = re fa (u, λ) [v, μ] = \int_{Ω} \nabla u \cdot \nabla v re x - \int_{Ω} fa v re x - \int_{Γ_{N.}} sol v re S. + \int_{Γ_{N.}} λ v re S. + \int_{Γ_{N.}} μ (u - u_{re}) re S. \forall (v, μ) \in V. \times {L.}^{2)} (Γ_{N.}),

$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}),$

- Δ u = f

$-\Delta u=f$

\frac{\partial u}{\partial n} = g - λ

$\frac{\partial u}{\partial\mathbf{n}} = g-\lambda$

Γ_{N}

$\Gamma_\mathrm{N}$

Γ_{N}

$\Gamma_\mathrm{N}$

$\lambda\ll|g|$

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

Wniosek jest taki, że nie można oczekiwać, że PDE drugiego rzędu zaakceptuje dwa niezależne warunki brzegowe.

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

— Jan Blechta
źródło

\nabla^{2} u - f

$\nabla ^2 u - f$

u \in H^{1}

$u\in H^1$

Moja matematyka niestety nie jest do zera i nie jestem pewien co do matematycznych implikacji przestrzeni Banacha, ale staram się zrozumieć, dlaczego to równanie nie jest rozwiązaniem równania Poissona, gdy zanika termin mnożnika Lagrange'a. Z fizycznego punktu widzenia rozwiązanie (opisanego przeze mnie problemu praktycznego, nie mam na myśli rozwiązania matematycznego) musi istnieć, o ile widzę

— Markus

Jest to więc raczej odwrotny problem, znalezienie warunku granicznego dla dalekiego pola, które wraz z warunkiem Dirichleta, który można narzucić, daje obserwowany gradient normalny na granicy, na której mierzysz?

— Markus

Twoje podejście nie może działać, zdecydowanie z powodu wdrożenia i prawdopodobnie z powodu sformułowania.

Narzucając warunki Dirichleta w dolfinie, ostatecznie ustawia odpowiednie DOF twojego obszaru testowego na zero.

Oto fragment podręcznika fenics :

Rozdział 3.3.9 (koniec): Zastosowanie warunku brzegowego Dirichleta do układu liniowego określi wszystkie stopnie swobody, które należy ustawić na podaną wartość, i zmodyfikuje układ liniowy tak, aby jego rozwiązanie spełniało warunki brzegowe. Odbywa się to przez wyzerowanie i wstawienie 1 na przekątnej rzędów macierzy odpowiadających wartościom Dirichleta i wstawienie wartości Dirichleta w odpowiednim wpisie wektora po prawej stronie [...]

$v$ $\Gamma_m$ , co powoduje, że twój termin z mnożnikiem również zniknie.

Podsumowując, używając domyślnej procedury w dolfin, nie można zastosować zarówno Dirichleta, jak i innych warunków na tej samej granicy.

Zanim jednak spróbujesz naprawić to w swojej implementacji, znajdź odpowiednie przestrzenie testowe dla swojego sformułowania matematycznego (jak właśnie wspomniał @Jan Blechta).

— Jan
źródło

Rozumiem twój punkt widzenia - myślę, że moje sformułowanie może nie odzwierciedlać dokładnie tego, co wdrożyłem - moje przeprosiny. Zasada wariacyjna jest tylko mglistym wspomnieniem i staram się znów to zrozumieć. Przeczytałem instrukcję do przodu i do tyłu wraz z kilkoma przykładami kodu FEniCS implementującego mnożniki Lagrange'a. Myślałem, że podniesiony przez ciebie problem jest właśnie powodem, dla którego użyjesz drugiej funkcji testowej „d”.

— Markus

Zgadzam się z @JanBlechta. Najpierw musisz znaleźć odpowiednią przestrzeń dla mnożnika, co nie jest łatwe. Być może teksty na temat optymalizacji ograniczeń PDE, w których używa się mnożników do uwzględnienia warunków ubocznych, dostarczą pomocnych pomysłów. W tym artykule mnożnik służy do uwzględnienia zależnych od czasu warunków Dirichleta.

— stycznia