examples/example-3-sindypde-reactiondiffusion.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Demo: Reaction Diffusion Equation (PDE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### created by Yuying Liu, 11/30/18"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import sys\n",
    "module_path = os.path.abspath(os.path.join('..'))\n",
    "if module_path not in sys.path:\n",
    "    sys.path.append(module_path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pySINDy.sindypde import SINDyPDE\n",
    "import scipy.io as sio\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# this .mat file can be generated from two of our .m files in datasets directory,\n",
    "# but since it's too large, we'll leave the user to generate the .mat file by themselves\n",
    "data = sio.loadmat('../datasets/reaction_diffusion.mat')\n",
    "data.keys()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "U = np.real(data['u'])\n",
    "V = np.real(data['v'])\n",
    "t = np.real(data['t'].flatten())\n",
    "x = np.real(data['x'].flatten())\n",
    "y = np.real(data['y'].flatten())\n",
    "dt = t[1] - t[0]\n",
    "dx = x[1] - x[0]\n",
    "dy = y[1] - y[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = SINDyPDE(name='SINDyPDE model for Reaction-Diffusion Eqn')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "U1 = U[100:200, 100:200, 200:230]\n",
    "V1 = V[100:200, 100:200, 200:230]\n",
    "model.fit({'u': U1, 'v': V1}, dt, [dx, dy], space_deriv_order=2, poly_degree=2, sample_rate=0.01, cut_off=0.05, deriv_acc=5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "activated1 = [model.descriptions[i] for i in np.arange(model.coefficients.shape[0]) if model.coefficients[i, 0] != 0]\n",
    "activated2 = [model.descriptions[i] for i in np.arange(model.coefficients.shape[0]) if model.coefficients[i, 1] != 0]\n",
    "\n",
    "print(activated1)\n",
    "print(activated2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model.coefficients"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model.descriptions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from findiff import FinDiff\n",
    "deriv_acc = 5\n",
    "\n",
    "U1 = U[100:200, 100:200, 200:230]\n",
    "V1 = V[100:200, 100:200, 200:230]\n",
    "\n",
    "d1_dt = FinDiff(U1.ndim-1, dt, 1, acc=deriv_acc)\n",
    "d2_xx = FinDiff(0, dx, 2, acc=deriv_acc)\n",
    "d2_yy = FinDiff(1, dy, 2, acc=deriv_acc)\n",
    "\n",
    "u_t = d1_dt(U1).flatten()\n",
    "v_t = d1_dt(V1).flatten()\n",
    "x_t = np.vstack([u_t, v_t]).T\n",
    "print('finished time derivative computation!')\n",
    "\n",
    "u_xx = d2_xx(U1).flatten()\n",
    "u_yy = d2_yy(U1).flatten()\n",
    "v_xx = d2_xx(V1).flatten()\n",
    "v_yy = d2_yy(V1).flatten()\n",
    "u = U1.flatten()\n",
    "v = V1.flatten()\n",
    "uv2 = (U1*V1*V1).flatten()\n",
    "u2v = (U1*U1*V1).flatten()\n",
    "u3 = (U1*U1*U1).flatten()\n",
    "v3 = (V1*V1*V1).flatten()\n",
    "\n",
    "lib = np.vstack([u_xx, u_yy, v_xx, v_yy, u, v, uv2, u2v, u3, v3]).T\n",
    "\n",
    "print(np.linalg.lstsq(lib, x_t, rcond=None)[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}