← Back to branch summary

~jrjohansson/qutip/master

« back to all changes in this revision

Viewing changes to qutip/steady.py

  • Committer: Paul Nation
  • Date: 2011-04-21 04:46:56 UTC
  • Revision ID: git-v1:dd4c966b490aa468dfbd28cef66694df4bf235c8

Show diffs side-by-side

added added

removed removed

Lines of Context:
qutip/steady.py
 
1
#This file is part of QuTIP.
 
2
#
 
3
#    QuTIP is free software: you can redistribute it and/or modify
 
4
#    it under the terms of the GNU General Public License as published by
 
5
#    the Free Software Foundation, either version 3 of the License, or
 
6
#   (at your option) any later version.
 
7
#
 
8
#    QuTIP is distributed in the hope that it will be useful,
 
9
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
 
10
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
11
#    GNU General Public License for more details.
 
12
#
 
13
#    You should have received a copy of the GNU General Public License
 
14
#    along with QuTIP.  If not, see <http://www.gnu.org/licenses/>.
 
15
#
 
16
# Copyright (C) 2011, Paul D. Nation & Robert J. Johansson
 
17
#
 
18
###########################################################################
 
19
from scipy import prod, finfo
 
20
import scipy.sparse as sp
 
21
import scipy.linalg as la
 
22
from qobj import *
 
23
from istests import *
 
24
from scipy.sparse.linalg import spsolve
 
25
def steady(L):
 
26
        tol=1e-6
 
27
        maxiter=20
 
28
        eps=finfo(float).eps
 
29
        if (not isoper(L)) & (not issuper(L)):
 
30
                raise TypeError('Steady states can only be found for oeprators or superoperators.')
 
31
        rhoss=qobj()
 
32
        sflag=issuper(L)
 
33
        if sflag:
 
34
                rhoss.dims=L.dims[0]
 
35
                rhoss.shape=[prod(rhoss.dims[0]),prod(rhoss.dims[1])]
 
36
        else:
 
37
                rhoss.dims=[l.dims[0],1]
 
38
                rhoss.shape=[prod(rhoss.dims[0]),1]
 
39
        n=len(L.data.todense())
 
40
        L1=L.data+eps*la.norm(L.data.todense(),inf)*sp.eye(n,n)
 
41
        v=randn(n,1)
 
42
        it=0
 
43
        while (la.norm(L.data*v,inf)>tol) & (it<maxiter):
 
44
                v=spsolve(L1,v)
 
45
                v=v/la.norm(v,inf)
 
46
                it=it+1
 
47
        if it>=maxiter:
 
48
                raise ValueError('Failed to find steady state after ' + str(maxiter) +' iterations')
 
49
        rhoss.data=v
 
50
        #normalise according to type of problem
 
51
        if sflag:
 
52
                trow=eye(rhoss.shape[0],rhoss.shape[0])
 
53
                trow=reshape(trow,(1,n))
 
54
                rhoss.data=rhoss.data/sum(dot(trow,rhoss.data))
 
55
        else:
 
56
                rhoss.data=rhoss.data/la.norm(rhoss.data)
 
57
        rhoss.data=reshape(rhoss.data,(rhoss.shape[0],rhoss.shape[1])).T
 
58
        out=qobj(rhoss.data)
 
59
        out.dims=rhoss.dims
 
60
        out.shape=rhoss.shape
 
61
        return out
 
62
        
 
63
        
 
64