Finite Elemente Programm < DGL < Numerik < Hochschule < Mathe < Vorhilfe
|
| Status: |
(Frage) überfällig  | | Datum: | 14:49 Sa 08.07.2006 | | Autor: | Joergi |
Hallo zusammen,
wir sind eine kleine Truppe, die die Numerik noch überlebt hat und der man jetzt so kurz vor Semesterende noch eine dicke Aufgabe "ans Herz gelegt" hat, bei der wir nun nicht mehr weiter wissen und hoffen, dass man uns helfen kann.
Wir wollen das folgende inhomogene Randwertproblem mit der Finite Element Methode (FEM) lösen:
Gegeben seien [mm]\Omega \subset \IR^{2}[/mm] und Randdaten [mm]r \in \gamma(H^{1}(\Omega)) \subset L_{2}( \partial \Omega)[/mm]. Dann suchen wir [mm]u \in H^{1}(\Omega)[/mm], so dass
[mm](\nabla u,\nabla v) = 0[/mm] [mm] \forall v \in H^{1}_{0}(\Omega)[/mm] (*)
[mm]\gamma(u) = r[/mm].
Dann lässt sich (*) wie folgt umwandeln: Sei [mm] \overline{r}\in H^{1}(\Omega)[/mm] mit [mm]\gamma(\overline{r}) = r[/mm] und [mm]u_{0}\in H^{1}_{0}(\Omega)[/mm], dann gilt:
[mm](\nabla u_{0},\nabla v) = -(\nabla \overline{r},\nabla v) [/mm] [mm] \forall v \in H^{1}_{0}(\Omega)[/mm] .
Die Lösung von (*) ergibt sich dann zu [mm]u \equiv u_{0} + \overline{r}[/mm].
Weiter seien zu [mm]\Omega[/mm] eine Triangulierung, die Knoten und die Konnektivität gegeben. Wir nehmen an, dass zusätzlich [mm]r \in C/(\partial \Omega)[/mm]. Dann approximieren wir u, [mm] u_{0},[/mm] [mm]\overline{r}[/mm] durch [mm] u_{h}, u_{h,0},[/mm] [mm] \overline{r}_{h}[/mm]:
[mm] \overline{r}_{h} \equiv \summe_{i=M_{I}+1}^{M_{N}}r(N_{i})* \phi_{i}[/mm], [mm] u_{h,0}\in V_{h,0}[/mm], so dass:
[mm](\nabla u_{h},\nabla v_{h}) = -(\nabla \overline{r_{h}},\nabla v_{h}) [/mm] (**)
[mm] \forall v_{h} \in V_{h,0}[/mm], und
[mm]u_{h} \equiv u_{h,0} + \overline{r_{h}}\in V_{h}[/mm].
a) Dann ist (**) äquivalent zu :
[mm] (\nabla u_{h,0},\nabla \phi_{i}) = -(\nabla \overline{r_{h}},\nabla \phi_{i}) [/mm] [mm] \forall i \in \{1,2, \ldots ,M_{I}\}[/mm]. (***)
Schreibe den Algorithmus aus der Vorlesung auf und modifiziere ihn so, dass die Matrix und die rechte Seite des LGS zu (***) aufgebaut werden.
b) Die Gesamtzahl der Knoten, die Anzahl der inneren Knoten und die Anzahl der Elemente seien in Matlab-Files durch [mm]M_{nodes}[/mm], [mm]M_{inner}[/mm], [mm]M_{elements}[/mm] gegeben. Weiter seien dort die x- und y-Koordinaten der Knoten in den Feldern x und y gegeben und die Konnektivität in dem Feld n. Schreibe ein Matlab-Programm, das damit das LGS aus a) aufbaut, dann das LGS von Matlab lösen lässt und schließlich [mm]u_{h}[/mm] auf der Triangulierung visualisiert. Verwende zur Visualisierung den Befehl trimesh . [Für Triangulierungen mit vielen Knoten sollte man die Matrix im sparsen Format behandeln und den passenden PCG-Löser einsetzen! Das ist hier nicht nötig.]
Als ersten Test überprüfe man, ob das Programm die exakte Lösung [mm]u(x,y) \equiv x+y[/mm] auf den gegebenen Triangulierungen findet.
Wir haben dazu Matlab-Files mit Triangulierungen erhalten!
Melements = 182;
n = zeros(182,3);
n(1,1) = 1; n(1,2) = 112; n(1,3) = 111;
n(2,1) = 112; n(2,2) = 1; n(2,3) = 2;
n(3,1) = 111; n(3,2) = 3; n(3,3) = 1;
n(4,1) = 2; n(4,2) = 1; n(4,3) = 4;
n(5,1) = 5; n(5,2) = 1; n(5,3) = 3;
n(6,1) = 112; n(6,2) = 2; n(6,3) = 6; usw.
Visualisiere nun die numerische Lösung zu den Randdaten:
[mm]r(x,y)[/mm] [mm] \equiv[/mm] [mm]\begin{cases} 0, \\ \bruch{1+y-x}{|x|+|y|}, sonst. \end{cases}[/mm]
Wir haben zunächst mal die Aufgabe a) gelöst und uns an die Vorgaben der Vorlesung gehalten:
[mm] (\nabla u_{h,0},\nabla \phi_{i}) = -(\nabla F_{h},\nabla \phi_{i}) [/mm]
[mm] \gdw a (u_{h,0}, \phi_{i}) = - a (\overline{r_{h}}, \phi_{i})[/mm]
[mm] \gdw \summe_{|\alpha|, |\beta| \le k }^{ } \integral_{\Omega}^{ }{a_{\alpha \beta}(x)*D^{\beta}u_{h,0}(x)* D^{\alpha}\phi_{i}(x) dx} = \summe_{|\alpha|, |\beta| }^{ } \integral_{\Omega}^{ }{a_{\alpha \beta}(x)*D^{\beta}\overline{r_{h}} (x)* D^{\alpha}\phi_{i}(x) dx} [/mm]
[mm] \gdw \integral_{\Omega}^{ }{a_{\alpha \beta}(x)* (\nabla u_{h,0} (x)^{T}* \nabla \phi_{i}(x) dx} = - \integral_{\Omega}^{ }{a_{\alpha \beta}(x)* \nabla \overline{r_{h}} (x)^{T}* \nabla \phi_{i}(x) dx} [/mm]
[mm] \gdw \integral_{\Omega}^{ }{a_{\alpha \beta}(x)* (\nabla (\summe_{j=1}^{M_{I}}u_{j} \phi_{j})* \nabla \phi_{i}(x) dx} = - \integral_{\Omega}^{ }{a_{\alpha \beta}(x)* \nabla (\summe_{j=M_{I}+1}^{M_{N}}r(N_{j})*\Phi_{j})}* \nabla \phi_{i}(x) dx} [/mm]
[mm] \gdw \summe_{j=1}^{M_{I}}u_{j} \integral_{\Omega}^{ }{ a_{\alpha \beta}(x)* \nabla \phi_{j}(x)^{T}\nabla \phi_{i}(x) } = - \summe_{j=M_{I}+1}^{M_{N}}r(N_{j})* \integral_{\Omega}^{ }{a_{\alpha \beta}(x) \nabla \phi_{j}(x)^{T} \nabla \phi_{i}(x) dx [/mm]
[mm] \gdw \summe_{j=1}^{M_{I}}u_{j}*a(\phi_{j},\phi_{i}) = - \summe_{j=M_{I}+1}^{M_{N}}r(N_{j})* a(\phi_{j},\phi_{i})[/mm]
[mm] \gdw \summe_{j=1}^{M_{I}}u_{j}*(\nabla \phi_{j},\nabla \phi_{i}) = - \summe_{j=M_{I}+1}^{M_{N}}r(N_{j})* (\nabla \phi_{j},\nabla \phi_{i}) [/mm]
[mm] \gdw AU = - B[/mm] wobei:
[mm]A = (a_{ij}) \in \IR^{M_{I}\times M_{I}}[/mm]
[mm] (a_{ij}) = (\nabla \phi_{j},\nabla \phi_{i}) = \integral_{\Omega}^{ }{\nabla \phi_{j}^{T}* \nabla \phi_{i}dxdy} [/mm]
[mm]U = (u_{1}, \cdots, u_{M_{I}})^{T} \in \IR^{M_{I}\times 1}[/mm], also die rechte Seite [mm] \in \IR^{M_{I}\times 1}[/mm]
[mm]B = (b_{1}, \cdots, b_{M_{I}})^{T}[/mm]
[mm] b_{i} = \summe_{j=M_{I}+1}^{M_{N}}r(N_{j})* (\nabla \phi_{j},\nabla \phi_{i})[/mm] für [mm]i \in \{1, \cdots, M_{I}\}[/mm] und
[mm] \summe_{j=1}^{M_{I}}u_{j}*(\nabla \phi_{j},\nabla \phi_{i}) = - \summe_{j=M_{I}+1}^{M_{N}}r(N_{j})* (\nabla \phi_{j},\nabla \phi_{i}) [/mm].
Nun unsere Frage. Wie sieht die rechte Seite genau aus? Wie muss das Programm aussehen, dass das LGS aufbaut!? Wir hängen also bei der b) ohne weitere Ideen völlig fest!
Wir sind echt etwas überfordert und hoffen, dass sich jemand erbarmt und gemeinsam mit uns versucht das Problemchen zu lösen.
Vielen lieben Dank im Voraus für alle Mühen, Tipps und andere Hilfen
|
|
| |
|
| Status: |
(Mitteilung) Reaktion unnötig  | | Datum: | 15:20 Fr 14.07.2006 | | Autor: | matux |
$MATUXTEXT(ueberfaellige_frage)
|
|
|
|
