Merge pull request #415 from DiODeProject/new-version-variance-notebook

Updated the demo notebook Variance Suppression with new analysis of the stochastic noise.
DiODeProject · Mar 18, 2022 · 470af4a · 470af4a
2 parents 223fbf6 + 58be0b6
commit 470af4a
Showing 1 changed file with 159 additions and 0 deletions.
diff --git a/DemoNotebooks/Variance_suppression.ipynb b/DemoNotebooks/Variance_suppression.ipynb
@@ -259,6 +259,165 @@
     "                          choose_yrange=[0,1], plotProportions=True, ylab=\"Subpopulations $x_i$\")"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Fokker-Planck equations and noise variance\n",
+    "\n",
+    "MuMoT allows you to easily compute the Fokker-Planck equations for both models and display the first-order and second-order moments of the system's noise."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "posmodel.showFokkerPlanckEquation()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "negmodel.showFokkerPlanckEquation()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Noise equations\n",
+    "\n",
+    "Having the Fokker-Planck equations allows us to derive the equations of motion for the fluctuations.\n",
+    "\n",
+    "### Model without negative social feedback"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "posmodel.showNoiseEquations()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Because MuMoT allows you to access the `sympy` equations, we apply some simplification to the equations in order to find the solution to $\\langle \\eta_{A}^2 \\rangle$ at convergence ($t\\rightarrow \\infty$). We assume that $a\\approx0$ and $\\Phi_U (t\\rightarrow \\infty) \\approx 0$."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "noiseEqPos=posmodel.getNoiseEquations()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#in order to handle sympy equations, some imports are necessary\n",
+    "from sympy import (\n",
+    "    collect,\n",
+    "    default_sort_key,\n",
+    "    Derivative,\n",
+    "    lambdify,\n",
+    "    latex,\n",
+    "    linsolve,\n",
+    "    numbered_symbols,\n",
+    "    preview,\n",
+    "    simplify,\n",
+    "    solve,\n",
+    "    Symbol,\n",
+    "    symbols,\n",
+    "    Function\n",
+    ")\n",
+    "from IPython.display import display, Math"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "simplifiedNoiseEqPos={}\n",
+    "\n",
+    "simplifiedNoiseEqPos['etaU2']=noiseEqPos[2][Function('M_2')(Symbol('eta_U')**2)].subs(Symbol('Phi_U'),0).subs(Symbol('a'),0)\n",
+    "display(Math(latex(simplifiedNoiseEqPos['etaU2'])+'=0'))\n",
+    "simplifiedNoiseEqPos['etaU2']=solve( simplifiedNoiseEqPos['etaU2'], Function('M_2')(Symbol('eta_U')**2) )[0]\n",
+    "display(Math(latex( Function('M_2')(Symbol('eta_U')**2) ) + '=' + latex(simplifiedNoiseEqPos['etaU2']) ))\n",
+    "\n",
+    "simplifiedNoiseEqPos['etaA*etaU']=noiseEqPos[2][Function('M_2')(Symbol('eta_A')*Symbol('eta_U'))].subs(\n",
+    "    Function('M_2')(Symbol('eta_U')**2),simplifiedNoiseEqPos['etaU2']).subs(Symbol('Phi_U'),0).subs(Symbol('a'),0)\n",
+    "display(Math(latex(simplifiedNoiseEqPos['etaA*etaU'])+'=0'))\n",
+    "simplifiedNoiseEqPos['etaA*etaU']=solve( simplifiedNoiseEqPos['etaA*etaU'], Function('M_2')(Symbol('eta_A')*Symbol('eta_U')) )[0]\n",
+    "display(Math(latex( Function('M_2')(Symbol('eta_A')*Symbol('eta_U')) ) + '=' + latex(simplifiedNoiseEqPos['etaA*etaU']) ))\n",
+    "\n",
+    "simplifiedNoiseEqPos['etaA2']=noiseEqPos[2][Function('M_2')(Symbol('eta_A')**2)].subs(\n",
+    "    Function('M_2')(Symbol('eta_A')*Symbol('eta_U')),simplifiedNoiseEqPos['etaA*etaU']).subs(Symbol('Phi_U'),0)\n",
+    "display(Math(latex(simplifiedNoiseEqPos['etaA2'])+'=0'))\n",
+    "simplifiedNoiseEqPos['etaA2']=solve( simplifiedNoiseEqPos['etaA2'], Function('M_2')(Symbol('eta_A')**2) )[0]\n",
+    "display(Math(latex( Function('M_2')(Symbol('eta_A')**2) ) + '=' + latex(simplifiedNoiseEqPos['etaA2']) ))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Model with negative social feedback"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "negmodel.showNoiseEquations()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "noiseEqNeg=negmodel.getNoiseEquations()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "For the system with negative feedback, we assume that $\\Phi_U (t\\rightarrow \\infty) \\approx 0$ and $\\langle \\eta_{A}\\eta_{U} \\rangle \\approx 0$ to find the solution to $\\langle \\eta_{A}^2 \\rangle$:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "simplifiedNoiseEqPos=noiseEqPos[2][Function('M_2')(Symbol('eta_A')**2)].subs(Symbol('Phi_U'),0).subs(\n",
+    "    Function('M_2')(Symbol('eta_A')*Symbol('eta_U')),0)\n",
+    "display(Math(latex(simplifiedNoiseEqPos)+'=0'))\n",
+    "simplifiedNoiseSizePos=solve( simplifiedNoiseEqPos, Function('M_2')(Symbol('eta_A')**2) )[0]\n",
+    "display(Math(latex( Function('M_2')(Symbol('eta_A')**2) ) + '=' + latex(simplifiedNoiseSizePos) ))"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},