{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# AM783 Applied Markov Processes | From ODEs to SDEs"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Hugo Touchette\n",
"\n",
"Last updated: 10 October 2020\n",
"\n",
"Python 3"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"from scipy.integrate import odeint"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"# Magic command for vectorised figures\n",
"%config InlineBackend.figure_format = 'svg'"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Solving ODEs in Python"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let us solve Newton's equations for the pendulum:\n",
"$$\n",
"m\\ddot\\theta +\\frac{mg}{\\ell}\\sin\\theta = 0.\n",
"$$\n",
"As seen in class, we can transform this 2nd-order ODE into a system of two coupled 1st-order ODEs:\n",
"\\begin{eqnarray*}\n",
"\\dot\\theta &=& \\varphi \\\\\n",
"\\dot\\varphi &=&\\ddot \\theta = -\\omega^2\\sin\\theta,\n",
"\\end{eqnarray*}\n",
"where\n",
"$$\n",
"\\omega = \\sqrt{\\frac{g}{\\ell}}\n",
"$$\n",
"is the natural frequency of the pendulum.\n",
"\n",
"The first task, for simulating this system in Python, is to define as a Python function the force appearing on the right-hand side of the coupled ODEs. We have two such forces, so the function will return a vector (array):"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"def f(x, t):\n",
" # x is our vector; we choose x[0]=theta and x[1]=phi\n",
" g = 9.8 # Gravitational constant\n",
" l = 1.0 # Pendulum length\n",
" w = np.sqrt(g/l)\n",
"\n",
" dx0dt = x[1]\n",
" dx1dt = -w**2 * np.sin(x[0])\n",
"\n",
" return([dx0dt, dx1dt])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Calling odeint, we can then do the simulation:"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"tfinal = 10.0\n",
"dt = 0.01\n",
"n = int(tfinal/dt)\n",
"tspan = np.linspace(0, tfinal, n)\n",
"\n",
"x0 = [np.pi/5, 0] # Initial angle of pi/5; initial 0 velocity\n",
"\n",
"x = odeint(f, x0, tspan) # Solution\n",
"\n",
"plt.plot(tspan, x[:,0], label='Angle')\n",
"plt.plot(tspan, x[:,1], label='Angular velocity')\n",
"plt.xlabel(r'$t$')\n",
"plt.legend(loc='upper right')\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Solving SDEs with the Euler-Maruyama scheme"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We solve now the linear ODE\n",
"$$\n",
"\\dot x(t) = -\\gamma x(t)\n",
"$$\n",
"by adding Gaussian white noise. This means that the SDE to solve is\n",
"$$\n",
"dX_t = -\\gamma X_t dt+\\sigma dW_t.\n",
"$$\n",
"Here $\\gamma>0$ is the linear restitution constant, while $\\sigma>0$ is the noise amplitude.\n",
"\n",
"The Euler-Maruyama scheme for this SDE is direct:"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"tfinal = 10.0\n",
"dt = 0.01\n",
"sqrtdt = np.sqrt(dt)\n",
"gamma = 1.0\n",
"sigma = 0.1\n",
"\n",
"n = int(tfinal/dt)\n",
"x = np.zeros(n+1)\n",
"\n",
"x[0] = 3.0 # Initial value\n",
"\n",
"for i in range(n):\n",
" x[i+1] = x[i] -gamma*x[i]*dt +sigma*sqrtdt*np.random.normal()\n",
"\n",
"tspan = np.linspace(0, tfinal, n+1)\n",
"plt.plot(tspan, x)\n",
"plt.xlabel(r'$t$')\n",
"plt.ylabel(r'$X_t$')\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Compare with the deterministic trajectory obtained with the Euler scheme (one could also use odeint):"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"tfinal = 20.0\n",
"dt = 0.01\n",
"sqrtdt = np.sqrt(dt)\n",
"gamma = 1.0\n",
"sigma = 0.5\n",
"\n",
"n = int(tfinal/dt)\n",
"x = np.zeros(n+1)\n",
"xdet = np.zeros(n+1)\n",
"\n",
"x[0] = 3.0 # Initial value\n",
"xdet[0] = x[0]\n",
"\n",
"for i in range(n):\n",
" x[i+1] = x[i] -gamma*x[i]*dt + sigma*sqrtdt*np.random.normal()\n",
" xdet[i+1] = xdet[i] -gamma*xdet[i]*dt\n",
"\n",
"tspan = np.linspace(0, tfinal, n+1)\n",
"plt.plot(tspan, x)\n",
"plt.plot(tspan, xdet, 'r-')\n",
"plt.xlabel(r'$t$')\n",
"plt.ylabel(r'$X_t$')\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Change the noise parameter $\\sigma$ to see the effect of the noise."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"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.9.12"
},
"varInspector": {
"cols": {
"lenName": 16,
"lenType": 16,
"lenVar": 40
},
"kernels_config": {
"python": {
"delete_cmd_postfix": "",
"delete_cmd_prefix": "del ",
"library": "var_list.py",
"varRefreshCmd": "print(var_dic_list())"
},
"r": {
"delete_cmd_postfix": ") ",
"delete_cmd_prefix": "rm(",
"library": "var_list.r",
"varRefreshCmd": "cat(var_dic_list()) "
}
},
"types_to_exclude": [
"module",
"function",
"builtin_function_or_method",
"instance",
"_Feature"
],
"window_display": false
}
},
"nbformat": 4,
"nbformat_minor": 2
}