Matroids Matheplanet Forum Index
Moderiert von matroid
Mathematik » Numerik & Optimierung » Residualer Fehlerschätzer
Druckversion
Druckversion
Antworten
Antworten
Autor
Universität/Hochschule Residualer Fehlerschätzer
Drgglbchr
Aktiv Letzter Besuch: im letzten Monat
Dabei seit: 15.11.2019
Mitteilungen: 153
Zum letzten BeitragZum nächsten BeitragZum vorigen BeitragZum erstem Beitrag  Themenstart: 2020-05-28


hallo!
ich bin gerade dabei die finite Elemente Methode für ein L-förmiges gebiet in Matlab zu implementieren.
für den residualen Fehlerschätzer soll ich diese Funktion schreiben:
Matlab
function [eta_res_loc,eta_res_edge_loc,eta_Diri_approx]=compute_residual_error_estimate(element,face,vertex, @(x,y)volumeFunction(x,y,case_flag), @(x,y,normal)func_NormalGrad_u(x,y,normal,case_flag), u_all, gqn_f,p)

und da liegt schon mein problem - ich weiß nicht, wie ich hier vorgehen soll :(
kann mir jemand einen Hinweis geben, wie so einen Fehlerschätzer berechnet wird? ich finde im internet nicht wirklich etwas dazu :/
lg



Eine Notiz zu diese Forumbeitrag schreiben Notiz   Profil  Quote  Link auf diesen Beitrag Link
Carmageddon
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 22.12.2009
Mitteilungen: 650
Zum letzten BeitragZum nächsten BeitragZum vorigen BeitragZum erstem Beitrag  Beitrag No.1, eingetragen 2020-05-29


2020-05-28 15:51 - Drgglbchr im Themenstart schreibt:
hallo!
ich bin gerade dabei die finite Elemente Methode für ein L-förmiges gebiet in Matlab zu implementieren.
für den residualen Fehlerschätzer soll ich diese Funktion schreiben:
Matlab
function [eta_res_loc,eta_res_edge_loc,eta_Diri_approx]=compute_residual_error_estimate(element,face,vertex, @(x,y)volumeFunction(x,y,case_flag), @(x,y,normal)func_NormalGrad_u(x,y,normal,case_flag), u_all, gqn_f,p)

und da liegt schon mein problem - ich weiß nicht, wie ich hier vorgehen soll :(
kann mir jemand einen Hinweis geben, wie so einen Fehlerschätzer berechnet wird? ich finde im internet nicht wirklich etwas dazu :/
lg


Guten Morgen,

Fehlerschätzer gibt es leider (oder zum Glück) nicht nur einen, sondern einen ganzen Schwung. Sie hängt auch extrem stark von der zu lösenden partiellen Differentialgleichung ab. Ausgehend von deinen anderen Threads nehme ich an, du möchtest die Poisson Gleichung lösen.

Welchen möchtest du denn genau verwenden? Der Standard-Ansatz für einen residualen Fehlerschätzer ist normalerweise das Residuum <math>\|\nabla (u - u_h)\|</math> für welche man eine (lokal!) berechenbare Abschätzung durch partielle Integration findet.

Sucht man z.B. nach "fem residual error estimator poisson equation" findet man relativ schnell z.b.:

hier, Equation 1.108

Die konkrete Implementierung hängt natürlich stark von deiner Sprache und deinem Framework ab.



Hier hat sogar jemand eine komplette Masterarbeit darüber geschrieben. Da habe ich jetzt allerdings nicht reingeschaut.



Woran genau scheiterst du? Hast du Probleme mit der mathematischen Theorie dahinter, oder mit der konkreten Implementierung? Ausgehend von deinen anderen Threads würde ich schätzen, du hast noch Probleme mit der Theorie dahinter - denn wenn man die Theorie verstanden hat, ist eigentlich auch klar wie man so etwas zumindest grob implementieren kann. (Das ist nicht böse gemeint, aber es fällt mir persönlich schwer dir zu helfen, wenn du einfach nur ein paar Codefragmente hinwirfst, mit der Bemerkung du weißt nicht wie das oder jenes geht... Je genauer du dein Problem beschreibst, umso leichter kann man dir Hilfestellung geben.)

Viele Grüße


-----------------
Zitat: "Es gibt einen Beweis aus der Physik: Er ist kurz, er ist elegant... und falsch"



Eine Notiz zu diese Forumbeitrag schreiben Notiz   Profil  Quote  Link auf diesen Beitrag Link
Drgglbchr
Aktiv Letzter Besuch: im letzten Monat
Dabei seit: 15.11.2019
Mitteilungen: 153
Zum letzten BeitragZum nächsten BeitragZum vorigen BeitragZum erstem Beitrag  Beitrag No.2, vom Themenstarter, eingetragen 2020-05-29


hallo carmageddon!
danke für deine Antwort! :)
also ich soll ein elliptisches problem auf einem L-förmigen gebiet $(-1,1)^2/[0,1]^2$ betrachten. entlang der einspringenden ecke gelten Neumann Randbedingungen und sonst dirichlet Randbedingungen.
also $-\triangle u = 0$ bzw $- \triangle u = f$ in $\Omega$
$u = 0$ auf $\Gamma_D$
$\partial/\partial_n u = g$ auf $\Gamma_N$

der Verfeinerungsindikator, der ja dann lt Skriptum für die Berechnung des Fehlerschätzers benötigt wird, ist:
$\eta_T^2(u_h) = h^2_T \|f\|_{L^2(T)}^2 + \sum_{S \in \mathcal{S}_h, S \subset \partial T} h_S \|\left[\nabla u_h \cdot n_S \right] \|_{L^2(s)}^2$
die weitere Theorie zum Fehlerschätzer ist mir eigentlich auch ziemlich klar...also a-posteriori Fehler, Herleitung usw..

aber für mich ist es sehr schwer das ganze in einen Matlab code zu gießen. ich habe mit Matlab noch nicht so viel zu tun gehabt und tue mir mit der Implementierung schwer.
ich wäre sehr dankbar, wenn du mir ein paar Grundgedanken bzw ein Konzept für die Implementierung des residualen Fehlerschätzers geben könntest :)

Dass du mir natürlich nur anhand eines Codefragments nicht die Lösung liefern kannst ist schon klar :) ich wollte nur zeigen, welche inputparameter ich für diese Funktion zur verfügung habe ;)
und ich würde den code auch gerne selber umgesetzt bekommen, bräuchte einfach nur eine "kleine Anleitung" bezüglich vorgangsweise, wenn das möglich ist :)

Lg Drgglbchr



Eine Notiz zu diese Forumbeitrag schreiben Notiz   Profil  Quote  Link auf diesen Beitrag Link
Carmageddon
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 22.12.2009
Mitteilungen: 650
Zum letzten BeitragZum nächsten BeitragZum vorigen BeitragZum erstem Beitrag  Beitrag No.3, eingetragen 2020-05-29


2020-05-29 14:40 - Drgglbchr in Beitrag No. 2 schreibt:
hallo carmageddon!
danke für deine Antwort! :)
also ich soll ein elliptisches problem auf einem L-förmigen gebiet <math>(-1,1)^2/[0,1]^2</math> betrachten. entlang der einspringenden ecke gelten Neumann Randbedingungen und sonst dirichlet Randbedingungen.
also <math>-\triangle u = 0</math> bzw <math>- \triangle u = f</math> in <math>\Omega</math>
<math>u = 0</math> auf <math>\Gamma_D</math>
<math>\partial/\partial_n u = g</math> auf <math>\Gamma_N</math>

Das dachte ich mir schon 🙂

2020-05-29 14:40 - Drgglbchr in Beitrag No. 2 schreibt:
der Verfeinerungsindikator, der ja dann lt Skriptum für die Berechnung des Fehlerschätzers benötigt wird, ist:
<math>\eta_T^2(u_h) = h^2_T \|f\|_{L^2(T)}^2 + \sum_{S \in \mathcal{S}_h, S \subset \partial T} h_S \|\left[\nabla u_h \cdot n_S \right] \|_{L^2(s)}^2</math>
die weitere Theorie zum Fehlerschätzer ist mir eigentlich auch ziemlich klar...also a-posteriori Fehler, Herleitung usw..

Alles klar, das sieht doch schon gut aus. Da sieht man ja alles was man braucht!

2020-05-29 14:40 - Drgglbchr in Beitrag No. 2 schreibt:
aber für mich ist es sehr schwer das ganze in einen Matlab code zu gießen. ich habe mit Matlab noch nicht so viel zu tun gehabt und tue mir mit der Implementierung schwer.
ich wäre sehr dankbar, wenn du mir ein paar Grundgedanken bzw ein Konzept für die Implementierung des residualen Fehlerschätzers geben könntest :)

Dass du mir natürlich nur anhand eines Codefragments nicht die Lösung liefern kannst ist schon klar :) ich wollte nur zeigen, welche inputparameter ich für diese Funktion zur verfügung habe ;)
und ich würde den code auch gerne selber umgesetzt bekommen, bräuchte einfach nur eine "kleine Anleitung" bezüglich vorgangsweise, wenn das möglich ist :)


Der erste Schritt ist es, eine funktionierende Implementierung hinzubekommen. Der zweite Schritt ist es, das ganze effizient hinzubekommen.
Was meine ich damit? In Matlab sollte man die Verwendung von for-Schleifen möglichst vermeiden und wann immer möglich mit Matrizen rechnen. Die Implementierung von Matrix-Operationen ist richtig gut und schnell in Matlab. Diese Art zu programmieren führt allerdings manchmal zu Code mit dem ein Anfänger (erstmal) nichts anfangen kann.

Deswegen konzentrieren wir uns auf den ersten Schritt: Wir rechnen <math>\eta_T^2(u_h)</math> für ein <math>T</math> aus. Anschließend kannst du dir überlegen wie man dies zu einer Funktion erweitert, welche mit Vektoren etc rechnet.


Wir haben zwei Terme, die wir berechnen müssen:

Erster Term:
<math>h_T^2 \|f\|_{L^2(T)}^2</math>
Das solllte Standard sein. Aus deinem Element <math>T</math> berechnest du <math>h_T</math>. Bei einer Dreieckszerlegung kannst du dies direkt berechnen. Meistens nimmt man hier die längste Seite oder den Umkreisdurchmesser. Ich habe die Herleitung auch gerade nicht mehr im Kopf, was das richtige ist.
Auch das Integral ist kein Problem:

<math>f</math> ist ein Vektor (zumeist) bestehend aus Funktions-Werten an deinen Knotenpunkten. Für lineare finite Elemente sind das meist die Eckpunkte deiner Dreieckszerlegung. Für lineare Elemente kann man sich hier direkt eine Funktion schreiben, welche die Eckpunkte <math>x = (x_1, x_2, x_3)</math> des Elements <math>T</math> und die Werte <math>f = (f_1, f_2, f_3)</math> an <math>x</math> bekommt:
Matlab
function [res] = MyIntVol(x,f)

Man sieht schon, wie man hier mit Vektoren arbeiten sollte. Das Erweitern auf Matrizen <math>x</math>, <math>f</math> sollte dann leicht möglich sein.

Solltest du höhere finite Elemente verwenden, transformierst du das ganze zunächst auf ein Referenzelement und führst dort eine numerische Integration durch. Aber ich nehme mal an, du arbeitest mit linearen Elementen.



Der zweite Term ist schon ein bisschen trickreicher

<math>\sum\limits_{S \in S_h, S \subset \partial T} h_S \|\nablda u_h \cdot n_S\|^2_{L^2(S)}</math>

Hier muss nun eine Integration auf einem Linienelement durchgeführt werden. D.h. wir benötigen geometrische Informationen über das Element <math>T</math>. Ich nehme mal an <math>S_h</math> ist die Menge aller Kanten. D.h. du musst für das Gitter zunächst die entsprechenden Normalenvektoren bestimmen. Beachte, dass sie für zwei benachbarte Elemente <math>T_1,T_2</math> das Vorzeichen des Normalenvektors <math>n_S</math> für die gemeinsame Kante <math>S</math> ändert!
Üblicherweise speichert man solche Informationen in Matrizen ab, welche a-priori berechnet werden.

Der nächste Schritt ist eine Funktion zu schreiben, welche für eine gegebene Kante <math>S</math> und Normalenvektor <math>n_S</math> das Linienintegral
<math>\int_S (\nabla u_h \cdot n_S)^2 ds</math>
ausrechnet. Auch hier gilt: Einfach auf eine Integral
 <math>\int\limits_0^1 ... ds</math>
transformieren und ausrechnen.

Anschließend muss man nur noch über alle Kanten loopen und fertig ist die Funktion
Matlab
function res = MyIntBorder(nabla_u, S, n)
    res = 0;
 
    % Loop over all edges
    for i=1:nEdges
        % Compute Integral over edge
        res := res + ....
    end
end


Jetzt am Ende noch alles aufsummieren und du bist fertig.


Wie ich bereits erwähnt habe: Die Verwendung von loops sollte tunlichst vermieden werden, aber für den Anfang ist das schon okay.

Wie genau deine Grid-Informationen abgespeichert sind, kann ich dir nicht sagen, dies musst du in deinem Framework nachschauen.

Ich hoffe ich konnte einen kleinen Denkanstoß für eine Implementierung geben.

Viele Grüße





-----------------
Zitat: "Es gibt einen Beweis aus der Physik: Er ist kurz, er ist elegant... und falsch"



Eine Notiz zu diese Forumbeitrag schreiben Notiz   Profil  Quote  Link auf diesen Beitrag Link
Drgglbchr
Aktiv Letzter Besuch: im letzten Monat
Dabei seit: 15.11.2019
Mitteilungen: 153
Zum letzten BeitragZum nächsten BeitragZum vorigen BeitragZum erstem Beitrag  Beitrag No.4, vom Themenstarter, eingetragen 2020-05-31


vielen dank!
dieses Grundgerüst hilft mir wirklich weiter :)
ich werde nun versuchen es in einen code zu gießen :)



Eine Notiz zu diese Forumbeitrag schreiben Notiz   Profil  Quote  Link auf diesen Beitrag Link
Drgglbchr hat die Antworten auf ihre/seine Frage gesehen.
Neues Thema [Neues Thema] Antworten [Antworten]    Druckversion [Druckversion]

 


Wechsel in ein anderes Forum:
 Suchen    
 
All logos and trademarks in this site are property of their respective owner. The comments are property of their posters, all the rest © 2001-2020 by Matroids Matheplanet
This web site was originally made with PHP-Nuke, a former web portal system written in PHP that seems no longer to be maintained nor supported. PHP-Nuke is Free Software released under the GNU/GPL license.
Ich distanziere mich von rechtswidrigen oder anstößigen Inhalten, die sich trotz aufmerksamer Prüfung hinter hier verwendeten Links verbergen mögen.
Lesen Sie die Nutzungsbedingungen, die Distanzierung, die Datenschutzerklärung und das Impressum.
[Seitenanfang]