{ "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.5.4" }, "name": "" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Finite Differenzen\n", "\n", "Ziel dieser Aufgabe ist es, das Poissonproblem\n", "\n", " -\u0394u = 20\u03c0\u00b2 sin(10\u03c0x) cos(20\u03c0y) in \u03a9 = (0,1)\u00b2\n", " u(x,y) = u\u2080(x) = sin(\u03c0x) auf \u2202\u03a9\n", "\n", "mit dem Finite-Differenzen-Verfahren auf einem 128x128-Rechtecksgitter zu l\u00f6sen.\n", "Die Gleichung soll dabei mit dem F\u00fcnf-Punkte-Stern\n", "\n", " $\\frac{1}{h^2} \\begin{bmatrix} & -1 & \\\\ -1 & 4 & -1 \\\\ & -1 & \\end{bmatrix}$\n", "\n", "diskretisiert werden.\n", "\n", "Die einzelnen Teilaufgaben behandeln jeweils Teile des kompletten Programms.\n", "\n", "*Hinweis:* Um Ihr Programm zu testen, k\u00f6nnen Sie auch zun\u00e4chst auch das Problem\n", "\n", " -\u0394u = 0 in \u03a9 = (0,1)\u00b2\n", " u(x,y) = x auf \u2202\u03a9\n", "\n", "betrachten.\n", "Dieses hat die einfache L\u00f6sung u(x,y) = x.\n", "\n", "*Hinweis:* Vor der Abgabe entfernen Sie bitte die Programmausgabe (Cell \u2192 All Output \u2192 Clear)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Ihr Programm sollte mit den folgenden Import-Anweisungen auskommen:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "%matplotlib inline\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import scipy.sparse as sparse\n", "from scipy.sparse.linalg import spsolve" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "Im Folgenden wird das Problem definiert und das Gitter erzeugt:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "n = 128\n", "f = lambda x: 20 * np.pi**2 * np.sin(10*np.pi*x[0]) * np.cos(20*np.pi*x[1])\n", "u0 = lambda x: np.sin(np.pi * x[0])\n", "\n", "xs, h = np.linspace(0, 1, n, retstep=True)\n", "ys = np.linspace(0, 1, n)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "`N` ist die Anzahl der Freiheitsgrade, einschlie\u00dflich der durch die Randbedingung festgehaltenen Punkte:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "N = n*n\n", "N" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Aufgabe 1\n", "\n", "Implementieren Sie eine Funktion `index(i,j)`, die jedem Paar (i,j) mit 0 \u2264 `i`, `j` < `n` bijektiv ein n(i,j) \u2208 {0,...,N-1} zuweist.\n", "Hierdurch werden die Freiheitsgrade durchgehend nummeriert.\n", "\n", "Implementieren Sie eine Funktion `coord(i,j)`, die jedem Paar (i,j) mit 0 \u2264 `i`, `j` < `n` die zugeh\u00f6rige Koordinate `x(i,j)` \u2208 \u211d\u00b2 zuweist." ] }, { "cell_type": "code", "collapsed": true, "input": [ "def index(i, j):\n", " # Hier muss Ihre Antwort hin\n", " pass\n", "\n", "def coord(i, j):\n", " # Hier muss Ihre Antwort hin\n", " pass" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Aufgabe 2\n", "\n", "Legen Sie eine d\u00fcnn besetzte Matrix `A` als `lil_matrix` geeigneter Gr\u00f6\u00dfe an; legen Sie einen Vektor `b` geeigneter Gr\u00f6\u00dfe an und initialisieren Sie ihn mit 0.\n", "\n", "Iterieren Sie \u00fcber die inneren Zeilen `i` und inneren Spalten `j` und berechnen Sie die zum Freiheitsgrad n(i,j) geh\u00f6rige Zeile der Systemmatrix `A` und der rechten Seite `b`.\n", "Benutzen Sie daf\u00fcr die Methode `index` und `coord` aus Aufgabe 1." ] }, { "cell_type": "code", "collapsed": false, "input": [ "# A = ...\n", "# b = ..." ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": null }, { "cell_type": "code", "collapsed": false, "input": [ "# for ..." ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Aufgabe 3\n", "\n", "Implementieren Sie eine Funktion `dirichlet(i,j)`, die in die zum n(i,j)-ten Freiheitsgrad geh\u00f6rige Zeile der Systemmatrix eine 1 auf die Diagonale schreibt und den n(i,j)-ten Eintrag der rechten Seite mit dem zugeh\u00f6rigen Randdaten f\u00fcllt.\n", "\n", "Iterieren Sie \u00fcber die \u00e4u\u00dferen Knoten und rufen Sie f\u00fcr jeden Knoten die `dirichlet`-Funktion auf." ] }, { "cell_type": "code", "collapsed": false, "input": [ "def dirichlet(i, j):\n", " # Hier muss Ihre Antwort hin\n", " pass" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": null }, { "cell_type": "code", "collapsed": false, "input": [ "# for ..." ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Aufgabe 4\n", "\n", "Lassen Sie sich den Speicherbedarf der d\u00fcnn besetzten Matrix ausgeben und vergleichen Sie mit dem Speicherbedarf einer dicht besetzten Matrix.\n", "Es gen\u00fcgt nur den Speicherbedarf f\u00fcr die Zahlen der Matrixeintr\u00e4ge selbst zu ber\u00fccksichtigen (8 Bytes f\u00fcr eine Flie\u00dfkommazahl mit doppelter Genauigkeit).\n", "\n", "Wie steigt der Speicherbedarf in Abh\u00e4ngigkeit von `n`?\n", "\n", "*Hinweis:* `A.nnz` liefert die Anzahl der Nicht-Null-Eintr\u00e4ge einer d\u00fcnn besetzen Matrix (**n**umber of **n**on-**z**eroes)." ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Hier muss Ihre Antwort hin" ], "language": "python", "metadata": { "scrolled": false }, "outputs": [], "prompt_number": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "(Schreiben Sie hier Ihre Antwort.)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Aufgabe 5 (optional)\n", "\n", "Benutzen Sie den foldenden programmcode, um f\u00fcr `n \u2264 32` die Besetzung der Systemmatrix `A` anzuschauen.\n", "K\u00f6nnen Sie eine Struktur erkennen?" ] }, { "cell_type": "code", "collapsed": false, "input": [ "if N <= 1024:\n", " Adense = np.array(abs(A).todense())\n", " plt.figure(figsize=(10,10))\n", " plt.pcolormesh(Adense, cmap=plt.cm.gray_r)\n", " plt.show()" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "(Schreiben Sie hier Ihre Antwort.)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Aufgabe 6\n", "\n", "Konvertieren Sie `A` in eine `csr_matrix` und l\u00f6sen Sie das Gleichungssystem `Ax = b`." ] }, { "cell_type": "code", "collapsed": false, "input": [ "A = A.tocsr()\n", "x = spsolve(A, b)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Aufgabe 7 (optional)\n", "\n", "Benutzen Sie den folgenden Programmcode um die L\u00f6sung zu betrachten.\n", "Entspricht Sie Ihren Erwartungen?" ] }, { "cell_type": "code", "collapsed": false, "input": [ "zs = np.zeros((n,n))\n", "for i in range(0,n):\n", " for j in range(0,n):\n", " zs[i,j] = x[index(i,j)]\n", "\n", "plt.figure(figsize=(10,10))\n", "plt.pcolormesh(ys, xs, zs)\n", "plt.show()" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "(Schreiben Sie hier Ihre Antwort.)" ] } ], "metadata": {} } ] }