Finite Element Methode Anwendung für Jacobi und Gauss-Seidel - Verfahren

Das Notebook dient aus Ausblick auf weiterführende numerische Methoden (im Wahlmodul HANA werden solche Themen aufgegriffen).

Das Beispiel basiert auf der C++ FEM Bibliothek NGSolve.

Im Beispiel werden die iterativen Löser

  • Jacobi - Verfahren

  • Gauss-Seidel - Verfahren

für lineare Gleichungssysteme vorgestellt.

[1]:

import netgen.geom2d as geom2d
from ngsolve import *
from ngsolve.webgui import Draw

[3]:

import numpy as np
import scipy.sparse as sp
import matplotlib.pyplot as plt

2D Geometrie

Eine klassische Geometrie für numerische Versuche ist der sogennante L-shaped domain. Das Gebiet ist aufgrund eines einspringenden Eckens nicht konvex, was in vielen Anwendungen die Konsequenz hat, dass Regularität der Lösung beeinträchtigt wird (dh. wie oft ist die Lösung stetig differenzierbar).

[4]:

geo = geom2d.SplineGeometry()
p1 = geo.AppendPoint (0,0)
p2 = geo.AppendPoint (1/2,0)
p3 = geo.AppendPoint (1/2,1/2)
p4 = geo.AppendPoint (1,1/2)
p5 = geo.AppendPoint (1,1)
p6 = geo.AppendPoint (0,1)

geo.Append (["line", p1, p2],bc='outer')
geo.Append (["line", p2, p3],bc='outer')
geo.Append (["line", p3, p4],bc='outer')
geo.Append (["line", p4, p5],bc='outer')
geo.Append (["line", p5, p6],bc='outer')
geo.Append (["line", p6, p1],bc='outer')

ngmesh = geo.GenerateMesh()
mesh = Mesh(ngmesh)

[5]:

Draw(mesh)

[5]:



[6]:

mesh.Refine()
mesh.Refine()
Draw(mesh)

[6]:

29e42793dd4b46bcbdd15273659b5dde

FEM Diskretisierung

Wir berechnen das Notebook mit zwei verschiedenen Konfigurationen:

  • order = 1 Das heist, wir benutzen nur lineare Basisfunktionen.

    • In dem Fall konvergieren beide Verfahren (Jacobi und Gauss-Seidel)

  • order = 2 (>1) Mit Hilfe von quadratisch oder höheren Basisfunktionen

    • In dem Fall konvergiert das Jacobi-Verfahren nicht mehr. Die System-Matrix genügt nicht mehr der Konvergenzbedingung.

    • Das Gauss-Seidel Verfahren funktioniert jedoch noch einwandfrei.

[7]:

order = 1

[8]:

fes = H1(mesh, order=order, dirichlet='outer')
u, v = fes.TnT()
print('ndof = ', fes.ndof)

ndof =  65

[9]:

gfu = GridFunction(fes)

Visualisierung 1. FEM Basisfunktion 1. Ordnung, daher Polynome 1. Grades

ef6a4eb2cee74cedb25b2cb05bb86efb

[10]:

gfu.vec.FV()[10] = 1

[11]:

Draw(gfu,mesh,'u')

[11]:
image.png

Matrix z.B. für Wärmeleitung

\[-\Delta u(x,y,z) = q(x,y,z)\]

analoge Verallgemeinerung der 1. Differentialgleichung

\[-u''(x) = q(x)\]

Für die Randbedingung nehmen wir homogene Dirichlet Randbedingung an. Daher für die Lösung gilt auf dem Rand \(u=0\):

\[u(\Gamma) = 0.\]

Diskretisiert, dh. mit einer endlichen Anzahl Basisfunktionen, folgt das lineare Gleichungssystem

\[A\cdot u = f.\]
[12]:

a = BilinearForm(fes)
a += grad(u)*grad(v)*dx

a.Assemble()

[12]:

<ngsolve.comp.BilinearForm at 0x15117fbb0>

[13]:

f = LinearForm(fes)
f += 10*v*dx

f.Assemble()

[13]:

<ngsolve.comp.LinearForm at 0x151180430>

Im FEM Code ist die Matrix als dünnbesetzte (sparse) Matrix gespeichert.

[13]:

rows,cols,vals = a.mat.COO()

[14]:

A = sp.csr_matrix((vals,(rows,cols)))

[15]:

plt.spy(A)
plt.show()
../_images/Gleichungen_Jacobi_Gauss-Seidel_FEM_Beispiel_27_0.png

Betrachte nur die freien Freiheitsgrade, dh. ohne Freiheitsgrade auf dem Rand:

[16]:

ind = np.arange(0,fes.ndof)[np.array(fes.FreeDofs())]
print(ind)

[ 8  9 10 14 15 21 22 23 24 25 26 27 28 29 30 31 32 36 37 39 40 41 42 43
 47 48 49 50 51 52 53 54 59]

[17]:

A

[17]:

<65x65 sparse matrix of type '<class 'numpy.float64'>'
        with 385 stored elements in Compressed Sparse Row format>

[18]:

65**2

[18]:

4225

Der Einfachheit wegen, benutzen wir nun eine vollbesetzte Matrix. (Was man natürlich in der Praxis nicht machen würde):

[19]:

B = []
for i in ind:
    B.append(A[i,ind].toarray()[0])
B = np.matrix(B)
B

[19]:

matrix([[ 4.01661837,  0.        ,  0.        , ...,  0.        ,
          0.        ,  0.        ],
        [ 0.        ,  4.0184236 ,  0.        , ..., -1.10325656,
          0.        ,  0.        ],
        [ 0.        ,  0.        ,  4.02997533, ...,  0.        ,
         -0.92242727,  0.        ],
        ...,
        [ 0.        , -1.10325656,  0.        , ...,  4.01932808,
          0.        ,  0.        ],
        [ 0.        ,  0.        , -0.92242727, ...,  0.        ,
          4.03203957,  0.        ],
        [ 0.        ,  0.        ,  0.        , ...,  0.        ,
          0.        ,  4.0305739 ]])

Rechte Seite für die freien Freiheitsgrade.

[20]:

fi=f.vec.FV().NumPy()[ind]

Die diskrete FEM Lösung der partellen Differentialgleichung ist gegeben durch die Lösung des linearen Gleichungssystems

\[B\cdot u = f\]

Das System wird nun mit den beiden Verfahren im Vergleich gelöst.

Jacobi Verfahren

Additive Zerlegung der vollbesetzten Matrix \(B\)

\[B = D-L-R\]
[21]:

d = np.diag(B)
D=np.diag(d)
L = -np.tril(B,-1)
R = -np.triu(B,1)

Kontrolle:

[22]:

np.linalg.norm(B-D+L+R,np.inf)

[22]:

0.0

Kontrolle: \(2D-B\) muss positiv definit sein

[23]:

min(np.linalg.eigvals(2*D-B))

[23]:

0.6065746969243466

Jacobi-Verfahren

\[D\cdot u_{k+1} = (L+R)\cdot u_k + b\]
[24]:

u0 = np.ones_like(fi)

[25]:

LR=L+R

[26]:

u1 = u0

tol = 1e-14
nIter = 300
k=0
res = 1
while res > tol and k < nIter:
    u2 = (LR.dot(u1)+fi)/d # sehr schnell, da nur eine Matrix Multiplikation
    res = np.linalg.norm(u2-u1,ord=np.inf)
    u1 = u2
    k += 1
    print(k,res)

1 0.4623395568427219
2 0.20894697354993874
3 0.14581905429998843
4 0.10427134839102636
5 0.08796066059054819
6 0.07390082493968464
7 0.0645655628700903
8 0.05629479691288519
9 0.04816904647814191
10 0.041744262201301896
11 0.035467215344990866
12 0.030525605181523807
13 0.025874051987301583
14 0.02214340740053322
15 0.01875850132572715
16 0.015983671972249758
17 0.013541958279954192
18 0.011500668494926203
19 0.009747289879043275
20 0.008257501476687579
21 0.007001524475260146
22 0.005920396592704269
23 0.005021972807222497
24 0.0042405929264732745
25 0.0035984413810953297
26 0.003035354734762241
27 0.0025765729491577516
28 0.002171655797589711
29 0.0018439514829060588
30 0.0015532257231091195
31 0.001319167311659708
32 0.001110668014623728
33 0.0009434940916211443
34 0.000794091932958152
35 0.000674683068944415
36 0.0005676956665418564
37 0.0004823973398002668
38 0.0004058201746827139
39 0.0003448824196931599
40 0.0002900917755961907
41 0.00024655296272124705
42 0.00020736155503819909
43 0.00017625078932620442
44 0.00014822355161964174
45 0.0001259911102753919
46 0.00010595120602763242
47 9.006195177674314e-05
48 7.573502924562714e-05
49 6.437820308208275e-05
50 5.413668151332818e-05
51 4.601876370757818e-05
52 3.869827596753028e-05
53 3.289512723358978e-05
54 2.7662893561430923e-05
55 2.3514198901675165e-05
56 1.9774713089004692e-05
57 1.680861494213204e-05
58 1.4136106008144633e-05
59 1.2015384545815966e-05
60 1.010547452406918e-05
61 8.589109181444954e-06
62 7.224223300794641e-06
63 6.139935828874954e-06
64 5.16455970400731e-06
65 4.389198822229989e-06
66 3.692183307257846e-06
67 3.1377097086116024e-06
68 2.6396181895260717e-06
69 2.2430898254932607e-06
70 1.8871523120855827e-06
71 1.6035674361702945e-06
72 1.3492140150672682e-06
73 1.1463960344904045e-06
74 9.646346154079666e-07
75 8.195759397811031e-07
76 6.896885300089473e-07
77 5.859370696414956e-07
78 4.931185945133976e-07
79 4.1890944252331863e-07
80 3.5258028158180466e-07
81 2.9950002361500694e-07
82 2.52100076092443e-07
83 2.1413183126650637e-07
84 1.802587265009592e-07
85 1.5309935902374505e-07
86 1.2889261541415564e-07
87 1.0946452277948993e-07
88 9.216546953405569e-08
89 7.826748421591034e-08
90 6.590478712764636e-08
91 5.5962540435050556e-08
92 4.712748380431009e-08
93 4.0014888735839804e-08
94 3.370079654230196e-08
95 2.861238312945602e-08
96 2.409987293061633e-08
97 2.045948765960759e-08
98 1.7234475402361227e-08
99 1.4629984002034746e-08
100 1.2325091891529638e-08
101 1.0461679245121758e-08
102 8.81436640343125e-09
103 7.481134545184176e-09
104 6.303778110616776e-09
105 5.349856169178935e-09
106 4.508371398426192e-09
107 3.825827765702172e-09
108 3.2243891445737916e-09
109 2.7360083643124256e-09
110 2.3061325626905216e-09
111 1.956672934788628e-09
112 1.6494159926772056e-09
113 1.3993548542146073e-09
114 1.1797373633548602e-09
115 1.0007979445170179e-09
116 8.438196252846808e-10
117 7.157706649607576e-10
118 6.035638100421181e-10
119 5.119298318589927e-10
120 4.3172398989099747e-10
121 3.6614750120733675e-10
122 3.088150846153326e-10
123 2.6188512469715874e-10
124 2.2090224094384325e-10
125 1.873159405363367e-10
126 1.580197639405867e-10
127 1.3398238074557867e-10
128 1.1304002178746941e-10
129 9.583644988708784e-11
130 8.086542546692499e-11
131 6.855244150116846e-11
132 5.784994705493318e-11
133 4.903699668545869e-11
134 4.1385950222405654e-11
135 3.507794055224167e-11
136 2.960820477682091e-11
137 2.5093149780275326e-11
138 2.118277775409183e-11
139 1.7950918529408e-11
140 1.5155210419948162e-11
141 1.2841783192385492e-11
142 1.0843104192304054e-11
143 9.187040017621939e-12
144 7.757960940324438e-12
145 6.572520305780927e-12
146 5.5508375673696264e-12
147 4.702183087346157e-12
148 3.971711848294035e-12
149 3.3641978092191493e-12
150 2.8418933872842445e-12
151 2.406908006236108e-12
152 2.0335955142058992e-12
153 1.7221224446473116e-12
154 1.4551138072249614e-12
155 1.2322365350314612e-12
156 1.0413336859471656e-12
157 8.817391261572993e-13
158 7.451261829771738e-13
159 6.308842337432452e-13
160 5.332401187274627e-13
161 4.5141668181258865e-13
162 3.815836535636663e-13
163 3.230193890146893e-13
164 2.731148640577885e-13
165 2.3120394487818885e-13
166 1.9551027463649007e-13
167 1.6542323066914832e-13
168 1.3994361225400098e-13
169 1.1840528557627295e-13
170 1.0014211682118912e-13
171 8.471001677889944e-14
172 7.166489623955385e-14
173 6.056266599330229e-14
174 5.129230373768223e-14
175 4.3298697960381105e-14
176 3.6692870963861424e-14
177 3.108624468950438e-14
178 2.6256774532384952e-14
179 2.225997164373439e-14
180 1.887379141862766e-14
181 1.587618925213974e-14
182 1.3489209749195652e-14
183 1.149080830487037e-14
184 9.71445146547012e-15

[27]:

uset = np.zeros(fes.ndof)
uset[ind]=u2
gfu.vec.FV().NumPy()[:] = uset

[28]:

Draw(gfu,mesh,'u')

[28]:

Vergleich mit sparse Cholesky Solver:

[29]:

gfu1 = GridFunction(fes)
gfu1.vec.data = a.mat.Inverse(freedofs = fes.FreeDofs(),inverse='sparsecholesky')*f.vec

[30]:

Draw(gfu1,mesh,'u1')

[30]:

db0256fc679f4c76975ada3aeebd55be

[31]:

ufem = gfu1.vec.FV().NumPy()[ind]

[32]:

plt.plot(u2-ufem,'--')
plt.show()
../_images/Gleichungen_Jacobi_Gauss-Seidel_FEM_Beispiel_54_0.png

Gauss-Seidel Verfahren

Gauss-Seidel Verfahren

\[D\cdot u_{k+1} = L\cdot u_{k\color{red}+\color{red}1} + R\cdot u_k + b\]

Damit folgt:

\[(D-L)\cdot u_{k+1} = R\cdot u_k + b\]
[33]:

DL = D-L

[34]:

def forward(A, b):
    x = np.zeros_like(b) # Lösungsvektor
    x[0] = b[0]/A[0,0]
    for k in range(1,len(b)):
        x[k] = (b[k]-A[k,:k].dot(x[:k]))/A[k,k]
    return x

[35]:

u1 = u0

tol = 1e-14
nIter = 300
k=0
res = 1
while res > tol and k < nIter:
    u2 = forward(DL,R.dot(u1)+fi)
    res = np.linalg.norm(u2-u1,ord=np.inf)
    u1 = u2
    k += 1
    print(k,res)

1 0.4617640583386501
2 0.3122143418744079
3 0.18366149241260532
4 0.11779585960680106
5 0.08867277978823018
6 0.06641312161397506
7 0.048772809847041465
8 0.03543837249387738
9 0.025588499051642455
10 0.01840393316927824
11 0.013202724688611
12 0.009455183585547589
13 0.00676343026030346
14 0.004834070599857743
15 0.0034531509772105506
16 0.0024657491513928598
17 0.0017602095637974657
18 0.0012563128792726919
19 0.0008965492930213537
20 0.000639750886972823
21 0.0004564783277971052
22 0.00032569453565517437
23 0.00023237411491827675
24 0.00016578915431364472
25 0.00011828193987850621
26 8.43872235185561e-05
27 6.020495864478281e-05
28 4.295226498313154e-05
29 3.064352784071245e-05
30 2.1862044369547196e-05
31 1.55970472752176e-05
32 1.1127402693267285e-05
33 7.938622402614559e-06
34 5.6636510446517185e-06
35 4.0406187778452285e-06
36 2.8826993039854365e-06
37 2.056605042599635e-06
38 1.4672446712604703e-06
39 1.0467772503797157e-06
40 7.468030715496177e-07
41 5.327923669473655e-07
42 3.8011058595133207e-07
43 2.711827117862775e-07
44 1.9347020463467146e-07
45 1.3802769410764526e-07
46 9.847328347190043e-08
47 7.025393311543837e-08
48 5.012136627957631e-08
49 3.5758162097998536e-08
50 2.551100103387327e-08
51 1.820035394617392e-08
52 1.2984708031016368e-08
53 9.26370152809497e-09
54 6.609017999537059e-09
55 4.715082879691579e-09
56 3.3638896113075134e-09
57 2.3999056075751923e-09
58 1.7121687401200347e-09
59 1.221515555371866e-09
60 8.714679533561309e-10
61 6.217329984181674e-10
62 4.435639078259612e-10
63 3.164527528909389e-10
64 2.2576740477120438e-10
65 1.6106960210038324e-10
66 1.1491213536274358e-10
67 8.198203227394174e-11
68 5.848854733869757e-11
69 4.1727621358234046e-11
70 2.976974222690387e-11
71 2.12388440168354e-11
72 1.51524348623866e-11
73 1.0810075057321455e-11
74 7.712330774012344e-12
75 5.502265310042276e-12
76 3.9255265704696285e-12
77 2.800482068465726e-12
78 1.9980683774178942e-12
79 1.4253598301650072e-12
80 1.0170198017078746e-12
81 7.254752354413085e-13
82 5.17585974080248e-13
83 3.6926017799032707e-13
84 2.635114348947809e-13
85 1.8801626922027026e-13
86 1.3405943022348765e-13
87 9.564571357145724e-14
88 6.822320486321587e-14
89 4.8683279629813114e-14
90 3.480549182199866e-14
91 2.481348460037225e-14
92 1.765254609153999e-14
93 1.2656542480726785e-14
94 8.93729534823251e-15

[36]:

uset = np.zeros(fes.ndof)
uset[ind]=u2
gfu.vec.FV().NumPy()[:] = uset

Vergleich mit sparse Cholesky Solver:

[37]:

plt.plot(u2-ufem,'--')
plt.show()
../_images/Gleichungen_Jacobi_Gauss-Seidel_FEM_Beispiel_62_0.png

Bemerkung: Die Matrix ist für FEM Ansatzfunktionen höherer Ordnung (\(p\ge 2\)) nicht mehr diagonal dominant. Das Jacobi-Verfahren konvergiert daher nicht. Das Gauss-Seidel Vefahren jedoch schon.