2d finite Differenzen

[1]:

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D

[2]:

from scipy.linalg import cho_factor, cho_solve

Wir betrachten die Poisson-Gleichung

\[-\Delta u(x,y) = f(x,y)\quad (x,y)\in D\]

mit der Randbedingung

\[u(x,y) = g(x,y)\quad (x,y) \in \partial D.\]

Um die Lösung numerisch berechnen zu können, müssen wir die Gleichung diskretisieren. Dazu wenden wir den Differenzen-Quotienten in \(x\) und \(y\) Richtung für den Laplace Operator an:

\[\Delta u(x,y) = u_{xx}(x,y) + u_{yy}(x,y) \approx \frac{u(x-h,y)-2 u(x,y) + u(x+h,y)}{h^2} + \frac{u(x,y-h)-2 u(x,y) + u(x,y+h)}{h^2}.\]

Damit folgt für die inneren Punkte eines Gebiets

\[\Delta u(x,y) \approx \Delta_h u(x,y) = \frac{u(x-h,y)+u(x,y-h)-4 u(x,y) + u(x+h,y)+u(x,y+h)}{h^2}.\]

Für Randpunkte kann die Differenz nicht berechnet werden: so ist zum Beispiel am linken Rand eines Rechteck Gebiets \((x-h,y)\not\in D\).

Um die Diskretisierung einfach zu halten, betrachten wir ein quadratisches Gebiet \(D = [0,1]^2\).

Sei

\[u_{i,j} = u(x_i, y_j),\]

wobei mit

\[x_i = \frac{i}{n},\quad y_j = \frac{j}{n}\quad \text{für}\ i,j = 0,\ldots, n\]

von einer regelmässigen Zerlegung des Einheitsquadrats ausgegangen wird.

[3]:

n = 10
h = 1/n
y,x = np.mgrid[0:1:(n+1)*1j,0:1:(n+1)*1j]

[4]:

if n == 10:
    plt.scatter(x,y)
    plt.gca().set_aspect(1)
    plt.show()
../_images/EinfuehrungPDE_Beispiel2dFDMPoisson_6_0.png

Betrachten wir den Differenzen-Quotient für den Punkt \((0.4,0.4)\), so sind die Punkte

\[(0.3,0.4), (0.4,0.4), (0.5,0.4), (0.4,0.3), (0.4,0.5)\]

involviert.

[5]:

pkt = np.array([(0.3,0.4), (0.4,0.4), (0.5,0.4), (0.4,0.3), (0.4,0.5)])

[6]:

if n == 10:
    plt.scatter(x,y)
    plt.scatter(pkt[:,0],pkt[:,1])
    plt.gca().set_aspect(1)
    plt.show()
../_images/EinfuehrungPDE_Beispiel2dFDMPoisson_9_0.png

Den Differenzenoperator \(\Delta_h\) kann man auch mit einem sogenannten Differenzenstern darstellen:

\[\begin{split}(-\Delta_h) = \frac{1}{h^2} \begin{bmatrix} & -1 & \\ -1 & 4 & -1\\ & -1 & \end{bmatrix}\end{split}\]

Um das Gleichungssystem für die inneren Punkte aufstellen zu können, muss die Nummerierung der diskretisierten Punkte festgelegt werden. Man benutzt standardmässig die natürliche Nummerierung:

[7]:

pts = np.array([(xi,yi) for xi,yi in zip(x.flatten(),y.flatten())])

Damit folgt

[8]:

if n == 10:
    plt.figure(figsize=(8,8))
    plt.scatter(x,y)
    plt.scatter(pkt[:,0],pkt[:,1])
    k=0
    for p in pts:
        plt.text(*p,str(k))
        k+=1
    plt.gca().set_aspect(1)
    plt.show()
../_images/EinfuehrungPDE_Beispiel2dFDMPoisson_14_0.png

So folgt für den Punkt \((0.4,0.4)\) die Nummer 48 und die diskrete Gleichung lautet

\[-u_{37}-u_{47}+4\, u_{48}-u_{49}-u_{59} = f(x_{48}, y_{48})\]
[9]:

stencil = np.array([4,-1,-1,-1,-1])
indstencil = np.array([0,-1,1,-n-1,n+1])

[10]:

48+indstencil

[10]:

array([48, 47, 49, 37, 59])

Wir führen nun einen Index mit allen Punkten ein:

[11]:

ind = np.arange((n+1)**2)
ind

[11]:

array([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,  10,  11,  12,
        13,  14,  15,  16,  17,  18,  19,  20,  21,  22,  23,  24,  25,
        26,  27,  28,  29,  30,  31,  32,  33,  34,  35,  36,  37,  38,
        39,  40,  41,  42,  43,  44,  45,  46,  47,  48,  49,  50,  51,
        52,  53,  54,  55,  56,  57,  58,  59,  60,  61,  62,  63,  64,
        65,  66,  67,  68,  69,  70,  71,  72,  73,  74,  75,  76,  77,
        78,  79,  80,  81,  82,  83,  84,  85,  86,  87,  88,  89,  90,
        91,  92,  93,  94,  95,  96,  97,  98,  99, 100, 101, 102, 103,
       104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116,
       117, 118, 119, 120])

Die Randpunkt sind

[12]:

indBND = [k for k in range(n+1)] # Randpunkte unten
indBND += [k for k in range(n*(n+1),(n+1)**2)] # Randpunkte oben
indBND += [k for k in range(n+1,n*(n+1),n+1)] # Randpunkte links
indBND += [k for k in range(2*(n+1)-1,(n+1)**2,n+1)] # Randpunkte rechts
indBND = np.array(indBND)
indBND

[12]:

array([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,  10, 110, 111,
       112, 113, 114, 115, 116, 117, 118, 119, 120,  11,  22,  33,  44,
        55,  66,  77,  88,  99,  21,  32,  43,  54,  65,  76,  87,  98,
       109, 120])

Für die inneren Punkte haben wir

[13]:

indInt = list(set(ind) - set(indBND))
indInt = np.array(indInt)
indInt

[13]:

array([ 12,  13,  14,  15,  16,  17,  18,  19,  20,  23,  24,  25,  26,
        27,  28,  29,  30,  31,  34,  35,  36,  37,  38,  39,  40,  41,
        42,  45,  46,  47,  48,  49,  50,  51,  52,  53,  56,  57,  58,
        59,  60,  61,  62,  63,  64,  67,  68,  69,  70,  71,  72,  73,
        74,  75,  78,  79,  80,  81,  82,  83,  84,  85,  86,  89,  90,
        91,  92,  93,  94,  95,  96,  97, 100, 101, 102, 103, 104, 105,
       106, 107, 108])

Die Anzahl Freiheitsgrade unseres Problems ist gegeben durch

[14]:

indInt.shape[0]

[14]:

81

Die Werte der Funktion \(u\) ist für alle Randpunkte gegeben.

[15]:

if n == 10:
    plt.figure(figsize=(8,8))
    plt.scatter(pts[indInt,0],pts[indInt,1])
    plt.scatter(pkt[:,0],pkt[:,1])
    plt.scatter(pts[indBND,0],pts[indBND,1])
    k=0
    for p in pts:
        plt.text(*p,str(k))
        k+=1
    plt.gca().set_aspect(1)
    plt.show()
../_images/EinfuehrungPDE_Beispiel2dFDMPoisson_27_0.png

Für alle inneren Punkte haben wir das Gleichungssystem

\[-u_{i-n}-u_{i-1}+4\, u_{i}-u_{i+1}-u_{i+n} = h^2\, f(x_{i}, y_{i})\quad i \in \text{indInt}\]

Die Anzahl Zeilen der Systemmatrix wählen wir analog zu der Anzahl Gleichungen, während die Anzahl Spalten mit der Anzahl Punkte übereinstimmt:

[16]:

A = np.zeros((indInt.shape[0],ind.shape[0]))

Wir berechnen nun die Matrix zum Gleichungssystem oben:

[17]:

for i in range(A.shape[0]): # Daher über alle Zeilen
    A[i,indInt[i]+indstencil] += stencil # Stencil wird nur bei den inneren Freiheitsgrade dazu addiert.

Diese Matrix beinhaltet alle Freiheitsgrade. Beschränken wir die Matrix auf die inneren Freiheitsgrade so erhalten wir folgendes Bild (Bandmatrix):

[18]:

Aint = A[:,indInt]

[19]:

if n == 10:
    plt.matshow(Aint)
../_images/EinfuehrungPDE_Beispiel2dFDMPoisson_35_0.png

Nun soll das Gleichungssystem für alle internen Punkte extrahiert werden.

\[A_{\text{int}}\cdot u_{\text{int}} = b\]

Dazu wählen wir nur die Zeilen der inneren Freiheitsgrade und subtrahieren sämtliche Spalten der Freiheitsgrade des Rand von der Rechtenseite.

Lösungsvektor für alle Freiheitsgrade, dh. Punkte

[20]:

u = np.zeros_like(ind,dtype=float)

Funktion für die Randwerte

[21]:

def g(x,y):
    return x*y#np.ones_like(x)

Initialisieren der Rand Freiheitsgrade

[22]:

u[indBND] = g(pts[indBND,0],pts[indBND,1])

[23]:

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(x,y,u.reshape((n+1),(n+1)))

[23]:

<mpl_toolkits.mplot3d.art3d.Poly3DCollection at 0x7fce635dad60>
../_images/EinfuehrungPDE_Beispiel2dFDMPoisson_43_1.png

Funktion der rechten Seite der Gleichung. Wir wählen

\[f(x,y) = 1 \quad\forall (x,y)\in D\]
[24]:

def f(x,y):
    #return np.sin(np.pi*x)*np.sin(np.pi*y)
    return -10*np.ones_like(x)

[25]:

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(x,y,f(x,y))

[25]:

<mpl_toolkits.mplot3d.art3d.Poly3DCollection at 0x7fce63470280>
../_images/EinfuehrungPDE_Beispiel2dFDMPoisson_46_1.png

Berechnen der rechten Seite

[26]:

b = np.zeros_like(u[indInt])
b = h**2*f(pts[indInt,0],pts[indInt,0]) # f wird für die inneren Punkte berechnet
b -= A@u # im Vektor u sind nur Freiheitsgrade der Randpunkte nicht null (abhängig von g!).

[27]:

cl = cho_factor(Aint)
u[indInt] = cho_solve(cl,b)

[28]:

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(x,y,u.reshape((n+1),(n+1)),cmap=cm.coolwarm)

[28]:

<mpl_toolkits.mplot3d.art3d.Poly3DCollection at 0x7fce634f13a0>
../_images/EinfuehrungPDE_Beispiel2dFDMPoisson_50_1.png 
[29]:

plt.contourf(x,y,u.reshape((n+1,n+1)))
plt.gca().set_aspect(1)
plt.show()
../_images/EinfuehrungPDE_Beispiel2dFDMPoisson_51_0.png

Funktioniert nur mit NGSolve

Vergleich mit Referenzlösung (FEM 4. Ordnung):

[30]:

def FEMsolution(pts, n=100,order=4):
    from netgen.geom2d import unit_square
    from ngsolve import Mesh, H1, BilinearForm, LinearForm, GridFunction, grad, x, y, dx
    from ngsolve.webgui import Draw
    ngmesh = unit_square.GenerateMesh(maxh=1/n)
    mesh = Mesh(ngmesh)
    V = H1(mesh, order=order, dirichlet='bottom|right|top|left')
    u,v = V.TnT()
    a = BilinearForm(V)
    a += grad(u)*grad(v)*dx
    a.Assemble()
    b = LinearForm(V)
    b += -10*v*dx
    b.Assemble()
    gfu = GridFunction(V)
    gfu.Set(x*y, definedon=mesh.GetBoundaries())
    b.vec.data -= a.mat*gfu.vec
    gfu.vec.data += a.mat.Inverse(freedofs = V.FreeDofs())*b.vec
    Draw(gfu)
    return np.array([gfu(mesh(*p)) for p in pts])

[31]:

np.linalg.norm(u-FEMsolution(pts,n=100))

[31]:

0.034050659459700795

[32]:

# 10, 0.034050659459700795
# 100, 0.0034719326551715442