![](\"images/cambridge.png\")
![](\"images/caplogo.svg\")
![](\"images/kavli.jpg\")
\n",
"\n",
"\n", "
\n", "
\n", "
\n", "
\n", "\n", "To run these slides in jupyter,\\\n", "head to https://github.com/fruzsinaagocs/oscode and\n", "\n", "
![](\"https://mybinder.org/badge_logo.svg\"){kind=icon}
\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Introduction"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"- pyoscode is a numerical solver for **osc**illatory **o**rdinary **d**ifferential **e**quations\n",
"- ODEs can be oscillatory through\n",
" - an oscillatory forcing term, $$\\ddot{x} + x = C\\cos{\\omega t},\\; \\omega \\gg 1$$\n",
" - oscillatory coefficients, $$\\dot{x} = xt + 10x\\cos{\\omega t},\\; \\omega \\gg 1$$\n",
" - slowly changing, but high frequency, $$\\ddot{x} + 2\\gamma(t)\\dot{x} + \\omega^2(t) x = 0,\\; \\omega \\gg 1$$\n",
"- oscode deals with the third type:\n",
" - one-dimensional\n",
" - ordinary\n",
" - $\\omega(t)$, $\\gamma(t)$ vary slowly for _at least part of the integration range_"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Motivation"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"- all extremely common in physics: pendula, suspension, circuits, celestial mechanics, ...\n",
"- even if $\\omega(t)$ changes slowly, if $\\omega \\gg 1$, numerical solution slow\n",
"- e.g. $\\ddot{x} + tx = 0$, the Airy equation:"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"slideshow": {
"slide_type": "-"
}
},
"outputs": [],
"source": [
"def f(t,x):\n",
" return np.array([x[1],-t*x[0]])\n",
"\n",
"ti = 1.0\n",
"tf = 1e3\n",
"x0 = np.array([sp.airy(-ti)[0], -sp.airy(-ti)[1]])"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"slideshow": {
"slide_type": "-"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"True\n",
"x(t=1000.0)=0.0302, dx/dt=-2.3852\n",
"analytic solution: x=0.0560, dx/dt=-2.6331\n",
"CPU times: user 1.63 s, sys: 0 ns, total: 1.63 s\n",
"Wall time: 1.64 s\n"
]
}
],
"source": [
"%%time\n",
"sol = scipy.integrate.solve_ivp(f,(ti,tf),x0,rtol=1e-3)\n",
"print(sol.success)\n",
"print(\"x(t={})={:.4f}, dx/dt={:.4f}\".format(tf,sol.y[0][-1],sol.y[1][-1]))\n",
"print(\"analytic solution: x={:.4f}, dx/dt={:.4f}\".format(sp.airy(-tf)[0],-sp.airy(-tf)[1]))"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"- why the noticeable wall time?\n",
" - `scipy.integrate.solve_ivp`'s default method is 4(5) order Runge-Kutta\n",
" - Runge-Kutta excellent for non-stiff, non-oscillatory equations/regions\n",
" - cannot deal with high-frequency oscillations (will see why)\n",
" - the Airy frequency $\\omega(t) = \\sqrt{t}$, steadily increases with time"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## A brief summary of the algorithm\n",
"\n",
"Based on [arxiv:1906.01421](https://arxiv.org/abs/1906.01421)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Runge-Kutta (RK) steps \n",
"\n",
"- $ \\dot{\\vec{x}} = F(\\vec{x}) $\n",
"- have (estimate of) solution $x$ at $t_i$, want to get $x$ at $t_i+h$\n",
"- Taylor expand: \n",
"\n",
"$$x(t_i+h) = x(t_i) + hF(x(t_i)) + \\frac{1}{2}h^2\\left.\\frac{dF}{dt}\\right\\vert_{t=t_i} + \\mathcal{O}(h^3)$$\n",
"\n",
"- and replace derivatives of $F$ with evaluations of $F$ at gridpoints:\n",
"\n",
"$$ x(t_i+h) = x(t_i) + \\sum_{j,\\;t_i![](\"images/cmb.png\")
\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"- power spectrum: $|$amplitude$|^2$ as a function of Fourier wavenumber $k$\n",
"- computation of spectrum _can be_ bottleneck in forward modelling, e.g.\n",
" - closed universes [arxiv:1907.08524](https://arxiv.org/abs/1907.08524)\n",
" - axion monodromy models [arxiv:1502.02114](https://arxiv.org/abs/1502.02114)\n",
" - initial conditions set at very early times [arXiv:2002.07042](https://arxiv.org/abs/2002.07042)\n",
"\n",
"