mkprop
by Tobias Jawecki
This package provides some Python routines for adaptive time integration based on Magnus-Krylov and Krylov methods.
This package is still under development. Future versions might have different API. The current version is not fully documented.
Getting started
Install the latest release from the pypi package manager python -m pip install mkprop
Mathematical problem description
Let
In a similar manner, we also consider non-autonomous system of ODE's.
The time propagation for the solution
Krylov methods
The routine expimv_pKry
provides an approximation to the action of the matrix exponential on a vector
usage: examples/basicKrylov.ipynb
import scipy.sparse
import scipy.linalg
import numpy as np
import mkprop
nrm = lambda x : np.linalg.norm(x,2)
n = 1000
e = np.ones(n)
e1 = np.ones(n-1)+1e-5*np.random.rand(n-1)
# M is finite difference discretization of 1d Laplace operator
M = scipy.sparse.diags([e1.conj(),-2*e,e1], [-1,0,1])
u=np.random.rand(n)
u=u/nrm(u)
dt = 50
tol = 1e-6
yref = scipy.sparse.linalg.expm_multiply(1j*dt*M,u)
y,_,_,mused = mkprop.expimv_pKry(M,u,t=dt,tol=tol)
print("approximation error = %.2e, tolerance = %.2e" % (nrm(yref-y)/dt, tol))
# output: approximation error = 5.06e-08, tolerance = 1.00e-06
Magnus integrators
The solution of the non-autonomous system of ODE's is given by a Magnus expansion
Commutator free Magnus (CFM) integrators are of the form
The following methods are available: CFM integrators with table of coefficients
Currently, only methods of order
Defining a problem with a time-dependent Hamiltonian H(t)
Define a simple problem, see also examples/exlaser.py
.
import numpy as np
import sympy as sym
import scipy.sparse
class doublewellproblem():
def __init__(self,n):
self.n = n
L = 5
self.x = np.linspace(-L,L,n+1)
# problem related inner product
dx = 2*L/(n+1)
self.inr = lambda x,y : dx*np.vdot(x,y)
# H_0
e = np.ones(n+1)
e1 = np.ones(n)
D2 = dx**(-2)*scipy.sparse.diags([e1.conj(),-2*e,e1], [-1,0,1])
self.H0 = D2
# a double well potential V_0
V0 = self.x**4 - 20*self.x**2
# add a time-dependent potential V(t)
# prepare V(t) and dV(t)
T = 5.0
omega = 10
x = sym.Symbol('x')
t = sym.Symbol('t')
exprV = 10*sym.sin((np.pi*t/T)**2)*sym.sin(omega*t)*x
exprdV = sym.diff(exprV, t)
evalVt = sym.lambdify([x,t], exprV, "numpy")
evaldVt = sym.lambdify([x,t], exprdV, "numpy")
self.Vt = lambda tau: (V0 + evalVt(self.x,tau))
self.dVt = lambda tau: evaldVt(self.x,tau)
def getprop(self):
# return nodes x and inner product
return self.x, self.inr
def getinitialstate(self):
# return initial state
s, x0 = 0.2, -2.5
u = (s*np.pi)**(-0.25)*np.exp(-(self.x-x0)**2/(2*s))
return u
def setupHamiltonian(self,t):
# return a routine to apply H(t) = H0 - (V0 + V(t))
# and dH(t) = - dV(t)
self.V = self.Vt(t)
self.dV = self.dVt(t)
H = self.H0 - scipy.sparse.diags([self.V], [0])
mv = lambda u : H.dot(u)
dH = -self.dV
dmv = lambda u : dH*u
return mv, dmv
def setupHamiltonianCFM(self,a,c,chat,t,dt):
# return a routine to apply sum_j aj*H(t + cj*dt)
# and dH(t) = sum_j (c_j+chat)*aj*dH(t + cj*dt)
# for some scalar chat
jexps = len(a)
V_CFM = a[0]*self.Vt(t+c[0]*dt)
dV_CFM = (c[0]+chat)*a[0]*self.dVt(t+c[0]*dt)
for j in range(jexps-1):
V_CFM += a[j+1]*self.Vt(t+c[j+1]*dt)
dV_CFM += (c[j+1]+chat)*a[j+1]*self.dVt(t+c[j+1]*dt)
self.V = V_CFM
self.dV = dV_CFM
if sum(a)!=0:
H = sum(a)*self.H0 - scipy.sparse.diags([self.V], [0])
mv = lambda u : H.dot(u)
else:
mv = lambda u : -self.V*u
dmv = lambda u : -self.dV*u
return mv, dmv
Magnus-Krylov methods
This package provides adaptive Magnus-Krylov methods, namely, using the adaptive midpoint rule and CFM integrators with error estimates based on symmetrized defects and works of Auzinger et al.. Again, Magnus-Krylov approximations
usage: examples/basicMagnusKrylov.ipynb
import numpy as np
import mkprop
from exlaser import doublewellproblem as prob
import matplotlib.pyplot as plt
# setup problem from exlaser.py
n=1200
Hamiltonian = prob(n)
x, inr = Hamiltonian.getprop()
nrm = lambda u : ((inr(u,u)).real)**0.5
u = Hamiltonian.getinitialstate()
# define initial and final time
tnow = 0
tend = 0.1
tau = tend-tnow
# some options for Krylov method, use Lanczos
ktype = 2
mmax = 60
# compute reference solution with adaptive fourth order CFM integrator
tolref=1e-6
dtinitref = 5e-2
yref,_,_,_,_,_,_ = mkprop.adaptiveCFMp4j2(u,tnow,tend,dtinitref,Hamiltonian,tol=tolref,
m=mmax,ktype=ktype,inr=inr)
# test adaptive midpoint rule
tol=1e-4
dtinit = 1e-3
y,_,_,_,_,_,_ = mkprop.adaptivemidpoint(u,tnow,tend,dtinit,Hamiltonian,tol=tol,
m=mmax,ktype=ktype,inr=inr)
print("approximation error = %.2e, tolerance = %.2e" % (nrm(yref-y)/tau, tol))
# output: approximation error = 6.40e-05, tolerance = 1.00e-04
Examples
examples/exlaser.py
: define Hamiltonian, etc..
examples/comparestepscosts.ipyn
: compare step sizes and computational cost of adaptive Magnus-Krylov methods. Step sizes:
Cost:
examples/testerrasymorder.ipynb
:
Errors and error estimates over the step size, asymptotic order.
Some references
A. Alvermann and H. Fehske. High-order commutator-free exponential time-propagation of driven quantum systems. J. Comput. Phys., 230(15):5930-5956, 2011. doi:10.1016/j.jcp.2011.04.006.
W. Auzinger, H. Hofstätter, and O. Koch. Symmetrized local error estimators for time-reversible one-step methods in nonlinear evolution equations. J. Comput. Appl. Math., 356:339-357, 2019. doi:10.1016/j.cam.2019.02.011.
W. Auzinger, H. Hofstätter, O. Koch, M. Quell, and M. Thalhammer. A posteriori error estimation for Magnus-type integrators. M2AN - Math. Model. Numer. Anal., 53(1):197-218, 2019. doi:10.1051/m2an/2018050.
W. Auzinger and O. Koch. An improved local error estimator for symmetric time-stepping schemes. Appl. Math. Lett., 82:106-110, 2018. doi:10.1016/j.aml.2018.03.001.
P. Bader, S. Blanes, and N. Kopylov. Exponential propagators for the Schrödinger equation with a time-dependent potential. J. Chem. Phys., 148(24):244109, 2018. doi:10.1063/1.5036838.
T. Jawecki, W. Auzinger, and O. Koch.
Computable upper error bounds for Krylov approximations to matrix exponentials and associated
T. Jawecki. Krylov techniques and approximations to the action of matrix exponentials. PhD thesis, TU Wien, Austria, 2022. doi:10.34726/hss.2022.45083.
T. Jawecki.
A study of defect-based error estimates for the Krylov approximation of
W. Magnus. On the exponential solution of differential equations for a linear operator. Comm. Pure Appl. Math., 7(4):649-673, 1954. doi:10.1002/cpa.3160070404.
S. Blanes, F. Casas, J. Oteo and J. Ros. The Magnus expansion and some of its applications. Phys. Rep., 470(5):151-238, 2009. doi:10.1016/j.physrep.2008.11.001.