Persamaan Poisson: Tetapkan gradien penuh sebagai kondisi batas melalui pengganda Lagrange

11

Saya memiliki masalah fisik yang diatur oleh persamaan Poisson dalam dua dimensi Saya memiliki pengukuran dua komponen gradienu /x danu /y sepanjang beberapa bagian batas, Γ m , jadi ingin memaksakan u

2u=f(x,y),inΩ
u/xu/yΓm dan menyebar ke bidang jauh.
uxi0=gm,onΓm

Komponen gradien tangensial, , saya hanya bisa mengintegrasikan dan kemudian menegakkan melalui kondisi Dirichlet, sehingga Γmuux(t,0) Untuk memaksakan komponen normal secara bersamaan,u

Γmux(t,0)ds=u0
, saya kumpulkan bahwa saya harus menggunakan pengganda Lagrange.ux(n,0)

Jadi saya pikir bentuk variasinya kemudian Saya menghabiskan waktu lama untuk mencoba menyatukannya dari informasi tentang masalah terkait seperti https://answers.launchpad.net/fenics/+question/212434https://answers.launchpad.net/fenics/+question / 216323

ΩuvdxλΓm(ux(n,0)gm)vds=Ωfvdx

tetapi masih tidak bisa melihat di mana saya salah. Upaya solusi saya sejauh ini adalah:

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)

yang berjalan tetapi memberikan hasil yang bising sama sekali tidak menyerupai solusi untuk persamaan Poisson. Tampaknya ada hubungannya dengan ruang fungsi gabungan, tetapi saya tidak dapat menemukan kesalahan.
Saya akan menghargai bantuan atau petunjuk ke arah yang benar - sudah banyak terima kasih!
Cheers
Markus

Markus
sumber
Biarkan saya mengerti: Anda memiliki data Dirichlet dan Neumann, tetapi hanya pada bagian dari batas?
Christian Clason
1
Seperti yang saya pahami OP, itu adalah gradien yang diberikan pada batas. Data Dirichlet digunakan untuk memaksakan turunan tangensial. Saya pikir itu aneh untuk memaksakan Dirichlet dan Neumann pada satu bagian dari batas tetapi mungkin dalam situasi khusus ini konsisten. Jadi, masalahnya bukan bagaimana menerapkan data gradien pada batas (melalui pengganda).
Jan
Benar, itu akan memberikan data yang konsisten, tetapi Anda masih memiliki masalah kurangnya stabilitas, dan fakta bahwa Anda memiliki kondisi batas hanya pada sebagian dari batas.
Christian Clason
Ok, izinkan saya memberi beberapa info lebih lanjut tentang masalah fisik spesifik yang saya coba selesaikan. Saya memiliki medan magnet statis yang bisa saya asumsikan sebagai rotasi simetris, sehingga 2D. Saya mengukur komponen radial dan aksial dari vektor kepadatan medan magnet sepanjang garis, cukup dekat dengan sumbu rotasi dan ingin melihat medan magnet ini pada jarak yang cukup jauh dari poros rotasi ini. Kombinasi Dirichlet dan Neumann BC hanyalah ide saya untuk mendekati masalah seperti yang dijelaskan Jan dengan fasih - menerapkan data gradien pada batas.
Markus
1
OK, itu mengubah banyak hal secara signifikan. Jadi, Anda memiliki domain tak terbatas, dan informasi turunan pada seluruh bagian "terbatas" batas?
Christian Clason

Jawaban:

8

Pertama, poin umum: Anda tidak dapat menentukan kondisi batas arbitrer untuk operator diferensial parsial dan berharap bahwa persamaan diferensial parsial (yang selalu mencakup kondisi operator dan batas) ditempatkan dengan baik, yaitu, mengakui solusi unik yang bergantung terus menerus pada data - yang semuanya merupakan syarat yang diperlukan untuk benar-benar mencoba menghitung sesuatu.

Bergantung pada operator, sering kali ada sejumlah kondisi valid yang dapat Anda terapkan (untuk mendapatkan rasa, lihat monograf tiga volume oleh Lions and Magenes). Namun, apa yang Anda coba lakukan (tentukan gradien penuh, yang setara dengan kondisi Dirichlet dan Neumann pada batas (bagian dari) yang sama untuk PDE elips orde dua) tidak termasuk di antara mereka - ini dikenal sebagai masalah lateral Cauchy, dan salah posisi: tidak ada jaminan bahwa pasangan data batas yang diberikan mengakui solusi, dan bahkan jika ada, tidak ada stabilitas sehubungan dengan gangguan kecil dalam data. (Sebenarnya, ini adalah masalah keliru yang asli dalam arti Hadamard, dan contoh klasik mengapa Anda tidak dapat mengabaikan kondisi batas saat membahas posisi yang baik. Anda dapat menemukan contoh eksplisit dalam kuliahnya tentang masalah Cauchy dalam diferensial parsial linear persamaan dari tahun 1920-an.)

(r,R)×(a,b)x=rRxy=ay=b

  1. Jika Anda dapat memaksakan kondisi batas (Neumann, Robin, Dirichlet - yang Anda perlukan untuk memperbaiki konstanta dalam integrasi turunan tangensial, omong-omong), maka cukup menggunakan komponen normal gradien Anda sebagai kondisi Neumann (jika Anda dapat memperbaiki mode konstan) atau mengintegrasikan komponen tangensial sebagai kondisi Dirichlet. Karena kedua kondisi tersebut mungkin berhubungan dengan fungsi yang sama, Anda juga mendapatkan solusi yang sama.

  2. y=ay=bΔu=fΔu+εΔ2u=fε>0H2uuεuε0

    H2

Christian Clason
sumber
Untuk implementasi oleh elemen campuran di FEniCS lihat biharmonicdemo. Ini mungkin tanpa istilah Laplace tapi saya rasa itu dapat dengan mudah ditambahkan.
Jan Blechta
Hai Christian, terima kasih atas saran Anda! Saya mendapat kesan bahwa persamaan Poisson tidak berbahaya sejauh menyangkut stabilitas numerik - terima kasih telah menunjukkannya. Saya akan membacanya seperti yang Anda sarankan. Tidak yakin apakah ini mengubah banyak hal, tetapi seperti yang disebutkan dalam komentar lebih lanjut, kombinasi Dirichlet-Neumann mungkin menyesatkan. 'Semua' yang saya cari, adalah cara untuk memaksakan data gradien pada batas.
Markus
2
Persamaan Poisson adalah jinak, tetapi itu bukan persamaan yang Anda coba selesaikan :) (Kondisi batas adalah bagian integral dari persamaan; operator saja tidak cukup untuk memutuskan stabilitas.)
Christian Clason
Baiklah, itu memberi saya sesuatu untuk dikunyah. Terima kasih semuanya atas waktu, saran dan kesabaran Anda - dan permintaan maaf saya karena telah jatuh ke dalam perangkap XY ...
Markus
4

Anda tidak dapat berharap bahwa solusi untuk masalah Anda yang diubah akan menjadi solusi untuk masalah Poisson karena Anda perlu mengubah masalah itu entah bagaimana untuk membuatnya berpose dengan baik.

F(u,λ)=Ω12|u|2dxΩfudxΓNgudS+ΓNλ(uuD)dS
(u,λ)V×L2(ΓN)V={vH1;v|ΓD=0}ΓDuDΓNF(u)
0=DF(u)[v]=ΩuvdxΩfvdxΓNgvdSvV,
ΓNΓD
0=DF(u,λ)[v,μ]=ΩuvdxΩfvdxΓNgvdS+ΓNλvdS+ΓNμ(uuD)dS(v,μ)V×L2(ΓN),
Δu=fun=gλΓNΓN

λ|g|

ΓDvVΓD

Kesimpulannya adalah Anda bahwa Anda tidak dapat mengharapkan bahwa PDE urutan kedua akan mengakui dua kondisi batas independen.


F(u,λ)DF(u,λ)[v,μ]derivative()

F(u,λ)λL2(ΓN)λL2(Ω)λR

Jan Blechta
sumber
2ufuH1
Sayangnya, matematika saya tidak dapat tergores dan saya tidak yakin tentang implikasi matematis dari ruang Banach, tetapi saya berjuang untuk melihat mengapa persamaan tersebut bukanlah solusi untuk persamaan Poisson ketika istilah pengali Lagrange menghilang. Dari sudut pandang fisik solusi (untuk masalah praktis yang saya jelaskan, maksud saya bukan solusi dalam pengertian matematika) harus ada sejauh yang saya bisa lihat
Markus
Jadi ini adalah masalah terbalik, menemukan kondisi batas untuk bidang jauh yang, bersama dengan kondisi Dirichlet yang dapat Anda terapkan, menghasilkan gradien normal yang diamati pada batas di mana Anda mengukur?
Markus
3

Pendekatan Anda tidak dapat bekerja, pasti karena implementasi dan mungkin karena formulasi Anda.

Memaksakan kondisi Dirichlet di dolfin, akhirnya menetapkan DOF yang sesuai dari ruang pengujian Anda menjadi nol.

Ini adalah kutipan dari fenics-manual :

Bab 3.3.9 (akhir): Penerapan kondisi batas Dirichlet ke sistem linier akan mengidentifikasi semua derajat kebebasan yang harus ditetapkan ke nilai yang diberikan dan memodifikasi sistem linier sedemikian rupa sehingga solusinya menghormati kondisi batas. Ini dicapai dengan memusatkan dan menyisipkan 1 pada diagonal dari baris matriks yang sesuai dengan nilai Dirichlet, dan memasukkan nilai Dirichlet dalam entri yang sesuai dari vektor sisi kanan [...]

vΓm

Singkatnya, menggunakan rutin standar di dolfin Anda tidak dapat menerapkan Dirichlet dan kondisi lainnya pada batas yang sama.

Namun, sebelum Anda mencoba untuk memperbaiki ini dalam implementasi Anda, cari menemukan ruang tes yang tepat untuk formulasi matematika Anda (seperti @Jan Blechta baru saja disebutkan.)

Jan
sumber
Saya mengerti maksud Anda - saya pikir formulasi saya mungkin tidak persis mencerminkan apa yang telah saya terapkan - permintaan maaf saya. Prinsip variasional hanyalah memori kabur dan saya mencoba untuk mendapatkan kepalaku lagi. Saya telah membaca manual maju dan mundur bersama-sama dengan beberapa contoh kode FEniCS yang mengimplementasikan pengganda Lagrange. Saya pikir masalah yang Anda ajukan adalah alasan mengapa Anda akan menggunakan fungsi tes kedua 'd'.
Markus
Saya setuju dengan @JanBlechta. Pada awalnya, Anda perlu menemukan ruang yang tepat untuk pengali, yang tidak trivial. Mungkin teks tentang optimisasi kendala PDE, di mana orang menggunakan pengganda untuk memasukkan kondisi sisi, akan memberikan beberapa ide bermanfaat. Dalam hal ini kertas , multiplier digunakan untuk akun untuk waktu kondisi Dirichlet tergantung.
Jan