FEM: Singularität der Steifheitsmatrix

{(σ^{2} (x) u^{″} (x))}^{″} = f (x), 0 ⩽ x ⩽ 1

$\left( \sigma^{2}(x) u ''(x) \right)'' = f(x), \;\;\; 0 \leqslant x \leqslant 1$

u (0) = u (1) = 0

$u(0) = u(1) = 0$

u^{″} (0) = u^{″} (1) = 0

$u''(0) = u''(1) = 0$

σ (x) ⩾ σ_{0} > 0

$\sigma(x) \geqslant \sigma_{0} > 0$

A u = f

$Au = f$

A

$A$

Nach dem FEM-Schema reduziere ich mein Problem auf ein Optimierungsproblem

J (u) = (A u, u) - 2 (f, u) \to min_{u}

$J(u) = (Au,u) - 2(f,u) \to \min_{u}$ Ich

finite Elemente

h_{k} (x)

$h_{k}(x)$ als

v_{k} (x) = {\begin{aligned} 1 - {(\frac{x - x_{k}}{h})}^{2}, & x \in [x_{k - 1}, x_{k + 1}] \\ 0, & otherwise \end{aligned}

$v_{k}(x) = \left\{ \begin{array}{rl} 1 - \left( \frac{x-x_{k}}{h} \right)^2, & x \in [x_{k-1},x_{k+1}] \\ 0, & \text{otherwise} \end{array} \right.$ für jedes

k = 1, \dots, n - 1

$k = 1,\ldots,n-1$ , wobei

x_{k} = h k

$x_{k} = hk$ ,

h = \frac{1}{n}

$h = \frac{1}{n}$ . Finite Elemente

v_{0} (x)

$v_{0}(x)$ und

v_{n} (x)

$v_{n}(x)$ werden ähnlich eingeführt.

Ich versuche, den Vektor numerisch so zu finden, $\alpha$ dass $u(x) = \sum_{k=0}^{n} \alpha_{k} v_{k}(x)$ das Optimierungsproblem löst. Wir haben

J (u) = \sum_{i = 0}^{n} \sum_{j = 0}^{n} α_{i} α_{j} (A v_{i}, v_{j}) - \sum_{i = 0}^{n} 2 α_{i} (v_{i}, f) = α^{T} V α - 2 α^{T} b \to min_{α},

$J(u) = \sum\limits_{i=0}^{n} \sum\limits_{j=0}^{n} \alpha_{i} \alpha_{j} (Av_{i},v_{j}) - \sum\limits_{i=0}^{n} 2\alpha_{i} (v_{i},f) = \alpha^{T} V \alpha - 2\alpha^{T} b \to \min\limits_{\alpha},$ wobei

b_{i} = (f, v_{i})

$b_{i} = (f,v_{i})$ und

V_{i, j} = (A v_{i}, v_{j})

$V_{i,j} = (Av_{i},v_{j})$ . Nach Differenzierung in Bezug auf

α

$\alpha$ erhalte ich

V α = b,

$V\alpha = b,$ aber hier ist die Steifheitsmatrix

V

$V$ singulär. Also, was muss ich tun? Vielleicht muss ich andere finite Elemente wählen?

finite-element Applikationen
quelle

Hallo Nimza, hast du ein Testproblem, bei dem du die genaue Lösung kennst? Wenn ja, versuchen Sie zuerst,

V^{T} V α = V^{T} b

$V^T V \alpha = V^T b$ zu lösen, um zu testen, ob Ihre Basis innerhalb der Domäne korrekt ist. Wenn alles korrekt aussieht, ist es möglicherweise das falsch gestellte BC, das die Matrix singulär macht. Aber der BC scheint mir in Ordnung zu sein.

Shuhao Cao

Antworten:

In absteigender Reihenfolge der Wahrscheinlichkeit

Falsche Basis. Aus Ihrer Beschreibung geht hervor, dass Sie genau zwei quadratische Funktionen mit Unterstützung für jedes Element haben. Dieser Raum ist keine Teilung der Einheit und nicht (kontinuierliche erste Ableitungen). Um Ihr Problem vierter Ordnung direkt zu diskretisieren (anstatt es beispielsweise auf ein System von Gleichungen zweiter Ordnung zu reduzieren), benötigen Sie eine Basis. Beachten Sie, dass die Basis alle linearen Funktionen exakt wiedergeben kann. $C^1$ $C^1$ $C^1$
Unzureichende Randbedingungen. Dies ist offensichtlich, wenn Sie den Nullraum berechnen und zeichnen.
Falsche Montage. Überprüfen Sie die Zuordnung von Elementen zur zusammengesetzten Reihenfolge, um sicherzustellen, dass sie Ihren Erwartungen entspricht, z. B. dass die Ausrichtung der Elemente nicht umgekehrt wird.
Falsche lokale Montage. In 1D können Sie analytisch berechnen, wie die Elementsteifigkeitsmatrix aussieht (möglicherweise für einen vereinfachten Fall), und überprüfen, ob der Code sie reproduziert.

Jed Brown
quelle

Vielen Dank. 1. Ich denke, dass ich eine Basis brauchen werde, weil . Wenn ich dann nur Funktionen betrachte, die Randbedingungen erfüllen, dann ist .

C^{2}

$C^2$

(A u, v) = \int_{0}^{1} σ^{2} (x) u^{″} (x) v^{″} (x) d x

$(Au,v) = \int_{0}^{1} \sigma^{2}(x) u''(x) v''(x) dx$

\ker A = {0}

$\ker A = \{ 0 \}$

Applikation

Eine Basis reicht aus, der Integrand muss nicht stetig sein. Beachten Sie, dass die Randbedingungen für zweite Ableitungen zu einem Randintegral werden. Sie können eine Basis für die direkte Diskretisierung eines Problems vierter Ordnung verwenden, müssen jedoch Sprungterme wie bei diskontinuierlichen Galerkin-Methoden für Systeme erster und zweiter Ordnung integrieren. Es ist keine schlechte Methode, aber in 1D unnötig kompliziert, weil es so einfach ist, Basen mit einer beliebigen Reihenfolge der Kontinuität (z. B. Splines) zu konstruieren. Dieses Papier ist ein Beispiel für " DG".

C^{1}

$C^1$

C^{0}

$C^0$

C^{0}

$C^0$

Jed Brown

In Ordnung. Ich habe meine Basis korrigiert: jetzt auf und . Jetzt ist es . Aber die Methode funktioniert immer noch nicht.

v_{k} (x) = \cos^{2} (\frac{π}{2 h} (x - x_{i}))

$v_{k}(x) = \cos^2 \left(\frac{\pi}{2h}(x-x_{i}) \right)$

[x_{i - 1}, x_{i + 1}]

$[x_{i-1},x_{i+1}]$

i = 1, \dots, n - 1

$i = 1,\ldots,n-1$

C^{1}

$C^1$

Applikation

Die Basis sollte in der Lage sein, lineare Funktionen zu reproduzieren, dies kann jedoch nicht. Wenn Sie das behoben haben, überprüfen Sie, ob die Integrale korrekt ausgeführt werden, und überprüfen Sie dann die Randbedingungen.

C^{1}

$C^1$

Jed Brown

Das Problem hat eindeutig eine ODD-Ordnungsableitung. Insbesondere für größere Péclet-Zahlen behält die Steifheitsmatrix möglicherweise keine "feine" Form bei, wodurch während des Zusammenbaus Nullen erzeugt werden und daher eine singuläre oder manchmal sehr kleine Determinante erhalten wird, die durch die Schwingungen im Lösungsdiagramm wahrgenommen wird.

Die Lösung für diese Art von Problem ist unter anderem die Verwendung von Strafen. Insbesondere wird dies als Petrov-Galerkin-Methode bezeichnet .

Entschuldigung für mein schlechtes Englischverständnis.

Sohail Ahmed
quelle