manufactured_solution.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline\n",
    "\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "import sympy as sym"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solving \n",
    "\n",
    "$$\n",
    "-\\Delta u = f\n",
    "$$\n",
    "\n",
    "Example 1 (peak 2d): \n",
    "\n",
    "$$\n",
    "u(x,y) = x(x-1)y(y-1)e^{-100\\big( (x-0.5)^2+(y-0.117)^2\\big)}\n",
    "$$\n",
    "\n",
    "Example 2 (for L-shaped domain, no-rhs, only boundary conditions, in polar coordinates):\n",
    "\n",
    "$$\n",
    "u(r,\\theta) = r^{2/3} \\sin\\left(\\frac{2\\theta+5\\pi}{3}\\right)\n",
    "$$\n",
    "\n",
    "etc.\n",
    "\n",
    "Play around with your own examples!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sympy import atan2, sin, cos\n",
    "\n",
    "def u(x,y):\n",
    "    return x*(x-1)*y*(y-1)*exp(-100*(x-0.5)**2-100*(y-0.117)**2)\n",
    "\n",
    "def u1(x,y):\n",
    "    r = (x*x+y*y)**(1/2)\n",
    "    theta = atan2(x,y)\n",
    "    return r**(2/3)*sin(theta*2/3)\n",
    "\n",
    "def u2(x,y,z):\n",
    "    return x*(x-1)*y*(y-1)*z*(z-1)*exp(-100*(x-0.5)**2-100*(y-0.117)**2-100*(z-0.331)**2)\n",
    "\n",
    "def u3(x,y,z):\n",
    "    return (x*x + y*y + z*z)**(1/4);\n",
    "\n",
    "def u4(x,y):\n",
    "    return sin(2*sym.pi*y)*sin(2*sym.pi*x)\n",
    "\n",
    "def u5(x,y):\n",
    "    return x*(x-1)*y*(y-1)*(10*exp(-100*(x-0.5)**2-100*(y-0.1)**2)+.1*exp(-100*(x-0.1)**2-100*(y-0.5)**2))\n",
    "\n",
    "all_2d_us = [u, u1, u4, u5]\n",
    "all_3d_us = [u2, u3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<function u at 0x7fb701623f28>\n",
      "<function u1 at 0x7fb701623ea0>\n",
      "<function u4 at 0x7fb701623b70>\n",
      "<function u5 at 0x7fb701623ae8>\n"
     ]
    }
   ],
   "source": [
    "for a in all_2d_us:\n",
    "    print(sym.latex(a))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def prm(symu, f):\n",
    "    print(\n",
    "        \"Exact solution expression     = \"+ str(symu).replace('**', '^')+\"\\n\" +\n",
    "        \"Forcing term expression       = \"+ str(f).replace('**', '^') + \"\\n\" +\n",
    "        \"Boundary condition expression = \"+ str(symu).replace('**', '^')+ \"\\n\" +\n",
    "        \"Problem constants             = \"+ \"pi:\"+str(pi)\n",
    "    )\n",
    "\n",
    "def generate_prm(symu):\n",
    "    f = -symu.diff(x,2)-symu.diff(y,2)\n",
    "    f = sym.simplify(f)\n",
    "    prm(symu, f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Exact solution expression     = x*y*(x - 1)*(y - 1)*exp(-100*(x - 0.5)^2 - 100*(y - 0.117)^2)\n",
      "Forcing term expression       = (x*(x - 1)*(-y*(y - 1)*(200*y - 23.4)^2 + 200*y*(y - 1) + 2*y*(200*y - 23.4) + 2*(y - 1)*(200*y - 23.4) - 2) + y*(y - 1)*(-x*(x - 1)*(200*x - 100.0)^2 + 200*x*(x - 1) + 2*x*(200*x - 100.0) + 2*(x - 1)*(200*x - 100.0) - 2))*exp(-100*(x - 0.5)^2 - 100*(y - 0.117)^2)\n",
      "Boundary condition expression = x*y*(x - 1)*(y - 1)*exp(-100*(x - 0.5)^2 - 100*(y - 0.117)^2)\n",
      "Problem constants             = pi:3.141592653589793\n"
     ]
    }
   ],
   "source": [
    "from sympy import exp\n",
    "x,y,z = sym.var('x,y,z')\n",
    "\n",
    "generate_prm(u(x,y))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Exact solution expression     = (x^2 + y^2)^0.333333333333333*sin(2*atan2(x, y)/3)\n",
      "Forcing term expression       = 0\n",
      "Boundary condition expression = (x^2 + y^2)^0.333333333333333*sin(2*atan2(x, y)/3)\n",
      "Problem constants             = pi:3.141592653589793\n"
     ]
    }
   ],
   "source": [
    "generate_prm(u1(x,y))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Exact solution expression     = x*y*z*(x - 1)*(y - 1)*(z - 1)*exp(-100*(x - 0.5)^2 - 100*(y - 0.117)^2 - 100*(z - 0.331)^2)\n",
      "Forcing term expression       = z*(z - 1)*(x*(x - 1)*(-y*(y - 1)*(200*y - 23.4)^2 + 200*y*(y - 1) + 2*y*(200*y - 23.4) + 2*(y - 1)*(200*y - 23.4) - 2) + y*(y - 1)*(-x*(x - 1)*(200*x - 100.0)^2 + 200*x*(x - 1) + 2*x*(200*x - 100.0) + 2*(x - 1)*(200*x - 100.0) - 2))*exp(-100*(x - 0.5)^2 - 100*(y - 0.117)^2 - 100*(z - 0.331)^2)\n",
      "Boundary condition expression = x*y*z*(x - 1)*(y - 1)*(z - 1)*exp(-100*(x - 0.5)^2 - 100*(y - 0.117)^2 - 100*(z - 0.331)^2)\n",
      "Problem constants             = pi:3.141592653589793\n"
     ]
    }
   ],
   "source": [
    "generate_prm(u2(x,y,z))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Exact solution expression     = (x^2 + y^2 + z^2)^0.25\n",
      "Forcing term expression       = (0.75*(x^2 + y^2)*(x^2 + y^2 + z^2)^0.75 - 1.0*(x^2 + y^2 + z^2)^1.75)*(x^2 + y^2 + z^2)^(-2.5)\n",
      "Boundary condition expression = (x^2 + y^2 + z^2)^0.25\n",
      "Problem constants             = pi:3.141592653589793\n"
     ]
    }
   ],
   "source": [
    "generate_prm(u3(x,y,z))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Exact solution expression     = sin(2*pi*x)*sin(2*pi*y)\n",
      "Forcing term expression       = 8*pi^2*sin(2*pi*x)*sin(2*pi*y)\n",
      "Boundary condition expression = sin(2*pi*x)*sin(2*pi*y)\n",
      "Problem constants             = pi:3.141592653589793\n"
     ]
    }
   ],
   "source": [
    "generate_prm(u4(x,y))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Exact solution expression     = x*y*(x - 1)*(y - 1)*(10*exp(-100*(x - 0.5)^2 - 100*(y - 0.1)^2) + 0.1*exp(-100*(x - 0.1)^2 - 100*(y - 0.5)^2))\n",
      "Forcing term expression       = -x*(x - 1)*(y*(y - 1)*(0.1*(200*y - 100.0)^2*exp(-100*(x - 0.1)^2 - 100*(y - 0.5)^2) + 10*(200*y - 20.0)^2*exp(-100*(x - 0.5)^2 - 100*(y - 0.1)^2) - 2000*exp(-100*(x - 0.5)^2 - 100*(y - 0.1)^2) - 20.0*exp(-100*(x - 0.1)^2 - 100*(y - 0.5)^2)) - 2*y*((20.0*y - 10.0)*exp(-100*(x - 0.1)^2 - 100*(y - 0.5)^2) + (2000*y - 200.0)*exp(-100*(x - 0.5)^2 - 100*(y - 0.1)^2)) - 2*(y - 1)*((20.0*y - 10.0)*exp(-100*(x - 0.1)^2 - 100*(y - 0.5)^2) + (2000*y - 200.0)*exp(-100*(x - 0.5)^2 - 100*(y - 0.1)^2)) + 20*exp(-100*(x - 0.5)^2 - 100*(y - 0.1)^2) + 0.2*exp(-100*(x - 0.1)^2 - 100*(y - 0.5)^2)) - y*(y - 1)*(x*(x - 1)*(10*(200*x - 100.0)^2*exp(-100*(x - 0.5)^2 - 100*(y - 0.1)^2) + 0.1*(200*x - 20.0)^2*exp(-100*(x - 0.1)^2 - 100*(y - 0.5)^2) - 2000*exp(-100*(x - 0.5)^2 - 100*(y - 0.1)^2) - 20.0*exp(-100*(x - 0.1)^2 - 100*(y - 0.5)^2)) - 2*x*((20.0*x - 2.0)*exp(-100*(x - 0.1)^2 - 100*(y - 0.5)^2) + (2000*x - 1000.0)*exp(-100*(x - 0.5)^2 - 100*(y - 0.1)^2)) - 2*(x - 1)*((20.0*x - 2.0)*exp(-100*(x - 0.1)^2 - 100*(y - 0.5)^2) + (2000*x - 1000.0)*exp(-100*(x - 0.5)^2 - 100*(y - 0.1)^2)) + 20*exp(-100*(x - 0.5)^2 - 100*(y - 0.1)^2) + 0.2*exp(-100*(x - 0.1)^2 - 100*(y - 0.5)^2))\n",
      "Boundary condition expression = x*y*(x - 1)*(y - 1)*(10*exp(-100*(x - 0.5)^2 - 100*(y - 0.1)^2) + 0.1*exp(-100*(x - 0.1)^2 - 100*(y - 0.5)^2))\n",
      "Problem constants             = pi:3.141592653589793\n"
     ]
    }
   ],
   "source": [
    "generate_prm(u5(x,y))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "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.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}