{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Test for different SME solvers against analytical solution for oscillator squeezing\n", "\n", "Manuel Grimm, Niels Lörch, and Denis V. Vasilyev\n", "\n", "30 August 2016\n", "\n", "Updated by Eric Giguere\n", "March 2018\n", "\n", "We solve the stochastic master equation for an oscillator coupled to a 1D field as discussed in [1]. There is a deterministic differential equation for the variances of the oscillator quadratures $\\langle\\delta X^2\\rangle$ and $\\langle\\delta P^2\\rangle$. This allows for a direct comparison between the numerical solution and the exact solution for a single quantum trajectory. In particular, we study scaling of deviations from analytical solution as a function of stepsize for different solvers:\n", "'euler-maruyama', 'pc-euler', 'milstein', 'milstein-imp', 'taylor15', 'taylor15-imp'.\n", "\n", "It is important to check the correct scaling since it is very easy to implement a higher order method in a wrong way such that it still works but as a lower order method.\n", "\n", "In this section we solve SME with a single Wiener increment: \n", "### $\\mathrm{d}\\rho = D[s]\\rho\\mathrm{d}t + H[s]\\rho \\mathrm{d}W + \\gamma D[a]\\rho\\mathrm{d}t$\n", "\n", "The steady state solution for the variance $V_{\\mathrm{c}} = \\langle X^2\\rangle - \\langle X\\rangle^2$ reads\n", "\n", "$V_{\\mathrm{c}} = \\frac1{4\\alpha^{2}}\\left[\\alpha\\beta - \\gamma + \\sqrt{(\\gamma-\\alpha\\beta )^{2} + 4\\gamma \\alpha^2}\\right]$\n", "\n", "where $\\alpha$ and $\\beta$ are parametrizing the interaction between light and the oscillator such that the jump operator is given by $s = \\frac{\\alpha+\\beta}2 a + \\frac{\\alpha-\\beta}2 a^{\\dagger}$\n", "\n", "[1] D. V. Vasilyev, C. a. Muschik, and K. Hammerer, Physical Review A 87, 053820 (2013). arXiv:1303.5888" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "%config InlineBackend.figure_formats = ['svg']\n", "\n", "from qutip import *\n", "from qutip.ui.progressbar import BaseProgressBar\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from scipy.integrate import odeint\n", "y_sse = None\n", "import time" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Just check that analytical solution coincides with the solution of ODE for the variance" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n" ], "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "def arccoth(x):\n", " return 0.5*np.log((1.+x)/(x-1.))\n", "\n", "############ parameters #############\n", "\n", "th = 0.1 # Interaction parameter\n", "alpha = np.cos(th)\n", "beta = np.sin(th)\n", "gamma = 1.\n", "def gammaf(t):\n", " return 0.25+t/12+t*t/6\n", "\n", "def f_gamma(t,*args):\n", " return (0.25+t/12+t*t/6)**(0.5)\n", "\n", "################# Solution of the differential equation for the variance Vc ####################\n", "T = 6.\n", "N_store = int(20*T+1)\n", "tlist = np.linspace(0,T,N_store)\n", "y0 = 0.5\n", "def func(y, t):\n", " return -(gammaf(t) - alpha*beta)*y - 2*alpha*alpha*y*y + 0.5*gammaf(t)\n", "y_td = odeint(func, y0, tlist)\n", "y_td = y_td.ravel()\n", "\n", "def func(y, t):\n", " return -(gamma - alpha*beta)*y - 2*alpha*alpha*y*y + 0.5*gamma\n", "y = odeint(func, y0, tlist)\n", "\n", "############ Exact steady state solution for Vc #########################\n", "Vc = (alpha*beta - gamma + np.sqrt((gamma-alpha*beta)**2 + 4*gamma*alpha**2))/(4*alpha**2)\n", "\n", "#### Analytic solution\n", "A = (gamma**2 + alpha**2 * (beta**2 + 4*gamma) - 2*alpha*beta*gamma)**0.5\n", "B = arccoth((-4*alpha**2*y0 + alpha*beta - gamma)/A)\n", "y_an = (alpha*beta - gamma + A / np.tanh(0.5*A*tlist - B))/(4*alpha**2)\n", "\n", "f, (ax, ax2) = plt.subplots(2, 1, sharex=True)\n", "\n", "ax.set_title('Variance as a function of time')\n", "ax.plot(tlist,y)\n", "ax.plot(tlist,Vc*np.ones_like(tlist))\n", "ax.plot(tlist,y_an)\n", "ax.set_ylim(0,0.5)\n", "\n", "ax2.set_title('Deviation of odeint from analytic solution')\n", "ax2.set_xlabel('t')\n", "ax2.set_ylabel(r'$\\epsilon$')\n", "ax2.plot(tlist,y_an - y.T[0]);" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "scrolled": true }, "outputs": [], "source": [ "####################### Model ###########################\n", "\n", "N = 30 # number of Fock states\n", "Id = qeye(N)\n", "a = destroy(N)\n", "s = 0.5*((alpha+beta)*a + (alpha-beta)*a.dag())\n", "x = (a + a.dag())/np.sqrt(2)\n", "H = Id\n", "c_op = [np.sqrt(gamma)*a]\n", "c_op_td = [[a,f_gamma]]\n", "sc_op = [s]\n", "e_op = [x, x*x]\n", "rho0 = fock_dm(N,0) # initial vacuum state\n", "#sc_len=1 # one stochastic operator\n", "\n", "\n", "############## time steps and trajectories ###################\n", "\n", "ntraj = 1 #100 # number of trajectories\n", "T = 6. # final time \n", "N_store = int(20*T+1) # number of time steps for which we save the expectation values/density matrix\n", "tlist = np.linspace(0,T,N_store)\n", "ddt = (tlist[1]-tlist[0])\n", "\n", "Nsubs = (10*np.logspace(0,1,10)).astype(int)\n", "stepsizes = [ddt/j for j in Nsubs] # step size is doubled after each evaluation \n", "Nt = len(Nsubs) # number of step sizes that we compare\n", "Nsubmax = Nsubs[-1] # Number of intervals for the smallest step size; \n", "dtmin = (tlist[1]-tlist[0])/(Nsubmax)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Test of different SME solvers" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plotting the figure" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "{'solver': 'euler-maruyama'}\n", "Total run time: 0.05s\n", "Total run time: 0.07s\n", "Total run time: 0.11s\n", "Total run time: 0.11s\n", "Total run time: 0.16s\n", "Total run time: 0.19s\n", "Total run time: 0.27s\n", "Total run time: 0.36s\n", "Total run time: 0.42s\n", "Total run time: 0.55s\n", "{'solver': 'platen'}\n", "Total run time: 0.16s\n", "Total run time: 0.21s\n", "Total run time: 0.27s\n", "Total run time: 0.36s\n", "Total run time: 0.45s\n", "Total run time: 0.61s\n", "Total run time: 0.74s\n", "Total run time: 0.95s\n", "Total run time: 1.24s\n", "Total run time: 1.64s\n", "{'solver': 'pred-corr'}\n", "Total run time: 0.14s\n", "Total run time: 0.16s\n", "Total run time: 0.18s\n", "Total run time: 0.25s\n", "Total run time: 0.32s\n", "Total run time: 0.41s\n", "Total run time: 0.54s\n", "Total run time: 0.68s\n", "Total run time: 0.89s\n", "Total run time: 1.15s\n", "{'solver': 'milstein'}\n", "Total run time: 0.08s\n", "Total run time: 0.09s\n", "Total run time: 0.17s\n", "Total run time: 0.16s\n", "Total run time: 0.22s\n", "Total run time: 0.28s\n", "Total run time: 0.37s\n", "Total run time: 0.48s\n", "Total run time: 0.62s\n", "Total run time: 0.79s\n", "{'solver': 'milstein-imp'}\n", "Total run time: 0.60s\n", "Total run time: 0.70s\n", "Total run time: 0.91s\n", "Total run time: 1.12s\n", "Total run time: 1.43s\n", "Total run time: 1.88s\n", "Total run time: 2.46s\n", "Total run time: 3.15s\n", "Total run time: 4.09s\n", "Total run time: 5.30s\n", "{'solver': 'pred-corr-2'}\n", "Total run time: 0.17s\n", "Total run time: 0.23s\n", "Total run time: 0.29s\n", "Total run time: 0.37s\n", "Total run time: 0.51s\n", "Total run time: 0.62s\n", "Total run time: 0.81s\n", "Total run time: 1.03s\n", "Total run time: 1.34s\n", "Total run time: 1.78s\n", "{'solver': 'explicit1.5'}\n", "Total run time: 0.37s\n", "Total run time: 0.43s\n", "Total run time: 0.66s\n", "Total run time: 0.73s\n", "Total run time: 0.93s\n", "Total run time: 1.30s\n", "Total run time: 1.57s\n", "Total run time: 2.01s\n", "Total run time: 2.66s\n", "Total run time: 3.40s\n", "{'solver': 'taylor1.5'}\n", "Total run time: 0.29s\n", "Total run time: 0.34s\n", "Total run time: 0.47s\n", "Total run time: 0.59s\n", "Total run time: 0.75s\n", "Total run time: 0.97s\n", "Total run time: 1.28s\n", "Total run time: 1.66s\n", "Total run time: 2.20s\n", "Total run time: 2.79s\n", "{'solver': 'taylor1.5-imp', 'args': {'tol': 1e-08}}\n", "Total run time: 0.83s\n", "Total run time: 1.02s\n", "Total run time: 1.31s\n", "Total run time: 1.70s\n", "Total run time: 2.03s\n", "Total run time: 2.48s\n", "Total run time: 3.22s\n", "Total run time: 4.10s\n", "Total run time: 5.28s\n", "Total run time: 6.94s\n", "{'solver': 'taylor2.0'}\n", "Total run time: 0.40s\n", "Total run time: 0.48s\n", "Total run time: 0.61s\n", "Total run time: 0.83s\n", "Total run time: 1.03s\n", "Total run time: 1.32s\n", "Total run time: 1.75s\n", "Total run time: 2.22s\n", "Total run time: 2.93s\n", "Total run time: 3.74s\n", "{'solver': 'taylor2.0', 'noiseDepth': 500}\n", "Total run time: 1.98s\n", "Total run time: 2.34s\n", "Total run time: 3.13s\n", "Total run time: 4.09s\n", "Total run time: 5.33s\n", "Total run time: 6.81s\n", "Total run time: 8.93s\n", "Total run time: 11.43s\n", "Total run time: 15.03s\n", "Total run time: 19.46s\n" ] } ], "source": [ "ntraj = 1\n", "\n", "def run_cte_cte(**kwargs):\n", " epsilon = np.zeros(Nt)\n", " std = np.zeros(Nt)\n", " print(kwargs)\n", " for j in range(0,ntraj):\n", " for jj in range(0,Nt):\n", " Nsub = Nsubs[jj]\n", " sol = smesolve(H, rho0, tlist, c_op, sc_op, e_op, nsubsteps=Nsub, **kwargs)\n", " epsilon_j = 1/T * np.sum(np.abs(y_an - (sol.expect[1]-sol.expect[0]*sol.expect[0].conj())))*ddt\n", " epsilon[jj] += epsilon_j\n", " std[jj] += epsilon_j\n", " epsilon/= ntraj\n", " std = np.sqrt(1/ntraj * (1/ntraj * std - epsilon**2))\n", " return epsilon\n", "\n", "def get_stats_cte_cte(**kw):\n", " start = time.time()\n", " y = run_cte_cte(**kw)\n", " tag = str(kw[\"solver\"])\n", " x = np.log(stepsizes)\n", " ly = np.log(y)\n", " fit = np.polyfit(x, ly, 1)[0]\n", " return y,tag,fit,time.time()-start\n", "\n", "stats_cte_cte = []\n", "stats_cte_cte.append(get_stats_cte_cte(solver='euler-maruyama'))\n", "\n", "stats_cte_cte.append(get_stats_cte_cte(solver='platen'))\n", "stats_cte_cte.append(get_stats_cte_cte(solver='pred-corr'))\n", "stats_cte_cte.append(get_stats_cte_cte(solver='milstein'))\n", "stats_cte_cte.append(get_stats_cte_cte(solver='milstein-imp'))\n", "stats_cte_cte.append(get_stats_cte_cte(solver='pred-corr-2'))\n", "\n", "stats_cte_cte.append(get_stats_cte_cte(solver='explicit1.5'))\n", "stats_cte_cte.append(get_stats_cte_cte(solver=\"taylor1.5\"))\n", "stats_cte_cte.append(get_stats_cte_cte(solver=\"taylor1.5-imp\", args={\"tol\":1e-8}))\n", "\n", "stats_cte_cte.append(get_stats_cte_cte(solver=\"taylor2.0\"))\n", "stats_cte_cte.append(get_stats_cte_cte(solver=\"taylor2.0\", noiseDepth=500))" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n" ], "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "fig = plt.figure()\n", "ax = plt.subplot(111)\n", "\n", "mark = \"o*vspx+^<>1hdD\"\n", "\n", "for i,run in enumerate(stats_cte_cte):\n", " ax.loglog(stepsizes, run[0], mark[i], label=run[1]+\": \" + str(run[2]))\n", "\n", "ax.loglog(stepsizes, 0.1*np.array(stepsizes)**0.5, label=\"$\\propto\\Delta t^{1/2}$\")\n", "ax.loglog(stepsizes, 0.1*np.array(stepsizes)**1, label=\"$\\propto\\Delta t$\")\n", "ax.loglog(stepsizes, 0.1*np.array(stepsizes)**1.5, label=\"$\\propto\\Delta t^{3/2}$\")\n", "ax.loglog(stepsizes, 0.2*np.array(stepsizes)**2.0, label=\"$\\propto\\Delta t^{2}$\")\n", "\n", "lgd=ax.legend(loc='center left', bbox_to_anchor=(1, 0.64), prop={'size':12})" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Deterministic part depend on time" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "{'solver': 'euler-maruyama'}\n", "Total run time: 0.25s\n", "Total run time: 0.13s\n", "Total run time: 0.18s\n", "Total run time: 0.23s\n", "Total run time: 0.29s\n", "Total run time: 0.38s\n", "Total run time: 0.52s\n", "Total run time: 0.60s\n", "Total run time: 0.80s\n", "Total run time: 0.99s\n", "{'solver': 'platen'}\n", "Total run time: 0.28s\n", "Total run time: 0.32s\n", "Total run time: 0.42s\n", "Total run time: 0.54s\n", "Total run time: 0.69s\n", "Total run time: 0.88s\n", "Total run time: 1.20s\n", "Total run time: 1.50s\n", "Total run time: 1.93s\n", "Total run time: 2.55s\n", "{'solver': 'pc-euler'}\n", "Total run time: 0.16s\n", "Total run time: 0.21s\n", "Total run time: 0.28s\n", "Total run time: 0.37s\n", "Total run time: 0.45s\n", "Total run time: 0.57s\n", "Total run time: 0.76s\n", "Total run time: 1.00s\n", "Total run time: 1.24s\n", "Total run time: 1.62s\n", "{'solver': 'milstein'}\n", "Total run time: 0.13s\n", "Total run time: 0.17s\n", "Total run time: 0.21s\n", "Total run time: 0.28s\n", "Total run time: 0.35s\n", "Total run time: 0.44s\n", "Total run time: 0.59s\n", "Total run time: 0.76s\n", "Total run time: 0.96s\n", "Total run time: 1.25s\n", "{'solver': 'milstein-imp'}\n", "Total run time: 0.89s\n", "Total run time: 1.07s\n", "Total run time: 1.36s\n", "Total run time: 1.79s\n", "Total run time: 2.20s\n", "Total run time: 2.77s\n", "Total run time: 3.65s\n", "Total run time: 4.65s\n", "Total run time: 6.09s\n", "Total run time: 7.90s\n", "{'solver': 'pc-euler-2'}\n", "Total run time: 0.28s\n", "Total run time: 0.33s\n", "Total run time: 0.46s\n", "Total run time: 0.57s\n", "Total run time: 0.73s\n", "Total run time: 0.93s\n", "Total run time: 1.23s\n", "Total run time: 1.61s\n", "Total run time: 2.04s\n", "Total run time: 2.71s\n", "{'solver': 'explicit1.5'}\n", "Total run time: 0.50s\n", "Total run time: 0.59s\n", "Total run time: 0.77s\n", "Total run time: 1.02s\n", "Total run time: 1.32s\n", "Total run time: 1.71s\n", "Total run time: 2.19s\n", "Total run time: 2.87s\n", "Total run time: 3.72s\n", "Total run time: 4.80s\n", "{'solver': 'taylor1.5'}\n", "Total run time: 0.50s\n", "Total run time: 0.56s\n", "Total run time: 0.86s\n", "Total run time: 0.98s\n", "Total run time: 1.29s\n", "Total run time: 1.64s\n", "Total run time: 2.12s\n", "Total run time: 2.72s\n", "Total run time: 3.52s\n", "Total run time: 4.57s\n", "{'solver': 'taylor1.5-imp', 'args': {'tol': 1e-08}}\n", "Total run time: 1.17s\n", "Total run time: 1.39s\n", "Total run time: 1.85s\n", "Total run time: 2.31s\n", "Total run time: 2.94s\n", "Total run time: 3.71s\n", "Total run time: 4.79s\n", "Total run time: 6.00s\n", "Total run time: 7.68s\n", "Total run time: 9.94s\n", "{'solver': 'taylor2.0'}\n", "Total run time: 0.67s\n", "Total run time: 0.74s\n", "Total run time: 0.97s\n", "Total run time: 1.28s\n", "Total run time: 1.66s\n", "Total run time: 2.15s\n", "Total run time: 2.78s\n", "Total run time: 3.56s\n", "Total run time: 4.66s\n", "Total run time: 6.11s\n", "{'solver': 'taylor2.0', 'noiseDepth': 500}\n", "Total run time: 2.18s\n", "Total run time: 2.63s\n", "Total run time: 3.46s\n", "Total run time: 4.55s\n", "Total run time: 5.83s\n", "Total run time: 7.59s\n", "Total run time: 10.06s\n", "Total run time: 12.80s\n", "Total run time: 16.65s\n", "Total run time: 21.63s\n" ] } ], "source": [ "ntraj = 1\n", "def run_(**kwargs):\n", " epsilon = np.zeros(Nt)\n", " std = np.zeros(Nt)\n", " print(kwargs)\n", " for j in range(0,ntraj):\n", " for jj in range(0,Nt):\n", " Nsub = Nsubs[jj]\n", " sol = smesolve(H, rho0, tlist, c_op_td, sc_op, e_op, nsubsteps=Nsub, **kwargs)\n", " epsilon_j = 1/T * np.sum(np.abs(y_td - (sol.expect[1]-sol.expect[0]*sol.expect[0].conj())))*ddt\n", " epsilon[jj] += epsilon_j\n", " std[jj] += epsilon_j\n", " epsilon/= ntraj\n", " std = np.sqrt(1/ntraj * (1/ntraj * std - epsilon**2))\n", " return epsilon\n", "\n", "def get_stats_d1(**kw):\n", " start = time.time()\n", " y = run_(**kw)\n", " tag = str(kw[\"solver\"])\n", " x = np.log(stepsizes)\n", " ly = np.log(y)\n", " fit = np.polyfit(x, ly, 1)[0]\n", " return y,tag,fit,time.time()-start\n", "\n", "stat_d1 = []\n", "stat_d1.append(get_stats_d1(solver='euler-maruyama'))\n", "\n", "stat_d1.append(get_stats_d1(solver='platen'))\n", "stat_d1.append(get_stats_d1(solver='pc-euler'))\n", "stat_d1.append(get_stats_d1(solver='milstein'))\n", "stat_d1.append(get_stats_d1(solver='milstein-imp'))\n", "stat_d1.append(get_stats_d1(solver='pc-euler-2'))\n", "\n", "stat_d1.append(get_stats_d1(solver='explicit1.5'))\n", "stat_d1.append(get_stats_d1(solver=\"taylor1.5\"))\n", "stat_d1.append(get_stats_d1(solver=\"taylor1.5-imp\", args={\"tol\":1e-8}))\n", "\n", "stat_d1.append(get_stats_d1(solver=\"taylor2.0\"))\n", "stat_d1.append(get_stats_d1(solver=\"taylor2.0\", noiseDepth=500))" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n" ], "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "fig = plt.figure()\n", "ax = plt.subplot(111)\n", "\n", "mark = \"o*vspx+^<>1hdD\"\n", "\n", "for i,run in enumerate(stat_d1):\n", " ax.loglog(stepsizes, run[0], mark[i], label=run[1]+\": \" + str(run[2]))\n", "\n", "ax.loglog(stepsizes, 0.1*np.array(stepsizes)**0.5, label=\"$\\propto\\Delta t^{1/2}$\")\n", "ax.loglog(stepsizes, 0.05*np.array(stepsizes)**1, label=\"$\\propto\\Delta t$\")\n", "ax.loglog(stepsizes, 0.12*np.array(stepsizes)**1.5, label=\"$\\propto\\Delta t^{3/2}$\")\n", "ax.loglog(stepsizes, 0.7*np.array(stepsizes)**2.0, label=\"$\\propto\\Delta t^{2}$\")\n", "\n", "\n", "lgd=ax.legend(loc='center left', bbox_to_anchor=(1, 0.64), prop={'size':12})" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Both d1 and d2 time-dependent\n", "Using a taylor simulation with large nsubsteps instead of analytical solution." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Total run time: 90.46s\n" ] }, { "data": { "text/plain": [ "[]" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/svg+xml": [ "\n", "\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n" ], "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "def f(t, args):\n", " return 0.5+0.25*t-t*t*0.125\n", "\n", "Nsubs = (15*np.logspace(0,0.8,10)).astype(int)\n", "stepsizes = [ddt/j for j in Nsubs] # step size is doubled after each evaluation \n", "Nt = len(Nsubs) # number of step sizes that we compare\n", "Nsubmax = Nsubs[-1] # Number of intervals for the smallest step size; \n", "dtmin = (tlist[1]-tlist[0])/(Nsubmax)\n", "\n", "sc_op_td = [[sc_op[0],f]]\n", "sol = smesolve(H, rho0, tlist, c_op_td, sc_op_td, e_op, nsubsteps=1000, method=\"homodyne\",solver=\"taylor15\")\n", "y_btd = sol.expect[1]-sol.expect[0]*sol.expect[0].conj()\n", "plt.plot(y_btd)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "{'solver': 'euler-maruyama'}\n", "Total run time: 0.24s\n", "Total run time: 0.39s\n", "Total run time: 0.35s\n", "Total run time: 0.43s\n", "Total run time: 0.62s\n", "Total run time: 0.63s\n", "Total run time: 0.79s\n", "Total run time: 0.92s\n", "Total run time: 1.14s\n", "Total run time: 1.38s\n", "{'solver': 'platen'}\n", "Total run time: 0.55s\n", "Total run time: 0.65s\n", "Total run time: 0.84s\n", "Total run time: 1.28s\n", "Total run time: 1.24s\n", "Total run time: 1.77s\n", "Total run time: 3.18s\n", "Total run time: 2.64s\n", "Total run time: 4.25s\n", "Total run time: 4.50s\n", "{'solver': 'pc-euler'}\n", "Total run time: 0.39s\n", "Total run time: 0.45s\n", "Total run time: 0.55s\n", "Total run time: 0.67s\n", "Total run time: 0.84s\n", "Total run time: 1.01s\n", "Total run time: 1.28s\n", "Total run time: 1.52s\n", "Total run time: 1.88s\n", "Total run time: 2.65s\n", "{'solver': 'milstein'}\n", "Total run time: 0.31s\n", "Total run time: 0.36s\n", "Total run time: 0.43s\n", "Total run time: 0.53s\n", "Total run time: 0.66s\n", "Total run time: 0.79s\n", "Total run time: 0.98s\n", "Total run time: 1.19s\n", "Total run time: 1.67s\n", "Total run time: 1.96s\n", "{'solver': 'milstein-imp'}\n", "Total run time: 1.65s\n", "Total run time: 1.97s\n", "Total run time: 2.39s\n", "Total run time: 2.99s\n", "Total run time: 3.64s\n", "Total run time: 4.49s\n", "Total run time: 5.42s\n", "Total run time: 6.54s\n", "Total run time: 8.10s\n", "Total run time: 10.49s\n", "{'solver': 'pc-euler-2'}\n", "Total run time: 0.75s\n", "Total run time: 0.72s\n", "Total run time: 0.95s\n", "Total run time: 1.07s\n", "Total run time: 1.35s\n", "Total run time: 1.63s\n", "Total run time: 2.03s\n", "Total run time: 2.63s\n", "Total run time: 3.05s\n", "Total run time: 3.71s\n", "{'solver': 'explicit1.5'}\n", "Total run time: 1.18s\n", "Total run time: 1.45s\n", "Total run time: 1.56s\n", "Total run time: 1.94s\n", "Total run time: 2.40s\n", "Total run time: 2.89s\n", "Total run time: 3.57s\n", "Total run time: 4.35s\n", "Total run time: 5.74s\n", "Total run time: 7.07s\n", "{'solver': 'taylor1.5'}\n", "Total run time: 1.02s\n", "Total run time: 1.22s\n", "Total run time: 1.48s\n", "Total run time: 1.81s\n", "Total run time: 2.27s\n", "Total run time: 2.73s\n", "Total run time: 3.39s\n", "Total run time: 4.13s\n", "Total run time: 5.10s\n", "Total run time: 6.24s\n", "{'solver': 'taylor1.5-imp', 'args': {'tol': 2e-09}}\n", "Total run time: 2.28s\n", "Total run time: 2.67s\n", "Total run time: 3.40s\n", "Total run time: 3.91s\n", "Total run time: 4.80s\n", "Total run time: 6.62s\n", "Total run time: 9.73s\n", "Total run time: 9.33s\n", "Total run time: 12.57s\n", "Total run time: 15.14s\n" ] } ], "source": [ "ntraj = 1\n", "\n", "def run_(**kwargs):\n", " epsilon = np.zeros(Nt)\n", " std = np.zeros(Nt)\n", " print(kwargs)\n", " for j in range(0,ntraj):\n", " for jj in range(0,Nt):\n", " Nsub = Nsubs[jj]\n", " sol = smesolve(H, rho0, tlist, c_op_td, sc_op_td, e_op, nsubsteps=Nsub, **kwargs)\n", " epsilon_j = 1/T * np.sum(np.abs(y_btd - (sol.expect[1]-sol.expect[0]*sol.expect[0].conj())))*ddt\n", " epsilon[jj] += epsilon_j\n", " std[jj] += epsilon_j\n", " epsilon/= ntraj\n", " std = np.sqrt(1/ntraj * (1/ntraj * std - epsilon**2))\n", " return epsilon\n", "\n", "def get_stats_d2(**kw):\n", " start = time.time()\n", " y = run_(**kw)\n", " tag = str(kw[\"solver\"])\n", " x = np.log(stepsizes)\n", " ly = np.log(y)\n", " fit = np.polyfit(x, ly, 1)[0]\n", " return y,tag,fit,time.time()-start\n", "\n", "stats_d2 = []\n", "stats_d2.append(get_stats_d2(solver='euler-maruyama'))\n", "\n", "stats_d2.append(get_stats_d2(solver='platen'))\n", "stats_d2.append(get_stats_d2(solver='pc-euler'))\n", "stats_d2.append(get_stats_d2(solver='milstein'))\n", "stats_d2.append(get_stats_d2(solver='milstein-imp'))\n", "stats_d2.append(get_stats_d2(solver='pc-euler-2'))\n", "\n", "stats_d2.append(get_stats_d2(solver='explicit1.5'))\n", "stats_d2.append(get_stats_d2(solver='taylor1.5'))\n", "stats_d2.append(get_stats_d2(solver='taylor1.5-imp', args={\"tol\":2e-9}))" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n" ], "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "fig = plt.figure()\n", "ax = plt.subplot(111)\n", "\n", "mark = \"o*vspx+^<>1hdD\"\n", "\n", "for i,run in enumerate(stats_d2):\n", " ax.loglog(stepsizes, run[0], mark[i], label=run[1]+\": \" + str(run[2]))\n", "\n", "ax.loglog(stepsizes, 0.1*np.array(stepsizes)**0.5, label=\"$\\propto\\Delta t^{1/2}$\")\n", "ax.loglog(stepsizes, 0.05*np.array(stepsizes)**1, label=\"$\\propto\\Delta t$\")\n", "ax.loglog(stepsizes, 0.1*np.array(stepsizes)**1.5, label=\"$\\propto\\Delta t^{3/2}$\")\n", "\n", "\n", "lgd=ax.legend(loc='center left', bbox_to_anchor=(1, 0.64), prop={'size':12})" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Multiple sc_ops with time dependence" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Total run time: 160.70s\n" ] }, { "data": { "text/plain": [ "[]" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/svg+xml": [ "\n", "\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n" ], "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "def f(t, args):\n", " return 0.5+0.25*t-t*t*0.125\n", "\n", "def g(t, args):\n", " return 0.25+0.25*t-t*t*0.125\n", "\n", "Nsubs = (20*np.logspace(0,0.6,8)).astype(int)\n", "stepsizes = [ddt/j for j in Nsubs] # step size is doubled after each evaluation \n", "Nt = len(Nsubs) # number of step sizes that we compare\n", "Nsubmax = Nsubs[-1] # Number of intervals for the smallest step size; \n", "dtmin = (tlist[1]-tlist[0])/(Nsubmax)\n", "\n", "sc_op2_td = [[sc_op[0],f],[sc_op[0],g]]\n", "sol = smesolve(H, rho0, tlist, c_op_td, sc_op2_td, e_op, nsubsteps=1000, method=\"homodyne\",solver=152)\n", "y_btd2 = sol.expect[1]-sol.expect[0]*sol.expect[0].conj()\n", "plt.plot(y_btd2)" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "{'solver': 'euler-maruyama'}\n", "Total run time: 0.53s\n", "Total run time: 0.58s\n", "Total run time: 0.70s\n", "Total run time: 0.87s\n", "Total run time: 1.08s\n", "Total run time: 1.28s\n", "Total run time: 1.57s\n", "Total run time: 2.01s\n", "{'solver': 'platen'}\n", "Total run time: 1.52s\n", "Total run time: 1.82s\n", "Total run time: 2.21s\n", "Total run time: 2.78s\n", "Total run time: 3.33s\n", "Total run time: 4.05s\n", "Total run time: 4.97s\n", "Total run time: 6.01s\n", "{'solver': 'pc-euler'}\n", "Total run time: 0.92s\n", "Total run time: 1.10s\n", "Total run time: 1.33s\n", "Total run time: 1.64s\n", "Total run time: 2.01s\n", "Total run time: 2.41s\n", "Total run time: 3.05s\n", "Total run time: 3.60s\n", "{'solver': 'milstein'}\n", "Total run time: 0.73s\n", "Total run time: 0.88s\n", "Total run time: 1.06s\n", "Total run time: 1.32s\n", "Total run time: 1.62s\n", "Total run time: 1.96s\n", "Total run time: 2.37s\n", "Total run time: 2.88s\n", "{'solver': 'milstein-imp'}\n", "Total run time: 3.11s\n", "Total run time: 3.61s\n", "Total run time: 4.36s\n", "Total run time: 5.51s\n", "Total run time: 6.58s\n", "Total run time: 9.64s\n", "Total run time: 10.74s\n", "Total run time: 12.29s\n", "{'solver': 'pc-euler-2'}\n", "Total run time: 1.54s\n", "Total run time: 1.86s\n", "Total run time: 2.50s\n", "Total run time: 2.79s\n", "Total run time: 3.36s\n", "Total run time: 4.03s\n", "Total run time: 4.98s\n", "Total run time: 6.56s\n", "{'solver': 'explicit1.5'}\n", "Total run time: 5.28s\n", "Total run time: 6.35s\n", "Total run time: 7.74s\n", "Total run time: 8.99s\n", "Total run time: 10.34s\n", "Total run time: 12.41s\n", "Total run time: 15.40s\n", "Total run time: 21.40s\n", "{'solver': 'taylor1.5'}\n", "Total run time: 3.30s\n", "Total run time: 3.85s\n", "Total run time: 4.52s\n", "Total run time: 5.63s\n", "Total run time: 7.22s\n", "Total run time: 8.72s\n", "Total run time: 9.51s\n", "Total run time: 11.84s\n", "{'solver': 'taylor1.5-imp'}\n", "Total run time: 5.43s\n", "Total run time: 6.81s\n", "Total run time: 7.60s\n", "Total run time: 10.05s\n", "Total run time: 11.88s\n", "Total run time: 12.88s\n", "Total run time: 15.63s\n", "Total run time: 22.31s\n" ] } ], "source": [ "ntraj = 1\n", "\n", "def run_multi(**kwargs):\n", " epsilon = np.zeros(Nt)\n", " std = np.zeros(Nt)\n", " print(kwargs)\n", " for j in range(0,ntraj):\n", " for jj in range(0,Nt):\n", " Nsub = Nsubs[jj]\n", " sol = smesolve(H, rho0, tlist, c_op_td, sc_op2_td, e_op, nsubsteps=Nsub, **kwargs)\n", " epsilon_j = 1/T * np.sum(np.abs(y_btd2 - (sol.expect[1]-sol.expect[0]*sol.expect[0].conj())))*ddt\n", " epsilon[jj] += epsilon_j\n", " std[jj] += epsilon_j\n", " epsilon/= ntraj\n", " std = np.sqrt(1/ntraj * (1/ntraj * std - epsilon**2))\n", " return epsilon\n", "\n", "def get_stats_multi(**kw):\n", " start = time.time()\n", " y = run_multi(**kw)\n", " tag = str(kw[\"solver\"])\n", " x = np.log(stepsizes)\n", " ly = np.log(y)\n", " fit = np.polyfit(x, ly, 1)[0]\n", " return y,tag,fit,time.time()-start\n", "\n", "stats_multi = []\n", "stats_multi.append(get_stats_multi(solver='euler-maruyama'))\n", "\n", "stats_multi.append(get_stats_multi(solver='platen'))\n", "stats_multi.append(get_stats_multi(solver='pc-euler'))\n", "stats_multi.append(get_stats_multi(solver='milstein'))\n", "stats_multi.append(get_stats_multi(solver='milstein-imp'))\n", "stats_multi.append(get_stats_multi(solver='pc-euler-2'))\n", "\n", "stats_multi.append(get_stats_multi(solver='explicit1.5'))\n", "stats_multi.append(get_stats_multi(solver=\"taylor1.5\"))\n", "stats_multi.append(get_stats_multi(solver=\"taylor1.5-imp\"))" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n" ], "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "fig = plt.figure()\n", "ax = plt.subplot(111)\n", "\n", "mark = \"o*vspx+^<>1hdD\"\n", "\n", "for i,run in enumerate(stats_multi):\n", " ax.loglog(stepsizes, run[0], mark[i], label=run[1]+\": \" + str(run[2]))\n", "\n", "ax.loglog(stepsizes, 0.1*np.array(stepsizes)**0.5, label=\"$\\propto\\Delta t^{1/2}$\")\n", "ax.loglog(stepsizes, 0.1*np.array(stepsizes)**1, label=\"$\\propto\\Delta t$\")\n", "ax.loglog(stepsizes, 0.2*np.array(stepsizes)**1.5, label=\"$\\propto\\Delta t^{3/2}$\")\n", "\n", "lgd=ax.legend(loc='center left', bbox_to_anchor=(1, 0.64), prop={'size':12})" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Versions" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
SoftwareVersion
QuTiP4.4.0.dev0+1cf1dd3e
Numpy1.16.0
SciPy1.2.0
matplotlib3.0.2
Cython0.29.2
Number of CPUs2
BLAS InfoOPENBLAS
IPython7.2.0
Python3.6.7 (default, Oct 22 2018, 11:32:17) \n", "[GCC 8.2.0]
OSposix [linux]
Wed Jan 16 16:53:54 2019 JST
" ], "text/plain": [ "" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from qutip.ipynbtools import version_table\n", "\n", "version_table()" ] } ], "metadata": { "anaconda-cloud": {}, "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.7" } }, "nbformat": 4, "nbformat_minor": 1 }