ritz-galerkin.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "# Wprowadzenie do obliczeń symbolicznych i numerycznych w Pythonie\n",
    "Będziemy potrzebowali bibliotek `numpy` oraz `sympy`. Jeżeli ich nie masz, to zainstaluj je za pomocą\n",
    "menedżera modułów (pip lub Conda w zależności czego używasz).\n",
    "\n",
    "Pierwsza to biblioteka do obliczeń numerycznych, druga to biblioteka do obliczeń symbolicznych. Obydwie przydadzą\n",
    "się nam do rozwiązywania równań różniczkowych metodą Ritza oraz Galerkina.\n",
    "\n",
    "Na wstępie przyjrzymy się modułowi do obliczeń symbolicznych."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "outputs": [],
   "source": [
    "from sympy import *\n",
    "# definijemy symbole\n",
    "x, y = symbols(\"x y\")"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n",
     "is_executing": false
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "outputs": [
    {
     "name": "stdout",
     "text": [
      "x**2*y**2/(x + y)**2\n",
      "Pochodna z z po x to:\n",
      "-2*x**2*y**2/(x + y)**3 + 2*x*y**2/(x + y)**2\n",
      "Całka z po y to:\n",
      "x**2*(-x**2/(x + y) - 2*x*log(x + y) + y)\n",
      "Całka podwójna oznaczona y od 1 do 2, x od 0 do 1:\n",
      "-17*log(3)/2 + 13/4 + 9*log(2)\n",
      "Wynik numeryczny:\n",
      "0.150120171360575\n"
     ],
     "output_type": "stream"
    }
   ],
   "source": [
    "z = (x * y / (x + y)) ** 2\n",
    "print(z)\n",
    "print('Pochodna z z po x to:')\n",
    "print(diff(z, x))\n",
    "print('Całka z po y to:')\n",
    "print(integrate(z, y))\n",
    "print('Całka podwójna oznaczona y od 1 do 2, x od 0 do 1:')\n",
    "result = integrate(z, (y, 1, 2), (x, 0, 1))\n",
    "print(result)\n",
    "print('Wynik numeryczny:')\n",
    "print(N(result))"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n",
     "is_executing": false
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Zadanie 1.\n",
    "Oblicz pochodną cząstkową funkcji $f(x,y)=\\frac{ln(x+1) \\cdot y}{\\left(x+y\\right)^2}$ po zmiennej x w punkcie $(x,y)=(4,4)$.\n",
    "\n",
    "Wskazówka: podstawienie - https://docs.sympy.org/latest/tutorial/basic_operations.html"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "outputs": [],
   "source": [],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n",
     "is_executing": false
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Rozwiązywanie układów równań liniowych\n",
    "Układ rówań liniowych można przedstawić w postaci równania macierzowego:\n",
    "\n",
    "$Ax=b$\n",
    "\n",
    "Zapewne znane są wam metody rozwiązywania układów równań. Poznamy tutaj tylko sposób jak je wywołać w bibliotece sympy."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "outputs": [
    {
     "name": "stdout",
     "text": [
      "{(-3, 3)}\n",
      "[-3.  3.]\n"
     ],
     "output_type": "stream"
    }
   ],
   "source": [
    "import numpy as np\n",
    "# rozwiązanie symboliczne za pomocą sympy\n",
    "system = Matrix(([1, 2, 3], [3, 3, 0]))\n",
    "print(linsolve(system, (x, y)))\n",
    "# rozwiązanie numeryczne za pomocą numpy\n",
    "A = np.array([[1,2], [3,3]])\n",
    "b = np.array([3,0])\n",
    "print(np.linalg.solve(A, b))"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n",
     "is_executing": false
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Numeryczne rozwiązywanie równań różniczkowych\n",
    "Mamy następujący problem:\n",
    "\n",
    "$Au=f$ w $D \\subset R^n$\n",
    "\n",
    "$u|_{\\partial D}=0$ - jednorodne warunki brzegowe\n",
    "\n",
    "W wypadku omawianych dzisiaj metod (Ritza oraz Galerkina) będziemy potrzebowali funkcji bazowych (zwanych też funkcjami kształtu, funkcjami podstawowymi).\n",
    "Następnie będziemy szukali takiej ich kombinacji liniowej, aby jak najlepiej przybliżyć rozwiązanie, czyli poszukujemy takich $\\left\\lbrace \\alpha_j\\right\\rbrace$, aby $\\left\\lbrace z_j\\right\\rbrace$: \n",
    "\n",
    "$v=\\sum_{i=1}^{n}\\alpha_jz_j$ najlepiej przybliżał analityczne rozwiązanie.\n",
    "\n",
    "Należy zwrócić uwagę na to, że podane metody nie mogą rozwiązać dowolnego równania różniczkowego czątkowego - to równanie musi się przedstawić w postaci dywergencyjnej (wykład 5).\n",
    "Przykłady w tym notebooku będą oczywiście spełniały podane założenia."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Przykład - Metoda Ritza / Reyleigha-Ritza\n",
    "Mamy dane równanie różniczkowe:\n",
    "\n",
    "$-u''(x)=-2$ dla $x\\in (0,1)$\n",
    "\n",
    "Warunki brzegowe: \n",
    "$u(0)=u(1)=0$\n",
    "\n",
    "Wykorzystamy następujące funkcje bazowe: $z_1(x)=x(1-x)$ oraz $z_2(x)=x^2(1-x).$\n",
    "\n",
    "Mamy równanie postaci $Au=f$. Aby uzyskać odpowiednie współczynniki $\\alpha$ musimy zminimalizować funkcjonał:\n",
    "$(Av,v)-2(f,v)$, co odpowiada funkcjonałowi:\n",
    "\n",
    "$F(v)=-\\int_0^1v''(x)v(x)dx +\\int_0^14v(x)dx = \\int_0^1\\left( u'(x)\\right)^2+4v(x)dx$"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "outputs": [
    {
     "name": "stdout",
     "text": [
      "[ 3.         -2.33333333]\n"
     ],
     "output_type": "stream"
    }
   ],
   "source": [
    "x, a1, a2 = symbols(\"x a1 a2\")\n",
    "z1 = x ** 2 * (1-x)\n",
    "z2 = x ** 3 * (1-x)\n",
    "v = a1 * z1 + a2 * z2\n",
    "F = integrate(diff(v,x) ** 2 + 4 * v, (x, 0, 1))\n",
    "\n",
    "r1 = [float(x) for x in Poly(diff(F, a1),(a1,a2)).coeffs()]\n",
    "r2 = [float(x) for x in Poly(diff(F, a2),(a1,a2)).coeffs()]\n",
    "\n",
    "# a1, a2, wyraz wolny\n",
    "A = np.array([r1[:-1], r2[:-1]])\n",
    "b = np.array([r1[-1], r2[-1]])\n",
    "print(np.linalg.solve(A, b))"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n",
     "is_executing": false
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Zadanie 2.\n",
    "Rozwiąż równanie analitycznie. Sporządź wykres błędu (różnicy między rozwiązaniem dokładnym, a przybliżonym).\n",
    "\n",
    "https://docs.sympy.org/latest/modules/plotting.html"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n",
     "is_executing": false
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Zadanie 3.\n",
    "Wypróbuj inny zestaw funkcji bazowych. Sporządź wykres zależności liczby funkcji bazowych w zależności od błędu \n",
    "aproksymacji dla liczby funkcji od 1 do 5. Niech funkcje będą postaci: $z_n=x^{n+1}(x-1)$.  "
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Metoda Galerkina\n",
    "Mamy dane równanie różniczkowe: $Au=f$.\n",
    "\n",
    "Przemnóżmy skalarnie to równanie przez u: $(Au,u)=(f,u)$.\n",
    "\n",
    "Wstawmy w miejce u nasze funkcje bazowe - stwórzmy macierz Grama G: \n",
    "- $(Az_1,z_1) = g_{11}$\n",
    "- $(Az_1,z_2) = g_{12}$\n",
    "- $(Az_m,z_n) = g_{mn}$\n",
    "- ...\n",
    "\n",
    "Aby uzyskać współczynniki $\\alpha$ przy funkcjach bazowych naszego przyblożonego rozwiązania należy rozwiązać układ $G\\alpha = b$.\n",
    "\n",
    "$\\alpha=[\\alpha_1,\\alpha_2,$...$]$\n",
    "\n",
    "Natomiast wektor wyrazów wolnych będzie składał się z iloczynów: $b=[(f,z_1),(f,z_2),$...$]$"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "outputs": [
    {
     "name": "stdout",
     "text": [
      "[ 3.         -2.33333333]\n"
     ],
     "output_type": "stream"
    }
   ],
   "source": [
    "def my_norm(a, b):\n",
    "    return integrate(diff(diff(a, x), x) * b, (x, 0, 1))\n",
    "base_functions = [z1, z2]\n",
    "G = np.array([[float(my_norm(x, y)) for x in base_functions] for y in base_functions])\n",
    "b = np.array([float(integrate(-2 * z, (x, 0, 1))) for z in base_functions])\n",
    "print(np.linalg.solve(G,b))"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n",
     "is_executing": false
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Jak można zauważyć otrzymujemy dokładnie takie samo rozwiązanie, jak to z wykorzystaniem metody Ritza.\n",
    "## Zadanie 4\n",
    "Wykorzystując metodę Galerkina rozwiąż równanie różniczkowe:\n",
    "\n",
    "$\\frac{d^2T}{dx^2} + 1000x^2 = 0$, $x\\in[0,1]$\n",
    "\n",
    "Przy następujących warunkach brzegowych: $T(0)=T(1)=0$\n",
    "\n",
    "Analizowane równanie jest formą równania przewodnictwa cieplnego w izolowanym pręcie, z wewnętrznym źródłem energii. Wydajność tego źródła jest proporcjonalna do $x^2$."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Zadanie 5\n",
    "Wykorzystując metodę Ritza rozwiąż równanie różniczkowe:\n",
    "\n",
    "$\\Delta u=1$ w obszarze $D = [-1,1]\\times[-1,1]$\n",
    "\n",
    "$u|_{\\partial D}=0$\n",
    "\n",
    "Jaka jest postać funkcjonału, który będziemy minimalizować? Jaka jest postać rozwiązania przy wykorzystaniu jednej funkcji kształtu $z(x,y)=(1-x^2)(1-y^2)$?\n",
    "\n"
   ],
   "metadata": {
    "collapsed": false
   }
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  },
  "pycharm": {
   "stem_cell": {
    "cell_type": "raw",
    "source": [],
    "metadata": {
     "collapsed": false
    }
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}