{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Test for different solvers for stochastic equation\n", "\n", "Based on development-smesolver-new-methods by Manuel Grimm, Niels Lörch, and Denis V. Vasilyev.\n", "\n", "Eric Giguere, March 2018" ] }, { "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", " \n", " \n", " \n", " \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 = 200\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", "\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": "markdown", "metadata": {}, "source": [ "## Test of different SME solvers" ] }, { "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 = 200 # 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 = list((13*np.logspace(0,1,10)).astype(np.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": [ "## Plotting the figure - Constant case" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Total run time: 164.50s\n" ] } ], "source": [ "# Analetical solution not available: \n", "# Compute the evolution with the best solver and very small step size and use it as the reference\n", "sol = ssesolve(H, fock(N), tlist, [sc_op[0]+c_op[0]], e_op, nsubsteps=2000, method=\"homodyne\",solver=\"taylor2.0\")\n", "y_sse = sol.expect[1]-sol.expect[0]*sol.expect[0].conj()" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "{'solver': 'euler-maruyama'}\n", "Total run time: 0.14s\n", "Total run time: 0.15s\n", "Total run time: 0.18s\n", "Total run time: 0.23s\n", "Total run time: 0.30s\n", "Total run time: 0.38s\n", "Total run time: 0.49s\n", "Total run time: 0.64s\n", "Total run time: 0.81s\n", "Total run time: 1.04s\n", "{'solver': 'platen'}\n", "Total run time: 0.44s\n", "Total run time: 0.50s\n", "Total run time: 0.63s\n", "Total run time: 0.82s\n", "Total run time: 1.06s\n", "Total run time: 1.34s\n", "Total run time: 1.75s\n", "Total run time: 2.33s\n", "Total run time: 3.06s\n", "Total run time: 3.86s\n", "{'solver': 'pred-corr'}\n", "Total run time: 0.26s\n", "Total run time: 0.31s\n", "Total run time: 0.41s\n", "Total run time: 0.57s\n", "Total run time: 0.82s\n", "Total run time: 0.90s\n", "Total run time: 1.16s\n", "Total run time: 1.49s\n", "Total run time: 1.95s\n", "Total run time: 2.48s\n", "{'solver': 'milstein'}\n", "Total run time: 0.14s\n", "Total run time: 0.18s\n", "Total run time: 0.22s\n", "Total run time: 0.29s\n", "Total run time: 0.37s\n", "Total run time: 0.47s\n", "Total run time: 0.63s\n", "Total run time: 0.78s\n", "Total run time: 1.02s\n", "Total run time: 1.33s\n", "{'solver': 'milstein-imp', 'tol': 1e-09}\n", "Total run time: 1.39s\n", "Total run time: 1.69s\n", "Total run time: 2.18s\n", "Total run time: 2.85s\n", "Total run time: 3.60s\n", "Total run time: 4.56s\n", "Total run time: 5.94s\n", "Total run time: 7.37s\n", "Total run time: 9.63s\n", "Total run time: 12.50s\n", "{'solver': 'pred-corr-2'}\n", "Total run time: 0.29s\n", "Total run time: 0.35s\n", "Total run time: 0.45s\n", "Total run time: 0.60s\n", "Total run time: 0.79s\n", "Total run time: 0.99s\n", "Total run time: 1.29s\n", "Total run time: 1.64s\n", "Total run time: 2.14s\n", "Total run time: 2.80s\n", "{'solver': 'explicit1.5'}\n", "Total run time: 0.77s\n", "Total run time: 0.94s\n", "Total run time: 1.22s\n", "Total run time: 1.60s\n", "Total run time: 2.13s\n", "Total run time: 2.63s\n", "Total run time: 3.45s\n", "Total run time: 4.39s\n", "Total run time: 5.84s\n", "Total run time: 7.88s\n", "{'solver': 'taylor1.5'}\n", "Total run time: 0.68s\n", "Total run time: 0.84s\n", "Total run time: 1.18s\n", "Total run time: 1.49s\n", "Total run time: 2.09s\n", "Total run time: 2.44s\n", "Total run time: 4.34s\n", "Total run time: 5.30s\n", "Total run time: 6.14s\n", "Total run time: 7.20s\n", "{'solver': 'taylor1.5-imp', 'tol': 1e-09}\n", "Total run time: 2.18s\n", "Total run time: 3.24s\n", "Total run time: 3.53s\n", "Total run time: 3.63s\n", "Total run time: 5.55s\n", "Total run time: 7.30s\n", "Total run time: 9.34s\n", "Total run time: 11.63s\n", "Total run time: 15.92s\n", "Total run time: 17.30s\n", "{'solver': 'taylor2.0'}\n", "Total run time: 0.85s\n", "Total run time: 1.05s\n", "Total run time: 1.38s\n", "Total run time: 1.82s\n", "Total run time: 2.76s\n", "Total run time: 4.01s\n", "Total run time: 6.84s\n", "Total run time: 8.19s\n", "Total run time: 8.84s\n", "Total run time: 11.23s\n", "{'solver': 'taylor2.0', 'noiseDepth': 500}\n", "Total run time: 8.60s\n", "Total run time: 9.90s\n", "Total run time: 13.52s\n", "Total run time: 17.00s\n", "Total run time: 21.31s\n", "Total run time: 27.48s\n", "Total run time: 36.43s\n", "Total run time: 45.66s\n", "Total run time: 61.22s\n", "Total run time: 75.65s\n" ] } ], "source": [ "ntraj = 1\n", "def run_sse(**kwargs):\n", " epsilon = np.zeros(Nt)\n", " std = np.zeros(Nt)\n", " print(kwargs)\n", " for jj in range(0,Nt):\n", " for j in range(0,ntraj):\n", " Nsub = Nsubs[jj]#int(Nsubmax/(2**jj))\n", " sol = ssesolve(H, fock(N), tlist, [sc_op[0]+c_op[0]], e_op, nsubsteps=Nsub, **kwargs)\n", " epsilon_j = 1/T * np.sum(np.abs(y_sse - (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(**kw):\n", " start = time.time()\n", " y = run_sse(**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 = []\n", "stats_cte.append(get_stats(solver='euler-maruyama'))\n", "\n", "stats_cte.append(get_stats(solver='platen'))\n", "stats_cte.append(get_stats(solver='pred-corr'))\n", "stats_cte.append(get_stats(solver='milstein'))\n", "stats_cte.append(get_stats(solver='milstein-imp', tol=1e-9))\n", "stats_cte.append(get_stats(solver='pred-corr-2'))\n", "\n", "stats_cte.append(get_stats(solver='explicit1.5'))\n", "stats_cte.append(get_stats(solver=\"taylor1.5\"))\n", "stats_cte.append(get_stats(solver=\"taylor1.5-imp\", tol=1e-9))\n", "\n", "stats_cte.append(get_stats(solver=\"taylor2.0\"))\n", "stats_cte.append(get_stats(solver=\"taylor2.0\", noiseDepth=500))" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "scrolled": false }, "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" ], "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "fig = plt.figure()\n", "ax = plt.subplot(111)\n", "mark = \"o*vspx+^<>1hdD\"\n", "\n", "for i,run in enumerate(stats_cte):\n", " ax.loglog(stepsizes, run[0], mark[i], label=run[1]+\": \" + str(run[2]))\n", "\n", "ax.loglog(stepsizes, 0.003*np.array(stepsizes)**0.5, label=\"$\\propto\\Delta t^{1/2}$\")\n", "ax.loglog(stepsizes, 0.01*np.array(stepsizes)**1, label=\"$\\propto\\Delta t$\")\n", "ax.loglog(stepsizes, 0.001*np.array(stepsizes)**1, label=\"$\\propto\\Delta t$\")\n", "ax.loglog(stepsizes, 0.01*np.array(stepsizes)**1.5, label=\"$\\propto\\Delta t^{3/2}$\")\n", "ax.loglog(stepsizes, 0.05*np.array(stepsizes)**2.0, label=\"$\\propto\\Delta t^{2}$\")\n", "\n", "ax.set_xlabel(r'$\\Delta t$ $\\left[\\gamma^{-1}\\right]$')\n", "ax.set_ylabel('deviation')\n", "\n", "lgd=ax.legend(loc='center left', bbox_to_anchor=(1, 0.64), prop={'size':12})\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Deterministic part time dependent" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Total run time: 258.38s\n" ] }, { "data": { "text/plain": [ "[]" ] }, "execution_count": 7, "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", " \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 H_f(t,args):\n", " return 0.125+t/12+t*t/72\n", "sol = ssesolve([H,[c_op[0].dag()*c_op[0]/2,H_f]], fock(N), tlist, sc_op, e_op, \n", " nsubsteps=2500, method=\"homodyne\",solver=\"taylor2.0\")\n", "y_sse_td = sol.expect[1]-sol.expect[0]*sol.expect[0].conj()\n", "plt.plot(y_sse_td)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "{'solver': 'euler-maruyama'}\n", "Total run time: 0.18s\n", "Total run time: 0.16s\n", "Total run time: 0.16s\n", "Total run time: 0.28s\n", "Total run time: 0.34s\n", "Total run time: 0.44s\n", "Total run time: 0.55s\n", "Total run time: 0.70s\n", "Total run time: 0.90s\n", "Total run time: 1.17s\n", "{'solver': 'platen'}\n", "Total run time: 0.43s\n", "Total run time: 0.51s\n", "Total run time: 0.67s\n", "Total run time: 0.89s\n", "Total run time: 1.11s\n", "Total run time: 1.42s\n", "Total run time: 1.86s\n", "Total run time: 2.34s\n", "Total run time: 3.09s\n", "Total run time: 4.04s\n", "{'solver': 'pred-corr'}\n", "Total run time: 0.26s\n", "Total run time: 0.32s\n", "Total run time: 0.41s\n", "Total run time: 0.54s\n", "Total run time: 0.69s\n", "Total run time: 0.89s\n", "Total run time: 1.12s\n", "Total run time: 1.45s\n", "Total run time: 1.90s\n", "Total run time: 2.51s\n", "{'solver': 'milstein'}\n", "Total run time: 0.19s\n", "Total run time: 0.22s\n", "Total run time: 0.31s\n", "Total run time: 0.33s\n", "Total run time: 0.42s\n", "Total run time: 0.65s\n", "Total run time: 0.96s\n", "Total run time: 1.17s\n", "Total run time: 1.33s\n", "Total run time: 1.45s\n", "{'solver': 'milstein-imp'}\n", "Total run time: 1.62s\n", "Total run time: 1.93s\n", "Total run time: 2.90s\n", "Total run time: 3.57s\n", "Total run time: 5.41s\n", "Total run time: 5.80s\n", "Total run time: 7.49s\n", "Total run time: 9.71s\n", "Total run time: 12.71s\n", "Total run time: 15.76s\n", "{'solver': 'pred-corr-2'}\n", "Total run time: 0.39s\n", "Total run time: 0.54s\n", "Total run time: 0.58s\n", "Total run time: 0.81s\n", "Total run time: 1.00s\n", "Total run time: 1.36s\n", "Total run time: 1.90s\n", "Total run time: 2.62s\n", "Total run time: 3.72s\n", "Total run time: 4.37s\n", "{'solver': 'explicit1.5'}\n", "Total run time: 1.26s\n", "Total run time: 2.17s\n", "Total run time: 2.90s\n", "Total run time: 2.62s\n", "Total run time: 3.31s\n", "Total run time: 4.27s\n", "Total run time: 4.98s\n", "Total run time: 6.04s\n", "Total run time: 8.11s\n", "Total run time: 10.18s\n", "{'solver': 'taylor1.5'}\n", "Total run time: 0.93s\n", "Total run time: 0.98s\n", "Total run time: 1.74s\n", "Total run time: 2.24s\n", "Total run time: 2.82s\n", "Total run time: 3.78s\n", "Total run time: 4.68s\n", "Total run time: 5.86s\n", "Total run time: 7.55s\n", "Total run time: 9.36s\n", "{'solver': 'taylor1.5-imp', 'tol': 1e-09}\n", "Total run time: 2.61s\n", "Total run time: 2.87s\n", "Total run time: 3.77s\n", "Total run time: 4.93s\n", "Total run time: 6.18s\n", "Total run time: 7.98s\n", "Total run time: 9.44s\n", "Total run time: 14.02s\n", "Total run time: 16.98s\n", "Total run time: 21.73s\n", "{'solver': 'taylor2.0'}\n", "Total run time: 1.43s\n", "Total run time: 1.78s\n", "Total run time: 1.93s\n", "Total run time: 2.44s\n", "Total run time: 3.37s\n", "Total run time: 4.20s\n", "Total run time: 5.96s\n", "Total run time: 7.64s\n", "Total run time: 9.52s\n", "Total run time: 13.98s\n", "{'solver': 'taylor2.0', 'noiseDepth': 500}\n", "Total run time: 8.99s\n", "Total run time: 9.84s\n", "Total run time: 12.45s\n", "Total run time: 17.20s\n", "Total run time: 21.36s\n", "Total run time: 28.13s\n", "Total run time: 38.48s\n", "Total run time: 49.10s\n", "Total run time: 65.57s\n", "Total run time: 79.84s\n" ] } ], "source": [ "ntraj = 1\n", "def run_sse_td(**kwargs):\n", " epsilon = np.zeros(Nt)\n", " std = np.zeros(Nt)\n", " print(kwargs)\n", " for jj in range(0,Nt):\n", " for j in range(0,ntraj):\n", " Nsub = Nsubs[jj]#int(Nsubmax/(2**jj))\n", " sol = ssesolve([H,[c_op[0].dag()*c_op[0]/2,H_f]], fock(N), tlist, sc_op, e_op, nsubsteps=Nsub, **kwargs)\n", " epsilon_j = 1/T * np.sum(np.abs(y_sse_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(**kw):\n", " y = run_sse_td(**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\n", "\n", "\n", "stats_td = []\n", "stats_td.append(get_stats(solver='euler-maruyama'))\n", "\n", "stats_td.append(get_stats(solver='platen'))\n", "stats_td.append(get_stats(solver='pred-corr'))\n", "stats_td.append(get_stats(solver='milstein'))\n", "stats_td.append(get_stats(solver='milstein-imp'))\n", "stats_td.append(get_stats(solver='pred-corr-2'))\n", "\n", "stats_td.append(get_stats(solver='explicit1.5'))\n", "stats_td.append(get_stats(solver=\"taylor1.5\"))\n", "stats_td.append(get_stats(solver=\"taylor1.5-imp\", tol=1e-9))\n", "\n", "stats_td.append(get_stats(solver=\"taylor2.0\"))\n", "stats_td.append(get_stats(solver=\"taylor2.0\", noiseDepth=500))" ] }, { "cell_type": "code", "execution_count": 9, "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", " \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_td):\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.5*np.array(stepsizes)**2.0, label=\"$\\propto\\Delta t^{2}$\")\n", "\n", "ax.set_xlabel(r'$\\Delta t$ $\\left[\\gamma^{-1}\\right]$')\n", "ax.set_ylabel('deviation')\n", "\n", "lgd=ax.legend(loc='center left', bbox_to_anchor=(1, 0.64), prop={'size':12})" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "## Both d1 and d2 time-dependent" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Total run time: 76.28s\n" ] }, { "data": { "text/plain": [ "[]" ] }, "execution_count": 17, "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", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \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 H_f(t,args):\n", " return 0.125+t/12+t*t/72\n", "def H_bf(t,args):\n", " return 0.125+t/10+t*t/108\n", "sc_op_td = [[sc_op[0],H_bf]]\n", "sol = ssesolve([H,[c_op[0].dag()*c_op[0]/2,H_f]], fock(N), tlist, sc_op_td, e_op, \n", " nsubsteps=2000, method=\"homodyne\",solver=\"taylor15\")\n", "y_sse_btd = sol.expect[1]-sol.expect[0]*sol.expect[0].conj()\n", "plt.plot(y_sse_btd)" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "{'solver': 'euler-maruyama'}\n", "Total run time: 0.25s\n", "Total run time: 0.15s\n", "Total run time: 0.28s\n", "Total run time: 0.24s\n", "Total run time: 0.33s\n", "Total run time: 0.41s\n", "Total run time: 0.52s\n", "Total run time: 0.71s\n", "Total run time: 0.85s\n", "Total run time: 1.05s\n", "{'solver': 'platen'}\n", "Total run time: 0.34s\n", "Total run time: 0.41s\n", "Total run time: 0.56s\n", "Total run time: 1.14s\n", "Total run time: 0.90s\n", "Total run time: 1.15s\n", "Total run time: 1.47s\n", "Total run time: 1.89s\n", "Total run time: 2.41s\n", "Total run time: 3.17s\n", "{'solver': 'pred-corr'}\n", "Total run time: 0.17s\n", "Total run time: 0.24s\n", "Total run time: 0.29s\n", "Total run time: 0.35s\n", "Total run time: 0.46s\n", "Total run time: 0.58s\n", "Total run time: 0.74s\n", "Total run time: 0.99s\n", "Total run time: 1.27s\n", "Total run time: 1.69s\n", "{'solver': 'milstein'}\n", "Total run time: 0.11s\n", "Total run time: 0.16s\n", "Total run time: 0.16s\n", "Total run time: 0.21s\n", "Total run time: 0.27s\n", "Total run time: 0.34s\n", "Total run time: 0.44s\n", "Total run time: 0.56s\n", "Total run time: 0.72s\n", "Total run time: 0.96s\n", "{'solver': 'milstein-imp'}\n", "Total run time: 0.79s\n", "Total run time: 0.98s\n", "Total run time: 1.27s\n", "Total run time: 1.74s\n", "Total run time: 2.13s\n", "Total run time: 3.02s\n", "Total run time: 3.65s\n", "Total run time: 4.57s\n", "Total run time: 5.86s\n", "Total run time: 7.27s\n", "{'solver': 'pred-corr-2'}\n", "Total run time: 0.25s\n", "Total run time: 0.27s\n", "Total run time: 0.37s\n", "Total run time: 0.47s\n", "Total run time: 0.59s\n", "Total run time: 0.75s\n", "Total run time: 1.02s\n", "Total run time: 1.26s\n", "Total run time: 1.61s\n", "Total run time: 2.07s\n", "{'solver': 'explicit1.5'}\n", "Total run time: 0.70s\n", "Total run time: 0.81s\n", "Total run time: 1.07s\n", "Total run time: 1.44s\n", "Total run time: 1.84s\n", "Total run time: 2.32s\n", "Total run time: 3.06s\n", "Total run time: 3.85s\n", "Total run time: 5.06s\n", "Total run time: 6.56s\n", "{'solver': 'taylor1.5'}\n", "Total run time: 0.46s\n", "Total run time: 0.56s\n", "Total run time: 0.82s\n", "Total run time: 0.96s\n", "Total run time: 1.26s\n", "Total run time: 1.61s\n", "Total run time: 2.06s\n", "Total run time: 2.67s\n", "Total run time: 3.43s\n", "Total run time: 4.49s\n", "{'solver': 'taylor1.5-imp', 'tol': 1e-09}\n", "Total run time: 1.35s\n", "Total run time: 1.39s\n", "Total run time: 1.81s\n", "Total run time: 2.25s\n", "Total run time: 2.87s\n", "Total run time: 3.76s\n", "Total run time: 4.67s\n", "Total run time: 6.00s\n", "Total run time: 8.13s\n", "Total run time: 11.06s\n" ] } ], "source": [ "ntraj = 1\n", "def run_sse_td(**kwargs):\n", " epsilon = np.zeros(Nt)\n", " std = np.zeros(Nt)\n", " print(kwargs)\n", " for jj in range(0,Nt):\n", " for j in range(0,ntraj):\n", " Nsub = Nsubs[jj]#int(Nsubmax/(2**jj))\n", " sol = ssesolve([H,[c_op[0].dag()*c_op[0]/2,H_f]], fock(N), tlist, sc_op_td, e_op, nsubsteps=Nsub, **kwargs)\n", " epsilon_j = 1/T * np.sum(np.abs(y_sse_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_b(**kw):\n", " y = run_sse_td(**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\n", "\n", "\n", "stats_d2_td = []\n", "stats_d2_td.append(get_stats_b(solver='euler-maruyama'))\n", "\n", "stats_d2_td.append(get_stats_b(solver='platen'))\n", "stats_d2_td.append(get_stats_b(solver='pred-corr'))\n", "stats_d2_td.append(get_stats_b(solver='milstein'))\n", "stats_d2_td.append(get_stats_b(solver='milstein-imp'))\n", "stats_d2_td.append(get_stats_b(solver='pred-corr-2'))\n", "\n", "stats_d2_td.append(get_stats_b(solver='explicit1.5'))\n", "stats_d2_td.append(get_stats_b(solver=\"taylor1.5\"))\n", "stats_d2_td.append(get_stats_b(solver=\"taylor1.5-imp\", tol=1e-9))" ] }, { "cell_type": "code", "execution_count": 19, "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" ], "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "fig = plt.figure()\n", "ax = plt.subplot(111)\n", "mark = \"o*vspx+^<>1hdD\"\n", "\n", "for i,run in enumerate(stats_d2_td):\n", " ax.loglog(stepsizes, run[0], mark[i], label=run[1]+\": \" + str(run[2]))\n", "\n", "ax.loglog(stepsizes, 0.03*np.array(stepsizes)**0.5, label=\"$\\propto\\Delta t^{1/2}$\")\n", "ax.loglog(stepsizes, 0.03*np.array(stepsizes)**1, label=\"$\\propto\\Delta t$\")\n", "ax.loglog(stepsizes, 0.03*np.array(stepsizes)**1.5, label=\"$\\propto\\Delta t^{3/2}$\")\n", "\n", "ax.set_xlabel(r'$\\Delta t$ $\\left[\\gamma^{-1}\\right]$')\n", "ax.set_ylabel('deviation')\n", "\n", "lgd=ax.legend(loc='center left', bbox_to_anchor=(1, 0.64), prop={'size':12})" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "## Multiple sc_ops, time-dependent" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Total run time: 343.58s\n" ] }, { "data": { "text/plain": [ "[]" ] }, "execution_count": 13, "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", " \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 H_f(t,args):\n", " return 0.125+t/12+t*t/36\n", "def H_bf(t,args):\n", " return 0.125+t/10+t*t/108\n", "sc_op_td = [[sc_op[0]],[sc_op[0],H_bf],[sc_op[0],H_f]]\n", "\n", "sol = ssesolve([H,[c_op[0].dag()*c_op[0]/2,H_f]], fock(N), tlist/3, sc_op_td, e_op, \n", " nsubsteps=2000, method=\"homodyne\",solver=\"taylor15\")\n", "y_sse_multi = sol.expect[1]-sol.expect[0]*sol.expect[0].conj()\n", "plt.plot(y_sse_multi)" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "{'solver': 'euler-maruyama'}\n", "Total run time: 0.41s\n", "Total run time: 0.56s\n", "Total run time: 0.40s\n", "Total run time: 0.75s\n", "Total run time: 0.70s\n", "Total run time: 0.90s\n", "Total run time: 1.62s\n", "Total run time: 2.15s\n", "Total run time: 3.35s\n", "Total run time: 3.28s\n", "{'solver': 'platen'}\n", "Total run time: 1.52s\n", "Total run time: 2.28s\n", "Total run time: 2.07s\n", "Total run time: 3.27s\n", "Total run time: 3.71s\n", "Total run time: 4.70s\n", "Total run time: 7.20s\n", "Total run time: 9.65s\n", "Total run time: 11.76s\n", "Total run time: 14.92s\n", "{'solver': 'pred-corr'}\n", "Total run time: 0.57s\n", "Total run time: 0.70s\n", "Total run time: 1.18s\n", "Total run time: 1.33s\n", "Total run time: 2.01s\n", "Total run time: 2.24s\n", "Total run time: 2.71s\n", "Total run time: 4.29s\n", "Total run time: 5.56s\n", "Total run time: 6.50s\n", "{'solver': 'milstein'}\n", "Total run time: 0.44s\n", "Total run time: 0.51s\n", "Total run time: 0.66s\n", "Total run time: 1.06s\n", "Total run time: 1.28s\n", "Total run time: 1.54s\n", "Total run time: 1.91s\n", "Total run time: 2.66s\n", "Total run time: 3.54s\n", "Total run time: 4.95s\n", "{'solver': 'milstein-imp'}\n", "Total run time: 1.31s\n", "Total run time: 1.58s\n", "Total run time: 2.31s\n", "Total run time: 2.92s\n", "Total run time: 4.14s\n", "Total run time: 5.03s\n", "Total run time: 5.96s\n", "Total run time: 9.10s\n", "Total run time: 10.52s\n", "Total run time: 11.72s\n", "{'solver': 'pred-corr-2'}\n", "Total run time: 0.81s\n", "Total run time: 1.08s\n", "Total run time: 1.33s\n", "Total run time: 1.92s\n", "Total run time: 2.73s\n", "Total run time: 3.27s\n", "Total run time: 5.21s\n", "Total run time: 6.45s\n", "Total run time: 7.51s\n", "Total run time: 10.05s\n", "{'solver': 'explicit1.5'}\n", "Total run time: 5.17s\n", "Total run time: 6.66s\n", "Total run time: 8.43s\n", "Total run time: 11.16s\n", "Total run time: 14.72s\n", "Total run time: 36.61s\n", "Total run time: 51.42s\n", "Total run time: 62.39s\n", "Total run time: 49.08s\n", "Total run time: 93.40s\n", "{'solver': 'taylor1.5'}\n", "Total run time: 2.36s\n", "Total run time: 2.72s\n", "Total run time: 5.01s\n", "Total run time: 10.42s\n", "Total run time: 11.30s\n", "Total run time: 18.45s\n", "Total run time: 14.06s\n", "Total run time: 13.55s\n", "Total run time: 19.14s\n", "Total run time: 21.60s\n", "{'solver': 'taylor1.5-imp', 'tol': 1e-09}\n", "Total run time: 2.45s\n", "Total run time: 2.94s\n", "Total run time: 4.47s\n", "Total run time: 5.56s\n", "Total run time: 7.15s\n", "Total run time: 9.24s\n", "Total run time: 12.75s\n", "Total run time: 15.05s\n", "Total run time: 20.00s\n", "Total run time: 27.10s\n" ] } ], "source": [ "ntraj = 1\n", "def run_sss_multi(**kwargs):\n", " epsilon = np.zeros(Nt)\n", " std = np.zeros(Nt)\n", " print(kwargs)\n", " for jj in range(0,Nt):\n", " for j in range(0,ntraj):\n", " Nsub = Nsubs[jj]#int(Nsubmax/(2**jj))\n", " sol = ssesolve([H,[c_op[0].dag()*c_op[0]/2,H_f]], fock(N), tlist/3, sc_op_td, e_op, nsubsteps=Nsub, **kwargs)\n", " epsilon_j = 1/T * np.sum(np.abs(y_sse_multi - (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", " \n", " return epsilon\n", "\n", "def get_stats_multi(**kw):\n", " y = run_sss_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)\n", "\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='pred-corr'))\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='pred-corr-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\", tol=1e-9))" ] }, { "cell_type": "code", "execution_count": 15, "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" ], "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "fig = plt.figure()\n", "ax = plt.subplot(111)\n", "mark = \"o*vspx+^<>Dd\"\n", "\n", "for run in stats_multi:\n", " ax.loglog(stepsizes, run[0], 'o', label=run[1]+\": \" + str(run[2]))\n", "\n", "ax.loglog(stepsizes, 0.05*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.05*np.array(stepsizes)**1.5, label=\"$\\propto\\Delta t^{3/2}$\")\n", "\n", "ax.set_xlabel(r'$\\Delta t$ $\\left[\\gamma^{-1}\\right]$')\n", "ax.set_ylabel('deviation')\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": 16, "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 15:58:27 2019 JST
" ], "text/plain": [ "" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from qutip.ipynbtools import version_table\n", "\n", "version_table()" ] }, { "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.7" } }, "nbformat": 4, "nbformat_minor": 2 }