foo.ipynb

foo.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 261,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$w_x(t_k) = \\frac{1}{h^2}\\cdot (x_{k+1} - 2\\cdot x_k + x_{k - 1}) - \\frac{dV}{dx}(x_k, y_k)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 262,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def Wx(dt, x, y):\n",
    "    return (x[2:] - 2 * x[1:-1] + x[:-2]) / (dt ** 2) - gvx(x[1:-1], y[1:-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$w_y(t_k) = \\frac{1}{h^2}\\cdot (y_{k+1} - 2\\cdot y_k + y_{k-1}) - \\frac{dV}{dy}(x_k, y_k)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 263,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def Wy(dt, x, y):\n",
    "    return (y[2:] - 2 * x[1:-1] + x[:-2]) / (dt ** 2) - gvy(x[1:-1], y[1:-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$W(t) = W_x^2(t) + W_y^2(t)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 264,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def W(h, x, y):\n",
    "    return Wx(h, x, y) ** 2 + Wy(h, x, y) ** 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\tau = \\sum_{k=1}^{n-1}2hW(t_k)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 265,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def tau(h, x, y, w):\n",
    "    return np.sum(2 * h * w(h, x[1:-1], y[1:-1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\frac{d\\tau}{dx_k} = 4h\\cdot \\bigg(W_x(t_k)\\cdot (-\\frac{2}{h}-\\frac{d^2V}{dx^2}(x_k,y_k)) + \\frac{W_x(t_{k-1}) + W_x(t_{k+1})}{h^2} + W_y(t_k)\\cdot \\frac{dV}{dxdy}(x_k, y_k)\\bigg)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 266,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def gtaux(h, x, y, wx, wy, d2vx, dvxy):\n",
    "    # without shift\n",
    "    wos = wx[1:-1] * (-2/(h ** 2) - d2vx(x[1:-1], y[1:-1]))\n",
    "    # with shift\n",
    "    ws = (wx[:-2] + wx[2:]) / (h ** 2)\n",
    "    # wy part\n",
    "    yp = wy[1:-1] * dvxy(x[1:-1], y[1:-1])\n",
    "\n",
    "    return 4 * h * (wos + ws + yp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\frac{d\\tau}{dy_k} = 4h\\cdot \\bigg(W_y(t_k)\\cdot (-\\frac{2}{h^2} - \\frac{d^2V}{dy^2}(x_k, y_k)) + \\frac{W_y(t_{k-1}) + W_y(t_{k + 1})}{h^2} + W_x(t_k)\\cdot \\frac{d^2V}{dydx}(x_k,y_k)\\bigg)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 267,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def gtauy(h, x, y, wx, wy, d2vy, dvxy):\n",
    "    # without shift\n",
    "    wos = wy[1:-1] * (-2/(h ** 2) - d2vy(x[1:-1], y[1:-1]))\n",
    "    # with shift\n",
    "    ws = (wy[:-2] + wy[2:]) / (h ** 2)\n",
    "    # wx part\n",
    "    xp = wx[1:-1] * dvxy(x[1:-1], y[1:-1])\n",
    "\n",
    "    return 4 * h * (wos + ws + xp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$V, \\frac{dV}{dx}, \\frac{dV}{dy}, \\frac{d^2V}{dx^2}, \\frac{d^2V}{dy^2}, \\frac{d^2V}{dxdy}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 268,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "###  V\n",
    "def v(x, y):\n",
    "    return np.sin(x) * np.sin(y)\n",
    "\n",
    "# \\frac{dv}{dx}\n",
    "def dvdx(x, y):\n",
    "    return np.cos(x) * np.sin(y)\n",
    "\n",
    "# \\frac{dv}{dy}\n",
    "def dvdy(x, y):\n",
    "    return np.sin(x) * np.cos(y)\n",
    "\n",
    "# \\frac{d^2 v}{dx^2}\n",
    "def d2vdx(x, y):\n",
    "    return -np.sin(x) * np.sin(y)\n",
    "\n",
    "# \\frac{d^2 v}{dy^2}\n",
    "def d2vdy(x, y):\n",
    "    return -np.sin(y) * np.sin(x)\n",
    "\n",
    "# \\frac{d^2 v}{dx \\cdot dy}\n",
    "def dvdxy(x, y):\n",
    "    return np.cos(x) * np.cos(y)\n",
    "\n",
    "### /V\n",
    "\n",
    "def gvx(x, y): return dvdx(x, y)\n",
    "def gvy(x, y): return dvdy(x, y)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 269,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "EPS = 1e-8\n",
    "DT = 1\n",
    "LEN = 5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 270,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "X0 = Y0 = 3 * np.pi / 2\n",
    "X1 = Y1 = 5 * np.pi / 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 271,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "x0 = np.linspace(X0, X1, LEN)\n",
    "y0 = np.linspace(Y0, Y1, LEN)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 272,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dx = [  0.00000000e+00   1.78144226e-08   1.78144226e-08   1.78144226e-08\n",
      "   0.00000000e+00]\n",
      "dy = [  0.00000000e+00   9.12327299e-09   9.12327299e-09   9.12327299e-09\n",
      "   0.00000000e+00]\n"
     ]
    }
   ],
   "source": [
    "dx = np.zeros(x0.shape)\n",
    "dy = np.zeros(y0.shape)\n",
    "\n",
    "dx[1:-1] = EPS * np.random.randn()\n",
    "dy[1:-1] = EPS * np.random.randn()\n",
    "\n",
    "print('dx =', dx)\n",
    "print('dy =', dy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 273,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "x1 = x0 + dx\n",
    "y1 = y0 + dy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 274,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tau0 = 2.399615652258972e-31, tau1 = 1.5715589623682999e-15\n"
     ]
    }
   ],
   "source": [
    "tau0 = tau(DT, x0, y0, W)\n",
    "tau1 = tau(DT, x1, y1, W)\n",
    "print(f'tau0 = {tau0}, tau1 = {tau1}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 275,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "gtau1x = [  0.00000000e+00  -1.06022293e-07  -1.40786892e-07  -1.06022291e-07\n",
      "   0.00000000e+00]\n",
      "gtau1y = [  0.00000000e+00  -7.12576913e-08  -3.64930922e-08  -7.12576932e-08\n",
      "   0.00000000e+00]\n"
     ]
    }
   ],
   "source": [
    "wx1 = Wx(DT, x1, y1)\n",
    "wy1 = Wy(DT, x1, y1)\n",
    "gtau1x = np.zeros(LEN)\n",
    "gtau1y = np.zeros(LEN)\n",
    "gtau1x[1:-1] = gtaux(DT, x1, y1, wx1, wy1, d2vdx, dvdxy)\n",
    "gtau1y [1:-1]= gtauy(DT, x1, y1, wx1, wy1, d2vdy, dvdxy)\n",
    "print('gtau1x =', gtau1x)\n",
    "print('gtau1y =', gtau1y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 276,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-7.91863222028e-15\n"
     ]
    }
   ],
   "source": [
    "foo = np.dot(gtau1x, dx) + np.dot(gtau1y, dy)\n",
    "print(foo)"
   ]
  }
 ],
 "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}