'''Basic Table, Matrix and Vector functions for Python 2.0 (von mir angepasst)
License: Public Domain Author: Raymond Hettinger email: python@rcn.com
Updates and documentation: http://users.rcn.com/python/download/python.htm
Revision In Use: 'File %n, Ver %v, Date %f' '''
import sys
from py2or3 import UserList
if sys.version_info[ 0 ] == 3: # python3
from py2or3 import reduce
Version = 'File MATFUNC.PY, Ver 183, Date 12-Dec-2002,14:33:42'
#==============================================================================
# Anpassung an Python2.0
# assert-statements mit float-vergleichen auskommentiert
# class : list,complex,Gleichung22,Gleichung3,Tensor
# function: real,imag
# Author: C.Huettel
# email : christian.huettel@volkswagen.de
# Datum : 23.2.05
#==============================================================================
import types
class list(UserList):
def __init__(self,data):
UserList.__init__(self,data)
#==============================================================================
def test( text, soll, ist, tol=1.e-6, fout=1 ):
if hasattr( soll,'__len__'):
ok = 1
for i in range( len( soll ) ):
ok = test( '', soll[ i ], ist[ i ], tol )
if not ok:
break
else:
ok = abs( soll - ist ) < tol
if text and fout:
print( text.ljust( 8, ' ' ), ok )
if not ok:
print( 'ist ', ist )
print( 'soll', soll )
return ok
#==============================================================================
def cumtest( cum, text, soll, ist, tol=1.e-6, fout=1 ):
( ok, n ) = cum
ok1 = test( text, soll, ist, tol, fout )
cum[:] = [ ok + ok1, n + 1 ]
return ok
#==============================================================================
def rk4( fkt, t, y0, fktpar=[], tol=1.e-2 ):
"""Integriert y'(t) = fkt( t, y(t)), y(t[0])=y0
mit Runge-Kutta und Schrittweitensteuerung
Fehlerschaetzer:
Bei Linearitaet: y1,2 = 0.5*( y0 + yi )
Fehler: | y1 + y2 - y0 - yi |/|yi + y0|
"""
o6 = 0.16666666666666666 # 1./6.
o3 = 0.33333333333333331 # 1./3.
t0 = t[ 0 ]
t1 = t[ 1 ]
eps = tol*( t[ -1 ] - t0 )/len( t )
ans = [ y0 ]
tout= [ t0 ]
h = t1 - t0
err = -1.
ok = 0
i = 1
while 1:
k0 = fkt( t0, y0, *fktpar )
nr = 0
while 1:
nr += 1
if ( t0 > t[ 0 ] and nr > 5 ) or nr > 20:
raise StopIteration(
'maximale Verfeinerung %d erreicht'%nr )
if t0 + h > t1 - eps:
h0 = h
h = t1 - t0
ok = 1
h2 = 0.5*h
y1 = y0 + h2*k0
k1 = fkt( t0 + h2, y1, *fktpar )
y2 = y0 + h2*k1
k2 = fkt( t0 + h2, y2, *fktpar )
k3 = fkt( t0 + h, y0 + h*k2, *fktpar )
yi = y0 + h*( o6*( k0 + k3 ) + o3*( k1 + k2 ) )
e1 = ( y1 + y2 - yi - y0 ).norm( )
n1 = ( yi + y0 ).norm( )
n2 = yi.norm( )
if n1 > 0.:
fe = e1/( tol*n1 )
elif n2 > 0.:
fe = e1/( tol*n2 )
else:
fe = e1/tol
if fe >= 1.:
if ok:
h = 0.5*h0
else:
h *= 0.5
ok = 0
else:
err = max( err, e1 )
y0 = yi
break
t0 += h
tout.append( t0 )
ans.append( y0 )
if ok:
ok = 0
i += 1
if i >= len( t ):
break
t1 = t[ i ]
h = min( t1 - t0, h0 )
elif fe < 0.7:
h *= 1.5
return ( tout, ans, err )
#==============================================================================
def rk4h( fkt, y0, zeit, fktpar=[] ):
"""Integriert y'(t) = fkt( t, y(t)), y(t[0])=y0
mit Runge-Kutta und ohne Schrittweitensteuerung"""
o6 = 0.16666666666666666 # 1./6.
o3 = 0.33333333333333331 # 1./3.
tol = 1.e-2
t0 = zeit[ 0 ]
t1 = zeit[ 1 ]
ans = [ y0 ]
h = t1 - t0
ok = 0
i = 1
#----------------------------------------------------------------------
def add( a, b ):
if hasattr( a, '__len__' ):
a = Vec( a )
b = Vec( b )
return a + b
#----------------------------------------------------------------------
def mul( a, b ):
if hasattr( b, '__len__' ):
b = Vec( b )
return a*b
#----------------------------------------------------------------------
for i in range( len( zeit ) - 1 ):
t = zeit[ i ]
h = zeit[ i + 1 ] - t
h2 = 0.5*h
k1 = fkt( y0, t, *fktpar )
k2 = fkt( add( y0, mul( h2, k1 ) ), t + h2, *fktpar )
k3 = fkt( add( y0, mul( h2, k2 ) ), t + h2, *fktpar )
k4 = fkt( add( y0, mul( h, k3 ) ), t + h, *fktpar )
y0 = add( y0, mul( h, add( mul( o6, add( k1, k4 ) ),
mul( o3, add( k2, k3 ) ) ) ) )
ans.append( y0 )
return ans
#==============================================================================
def real(x):
r=x
if isinstance(x,complex):
r=x.real
elif ( isinstance(x,Table) or type(x)==types.ListType
or type(x)==types.TupleType ):
r=map(lambda z: real(z),x)
if type(x)==types.TupleType:
r=tuple(r)
return r
#==============================================================================
def imag(x):
r=0.
if isinstance(x,complex):
r=x.imag
elif ( isinstance(x,Table) or type(x)==types.ListType
or type(x)==types.TupleType ):
r=map(lambda z: imag(z),x)
if type(x)==types.TupleType:
r=tuple(r)
return r
#==============================================================================
import operator, math, random
NPRE, NPOST = 0, 0 # Disables pre and post condition checks
#==============================================================================
def iszero(z): return abs(z) < .000001
def getreal(z):
try:
return z.real
except AttributeError:
return z
def getimag(z):
try:
return z.imag
except AttributeError:
return 0
def getconj(z):
try:
return z.conjugate()
except AttributeError:
return z
separator = [ '', '\t', '\n', '\n----------\n', '\n===========\n' ]
#==============================================================================
class Table(list):
dim = 1
concat = list.__add__ # A substitute for the overridden __add__ method
def __init__( self, elems ):
list.__init__( self, elems )
if len(self):
if hasattr(self[0], 'dim'):
self.dim = self[0].dim + 1
else:
self.dim = 1
def __getslice__( self, i, j ):
return self.__class__( list.__getslice__(self,i,j) )
def __str__( self ):
return separator[self.dim].join( map(str, self) )
def map( self, op, rhs=None ):
'''Apply a unary operator to every element in the matrix or a binary operator to corresponding
elements in two arrays. If the dimensions are different, broadcast the smaller dimension over
the larger (i.e. match a scalar to every element in a vector or a vector to a matrix).'''
if rhs is None: # Unary case
return self.dim==1 and self.__class__( map(op, self) ) or self.__class__( [elem.map(op) for elem in self] )
elif not hasattr(rhs,'dim'): # List / Scalar op
return self.__class__( [op(e,rhs) for e in self] )
elif self.dim == rhs.dim: # Same level Vec / Vec or Matrix / Matrix
assert NPRE or len(self) == len(rhs), 'Table operation requires len sizes to agree'
return self.__class__( map(op, self, rhs) )
elif self.dim < rhs.dim: # Vec / Matrix
return self.__class__( [op(self,e) for e in rhs] )
return self.__class__( [op(e,rhs) for e in self] ) # Matrix / Vec
def __mul__( self, rhs ):
return self.map( operator.mul, rhs )
def __div__( self, rhs ):
return self.map( operator.div, rhs )
def __sub__( self, rhs ): return self.map( operator.sub, rhs )
def __add__( self, rhs ): return self.map( operator.add, rhs )
def __rmul__( self, lhs ):
return self*lhs
def __rdiv__( self, lhs ): return self*(1.0/lhs)
def __rsub__( self, lhs ): return -(self-lhs)
def __radd__( self, lhs ): return self+lhs
def __abs__( self ): return self.map( abs )
def __neg__( self ): return self.map( operator.neg )
def conjugate( self ): return self.map( getconj )
def real( self ): return self.map( getreal )
def imag( self ): return self.map( getimag )
def flatten( self ):
if self.dim == 1: return self
return reduce( lambda cum, e: e.flatten().concat(cum), self, [] )
def prod( self ): return reduce(operator.mul, self.flatten(), 1.0)
def sum( self ): return reduce(operator.add, self.flatten(), 0.0)
def exists( self, predicate ):
for elem in self.flatten():
if predicate(elem):
return 1
return 0
def forall( self, predicate ):
for elem in self.flatten():
if not predicate(elem):
return 0
return 1
def __eq__( self, rhs ): return (self - rhs).forall( iszero )
#==============================================================================
class Vec(Table):
def __init__( self, elems ):
Table.__init__( self, elems )
#--------------------------------------------------------------------------
def dot( self, otherVec ): return reduce(
operator.add, map(operator.mul, self, otherVec), 0.0)
#--------------------------------------------------------------------------
def mmul( self, other ):
'Vec.Matrix'
return other.tr().mmul( self )
#--------------------------------------------------------------------------
def norm( self ): return math.sqrt(abs( self.dot(self) ))
#--------------------------------------------------------------------------
def dist( a, b ): return ( a - b ).norm( )
#--------------------------------------------------------------------------
def normalize( self ):
x = self.norm()
return Vec( [ y/x for y in self ] )
#--------------------------------------------------------------------------
def angle( self, other ):
aa = self.dot( self )
bb = other.dot( other )
ab = self.dot( other )
return math.acos( min( 1., max( -1., ab/math.sqrt( aa*bb ) ) ) )
#--------------------------------------------------------------------------
def outer( self, otherVec ): return Mat([otherVec*x for x in self])
#--------------------------------------------------------------------------
def diag( u ):
( x, y, z ) = u
return Tensor( [ [ x, 0., 0. ],
[ 0., y, 0. ],
[ 0., 0., z ] ] )
#--------------------------------------------------------------------------
def axialerTensor( u ):
( x, y, z ) = u
return Tensor( [ [ 0., -z, y ],
[ z, 0., -x ],
[ -y, x, 0. ] ] )
#--------------------------------------------------------------------------
def drehTensor( self ):
'''gibt den orthogonalen Tensor, der um diese Achse mit diesem Betrag
dreht. Ist nur in 3 Dimensionen definiert.'''
t = Tensor( eye( 3 ) )
alp = self.norm( )
#print 'alp2 =', alp*180./math.pi
if alp > 0.:
ud = self.normalize( ).axialerTensor( )
t = ( t + math.sin( alp )*ud
+ ( 1. - math.cos( alp ) )*ud.mmul( ud ) )
return t
#--------------------------------------------------------------------------
def cross( self, otherVec ):
'Compute a Vector or Cross Product with another vector'
assert len(self) == len(otherVec) == 3, 'Cross product only defined for 3-D vectors'
u, v = self, otherVec
if v.dim == 2: # v x T
E = eye( 3,3 )
return Tensor( sum( [ u.cross( E[ j ] ).outer( v[ j ] )
for j in ( 0,1,2 ) ], 0*E ) )
else:
return Vec([ u[1]*v[2]-u[2]*v[1],
u[2]*v[0]-u[0]*v[2],
u[0]*v[1]-u[1]*v[0] ])
def house( self, index ):
'Compute a Householder vector which zeroes all but the index element after a reflection'
v = Vec( Table([0]*index).concat(self[index:]) ).normalize()
t = v[index]
sigma = 1.0 - t**2
if sigma != 0.0:
t = v[index] = t<=0 and t-1.0 or -sigma / (t + 1.0)
v /= t
return v, 2.0 * t**2 / (sigma + t**2)
def polyval( self, x ):
'Vec([6,3,4]).polyval(5) evaluates to 6*x**2 + 3*x + 4 at x=5'
return reduce( lambda cum,c: cum*x+c, self, 0.0 )
def ratval( self, x ):
'Vec([10,20,30,40,50]).ratfit(5) evaluates to (10*x**2 + 20*x + 30) / (40*x**2 + 50*x + 1) at x=5.'
degree = len(self) / 2
num, den = self[:degree+1], self[degree+1:] + [1]
return num.polyval(x) / den.polyval(x)
#==============================================================================
class Matrix(Table):
__slots__ = ['size', 'rows', 'cols']
def __init__( self, elems ):
'Form a matrix from a list of lists or a list of Vecs'
Table.__init__( self, hasattr(elems[0], 'dot') and elems or map(Vec,map(list,elems)) )
self.size = self.rows, self.cols = len(elems), len(elems[0])
def tr( self ):
'Tranpose elements so that Transposed[i][j] = Original[j][i]'
return self.__class__(zip(*self))
def star( self ):
'Return the Hermetian adjoint so that Star[i][j] = Original[j][i].conjugate()'
return self.tr().conjugate()
def diag( self ):
'Return a vector composed of elements on the matrix diagonal'
return Vec( [self[i][i] for i in range(min(self.size))] )
def trace( self ): return self.diag().sum()
def spur( self ): return self.trace()
def invar2( self ): return 0.5*( self.spur()**2 - self.mmul( self ).spur())
def dot( self, other ):
'Matrix multiply by another matrix or a column vector '
return self.mmul( other.tr() ).spur()
def mmul( self, other ):
'Matrix multiply by another matrix or a column vector '
if other.dim==2: return self.__class__( map(self.mmul, other.tr()) ).tr()
assert NPRE or self.cols == len(other)
return Vec( map(other.dot, self) )
def augment( self, otherMat ):
'Make a new matrix with the two original matrices laid side by side'
assert self.rows == otherMat.rows, 'Size mismatch: %s * %s' % ("self.size", "otherMat.size")
return Mat( map(Table.concat, self, otherMat) )
def qr( self, ROnly=0 ):
'QR decomposition using Householder reflections: Q*R==self, Q.tr()*Q==I(n), R upper triangular'
R = self
m, n = R.size
for i in range(min(m,n)):
v, beta = R.tr()[i].house(i)
R -= v.outer( R.tr().mmul(v)*beta )
for i in range(1,min(n,m)): R[i][:i] = [0] * i
R = Mat(R[:n])
if ROnly: return R
Q = R.tr().solve(self.tr()).tr() # Rt Qt = At nn nm = nm
## self.qr = lambda r=0, c="self": not r and c=="self" and (Q,R) or Matrix.qr(self,r) #Cache result
## self.qr = (Q,R)
## assert NPOST or m>=n and Q.size==(m,n) and isinstance(R,UpperTri) or m<n and Q.size==(m,m) and R.size==(m,n)
# assert NPOST or Q.mmul(R)==self and Q.tr().mmul(Q)==eye(min(m,n))
return Q, R
def _solve( self, b ):
( L, U ) = self.lu( )
return U._solve( L._solve( b ) )
def _solveQR( self, b ):
'''General matrices (incuding) are solved using the QR composition.
For inconsistent cases, returns the least squares solution'''
#print 'Matrix._solve'
Q, R = self.qr()
return R.solve( Q.tr().mmul(b) )
def solve( self, b ):
'Divide matrix into a column vector or matrix'
return self._solve( b )
def solveOpt( self, b ):
'Divide matrix into a column vector or matrix and iterate to improve the solution'
if b.dim==2: return Mat( map(self.solve, b.tr()) ).tr()
assert NPRE or self.rows == len(b), 'Matrix row count %d must match vector length %d' % (self.rows, len(b))
x = self._solve( b )
diff = b - self.mmul(x)
maxdiff = diff.dot(diff)
for i in range(10):
xnew = x + self._solve( diff )
diffnew = b - self.mmul(xnew)
maxdiffnew = diffnew.dot(diffnew)
if maxdiffnew >= maxdiff: break
x, diff, maxdiff = xnew, diffnew, maxdiffnew
#print >> sys.stderr, i+1, maxdiff
# assert NPOST or self.rows!=self.cols or self.mmul(x) == b
return x
def rank( self ): return Vec([ not row.forall(iszero)
for row in self.qr(ROnly=1) ]).sum()
def norm( self ): return math.sqrt( sum( map( sum, self*self ) ) )
#==============================================================================
class Square(Matrix):
def lu( self ):
'Factor a square matrix into lower and upper triangular form such that L.mmul(U)==A'
n = self.rows
L, U = eye(n), Mat(self[:])
for i in range(n):
for j in range(i+1,U.rows):
# assert U[i][i] != 0.0, 'LU requires non-zero elements on the diagonal'
L[j][i] = m = 1.0 * U[j][i] / U[i][i]
U[j] -= U[i] * m
# assert NPOST or isinstance(L,LowerTri) and isinstance(U,UpperTri) and L*U==self
return LowerTri(L), UpperTri(U)
def __pow__( self, exp ):
'Raise a square matrix to an integer power (i.e. A**3 is the same as A.mmul(A.mmul(A))'
# assert NPRE or exp==int(exp) and exp>0, 'Matrix powers only defined for positive integers not %s' % exp
if exp == 1: return self
if exp&1: return self.mmul(self ** (exp-1))
sqrme = self ** (exp/2)
return sqrme.mmul(sqrme)
#--------------------------------------------------------------------------
def i( self ):
return Vec( [ self[ 1 ][ 2 ] - self[ 2 ][ 1 ],
self[ 2 ][ 0 ] - self[ 0 ][ 2 ],
self[ 0 ][ 1 ] - self[ 1 ][ 0 ] ] )
#--------------------------------------------------------------------------
def det( self ):
try:
return reduce( lambda p,x: p*x, self.lu( )[ 1 ].diag( ), 1. )
except ZeroDivisionError:
return 0.
#--------------------------------------------------------------------------
def detQR( self ):
return self.qr( ROnly=1 ).det()
#--------------------------------------------------------------------------
def inverse( self ): return self.solve( eye(self.rows) )
#--------------------------------------------------------------------------
def hessenberg( self ):
'''Householder reduction to Hessenberg Form (zeroes below the diagonal)
while keeping the same eigenvalues as self.'''
for i in range(self.cols-2):
v, beta = self.tr()[i].house(i+1)
self -= v.outer( self.tr().mmul(v)*beta )
self -= self.mmul(v).outer(v*beta)
return self
def hessenberg2( self ):
'''Householder reduction to Hessenberg Form (zeroes below the diagonal)
while keeping the same eigenvalues as self.'''
ans = self
for i in range( len( self ) - 2 ):
st = ans.tr( )
( v, beta ) = st[ i ].house( i + 1 )
ans = ans - v.outer( st.mmul( v )*beta )
ans = ans - ans.mmul( v ).outer( v*beta )
return ans
#--------------------------------------------------------------------------
def eigs( self ):
'Estimate principal eigenvalues using the QR with shifts method'
origTrace, origDet = self.trace(), self.det()
self = self.hessenberg()
eigvals = Vec([])
for i in range(self.rows-1,0,-1):
while not self[i][:i].forall(iszero):
shift = eye(i+1) * self[i][i]
q, r = (self - shift).qr()
self = r.mmul(q) + shift
eigvals.append( self[i][i] )
self = Mat( [self[r][:i] for r in range(i)] )
eigvals.append( self[0][0] )
# assert NPOST or iszero( (abs(origDet) - abs(eigvals.prod())) / 1000.0 )
# assert NPOST or iszero( origTrace - eigvals.sum() )
return Vec(eigvals)
#--------------------------------------------------------------------------
#==============================================================================
class Gleichung22(Square):
"definiert ein 2 x 2 Gleichungssystem"
#--------------------------------------------------------------------------
def __init__(self,elems):
Square.__init__(self,elems)
if self.rows!=2 or self.cols!=2:
raise IndexError( 'Gleichung22 mit rows, cols!=2 aufgerufen' )
#--------------------------------------------------------------------------
def det(self):
return self[0][0]*self[1][1]-self[0][1]*self[1][0]
#--------------------------------------------------------------------------
def solve(self,rhs=[]):
"solve(rhs=[])->list loest self*x=rhs, gibt x zurueck"
det=self.det()
x=(self[1][1]*rhs[0]-self[0][1]*rhs[1])/det
y=(-self[1][0]*rhs[0]+self[0][0]*rhs[1])/det
return [x,y]
#--------------------------------------------------------------------------
#==============================================================================
class Gleichung3:
"eine Gleichung 3.Grades"
o3=1./3.
d120=2.0943951023931953 # =2*pi/3
d240=d120+d120
#--------------------------------------------------------------------------
def __init__( self, elems ):
self.elems = elems
if len( elems )!=4:
raise( IndexError,'Gleichung3 mit n!=4 Koeffizienten aufgerufen' )
#--------------------------------------------------------------------------
def wert( self, x=0. ):
return self.elems[0] + x*( self.elems[1] + x*( self.elems[2]
+ x*self.elems[3] ) )
#--------------------------------------------------------------------------
def solve(self):
"gibt die drei complexen Nullstellen dieses Polynoms zurueck"
a0=self.elems[0]/self.elems[3]
a1=self.elems[1]/self.elems[3]
a2=self.elems[2]/self.elems[3]
y=-a2/3.
p=a1/3.-a2*a2/9.
q=0.5*a0-a2*a1/6.+a2*a2*a2/27.
z=q*q+p*p*p
if z>0.:
r=-q+math.sqrt(z)
if r<0.:
u=-(-r)**self.o3
else:
u=r**self.o3
r=-q-math.sqrt(z)
if r<0.:
v=-(-r)**self.o3
else:
v=r**self.o3
x1=-0.5*(u+v)+y
dif=math.sqrt(0.75)*(u-v)
x0=[complex(u+v+y,0.),complex(x1,dif), complex(x1,-dif)]
else:
# winkel phi = z
r=math.sqrt(-q*q/(p**3))
if q*r>0.:
r=-r
z=math.acos(r)
# radius r
r=2.*math.sqrt(-p)
x0=[complex(r*math.cos(z/3.)+y,0.),
complex(r*math.cos(z/3.+self.d120)+y,0.),
complex(r*math.cos(z/3.+self.d240)+y,0.)]
return x0
#--------------------------------------------------------------------------
#==============================================================================
class Gleichung4:
"eine Gleichung 4.Grades"
#--------------------------------------------------------------------------
def __init__( self, elems ):
self.elems = elems
if len( elems ) != 5:
raise IndexError( 'Gleichung4 mit n!=5 Koeffizienten aufgerufen' )
#--------------------------------------------------------------------------
def wert( self, x=0.):
a = self.elems
return a[0] + x*( a[1] + x*( a[2] + x*( a[3] + x*a[4] ) ) )
#--------------------------------------------------------------------------
def getReal( self ):
"gibt die reellen Nullstellen dieses Polynoms zurueck"
ans = [ ]
for z in self.solve( ):
if abs( z.imag ) < 1.e-6*abs( z.real ):
ans.append( z.real )
return ans
#--------------------------------------------------------------------------
def solve( self ):
"gibt die vier complexen Nullstellen dieses Polynoms zurueck"
aa = self.elems[ 4 ]
b = self.elems[ 2 ]/aa
c = self.elems[ 1 ]/aa
d = self.elems[ 0 ]/aa
aa = self.elems[ 3 ]/aa
y = -0.25*aa + 0j
a2 = aa*aa
p = b - 0.375*a2
q = c - ( 0.5*b - 0.125*a2 )*aa
r = d - 0.25*aa*c + ( b - 0.1875*a2 )*0.0625*a2
x0 = 4*[ 0. ]
if q == 0.: # biquadratische Gleichung
# x^4 + p*x^2 + r = 0 -> x^2 = -p/2 +- sqrt( p^2/4 - r )
w = complex( 0.25*p**2 - r, 0. )**0.5
x0[ 0 ] = ( -0.5*p + w )**0.5
x0[ 1 ] = -x0[ 0 ]
x0[ 2 ] = ( -0.5*p - w )**0.5
x0[ 3 ] = -x0[ 2 ]
else:
# Loesung der kubischen Resolvente
# ( x^2 - P )^2 = ( Qx + R )^2
# x^4 - ( 2P + Q^2 ) x^2 - 2QR x + P^2 - R^2 = 0
# p = - 2P - Q^2
# q = - 2QR
# r = P^2 - R^2
proots = Gleichung3( [ 0.125*( q*q - 4.*p*r ), -r, 0.5*p, 1. ]
).solve( )
P = proots[ 0 ].real
if P*P < r:
P = proots[ 1 ].real
if P*P < r:
P = proots[ 2 ].real
if P*P < r:
raise ArithmeticException(
"Gl4.solve(): keine Loesung" )
R = math.sqrt( P*P - r )
Q = -0.5*q/R
z = 0.25*Q*Q + P + R
# x^2 - P = +-( Qx + R )
# x^2 -+ Qx - P -+ R = 0
# x = +- Q/2 +- sqrt( Q^2/4 + P +- R )
if z < 0.:
x0[ 0 ] = complex( 0.5*Q, math.sqrt( -z ) )
x0[ 1 ] = x0[ 0 ].conjugate( )
else:
x0[ 0 ] = complex( 0.5*Q + math.sqrt( z ), 0. )
x0[ 1 ] = complex( 0.5*Q - math.sqrt( z ), 0. )
z = 0.25*Q*Q + P - R
if z < 0.:
x0[ 2 ] = complex( -0.5*Q, math.sqrt( -z ) )
x0[ 3 ] = x0[ 2 ].conjugate( )
else:
x0[ 2 ] = complex( -0.5*Q + math.sqrt( z ), 0. )
x0[ 3 ] = complex( -0.5*Q - math.sqrt( z ), 0. )
x0[ 0 ] += y
x0[ 1 ] += y
x0[ 2 ] += y
x0[ 3 ] += y
return x0
#--------------------------------------------------------------------------
#==============================================================================
class Tensor(Square):
"eine 3x3 Matrix"
#--------------------------------------------------------------------------
def __init__(self,elems):
Square.__init__(self,elems)
if self.rows!=3 or self.cols!=3:
raise IndexError("Tensor with more than 3x3 Elements called" )
#--------------------------------------------------------------------------
def det(self):
"det() -> float: Gibt die Determinante dieses Tensors zurueck"
return ( self[0][0]*( self[1][1]*self[2][2] - self[1][2]*self[2][1] )
+ self[0][1]*( self[1][2]*self[2][0] - self[1][0]*self[2][2] )
+ self[0][2]*( self[1][0]*self[2][1] - self[1][1]*self[2][0] ))
#--------------------------------------------------------------------------
def cof( A ):
'Kofaktor'
return 0.5*A.touter( A )
#--------------------------------------------------------------------------
def adj( A ):
'Adjunkte'
return A.cof( ).tr( )
#--------------------------------------------------------------------------
def sym( self ):
'gibt den symmetrischen Anteil dieses Tensors'
return 0.5*( self + self.tr( ) )
#--------------------------------------------------------------------------
def skw( self ):
'gibt den schief-symmetrischen Anteil dieses Tensors'
return 0.5*( self - self.tr( ) )
#--------------------------------------------------------------------------
def sph( self ):
'gibt den Kugel-Anteil dieses Tensors'
return Eye( self.spur( )/len( self ), dtype=self.dtype )
#--------------------------------------------------------------------------
def dev( self ):
'gibt den Deviator dieses Tensors'
return self - self.sph( )
#--------------------------------------------------------------------------
def cross( A, v ):
'''gibt das Kreuzprodukt zwischen Tensor A und Vektor v
( ei o ai ) x v = ei o ( ai x v )'''
( a1, a2, a3 ) = A # Zeilenvektoren
return Tensor( [ a1.cross( v ), a2.cross( v ), a3.cross( v ) ] )
#--------------------------------------------------------------------------
def touter( A, B ):
'gibt das aeussere Tensorprodukt dieses Tensors mit B'
E = Tensor( eye( 3,3 ) )
return ( ( A.spur( )*B.spur( ) - A.mmul( B ).spur( ) )*E
+ ( B.mmul( A ) + A.mmul( B ) - A.spur( )*B - B.spur( )*A ).tr( ) )
#--------------------------------------------------------------------------
def eigenVektorZu(self,x=0.):
"Berechnet den Eigenvektor zu einem gegebenen Eigenwert x"
t11=complex(real=self[0][0]); t12=complex(real=self[0][1])
t13=complex(real=self[0][2]); t21=complex(real=self[1][0])
t22=complex(real=self[1][1]); t23=complex(real=self[1][2])
t31=complex(real=self[2][0]); t32=complex(real=self[2][1])
t33=complex(real=self[2][2])
# bestes gleichungssystem auswaehlen
g1 = Gleichung22([[t22 - x, t23], [t32 , t33 - x]])
a1 = abs(g1.det())
g2 = Gleichung22([[t11 - x, t13], [t31 , t33 - x]])
a2 = abs(g2.det())
g3 = Gleichung22([[t11 - x, t12], [t21 , t22 - x]])
a3 = abs(g3.det())
# eigenvektor berechnen
v1 = complex(real=1.)
v2 = complex(real=1.)
v3 = complex(real=1.)
#print 'matfunc2'
if a1<a2:
if a2<a3:
# print "a3 am groessten"
[v1, v2] = g3.solve(rhs = [-t13, -t23])
else:
# print "a2 am groessten"
[v1, v3] = g2.solve(rhs = [-t12, -t32])
elif a1<a3: # a2 <= a1
# print " a3 am groessten"
[v1,v2] = g3.solve(rhs = [-t13, -t23])
else: # a3 <= a1 und a2 <= a1
# print " a1 am groessten"
[v2,v3] = g1.solve(rhs = [-t21, -t31])
# normieren
a=math.sqrt((v1*v1.conjugate() + v2*v2.conjugate() + v3*v3.conjugate()).real)
return Vec([v1/a,v2/a,v3/a])
#--------------------------------------------------------------------------
def eigenVectors(self,lam=[]):
return (self.eigenVektorZu(complex(lam[0])),
self.eigenVektorZu(complex(lam[1])),
self.eigenVektorZu(complex(lam[2])))
#--------------------------------------------------------------------------
def invar2(self):
"invar2() -> float: Gibt die zweite Invariante zurueck"
return 0.5*(self.trace()**2-self.mmul(self).trace())
#--------------------------------------------------------------------------
def eigenValues(self):
"berechnet die complexen Eigenwerte dieses Tensors"
return Gleichung3([self.det(),-self.invar2(),self.trace(),-1.]).solve()
#--------------------------------------------------------------------------
def eig(self):
m = len(self)
if m!=3:
raise IndexError('EigenSystem nur fuer 3x3-Matrizen definiert')
lam=self.eigenValues()
return (lam,self.eigenVectors(lam))
#--------------------------------------------------------------------------
def svd(self):
"""( U, S, V ) = T.svd() erstellt Tensoren U, S und V mit
1) T==U.tProd(S).tProd(V.tr())
2) U.tProd(U.tr())==1
3) V.tProd(V.tr())==1
4) S diagonal Gestalt"""
# rechts eigen system
x=Tensor(self.tr().mmul(self)).eigenSystem()
v=Tensor([x[1][0].real(),x[1][1].real(),x[1][2].real()]).tr()
# links eigen system
x=Tensor(self.mmul(self.tr())).eigenSystem()
u=Tensor([x[1][0].real(),x[1][1].real(),x[1][2].real()])
# eigenwerte
x=u.mmul(self).mmul(v)
s=Tensor([[x[0][0],0.,0.],[0.,x[1][1],0.],[0.,0.,x[2][2]]])
return (u.tr(),s,v)
#--------------------------------------------------------------------------
def apply(self,f=None):
"apply(f=None)->Tensor\n\
wendet die Funktion f auf diesen Tensor an"
(u,a,v)=self.svd()
s=Tensor([[f(a[0][0]),0.,0.],[0.,f(a[1][1]),0.],
[0.,0.,f(a[2][2])]])
return Tensor(u.mmul(s).mmul(v.tr()))
#--------------------------------------------------------------------------
def polarDecomposition( self ):
"t.polareZerlegung()->(u,v,r)\n\
berechnet Tensoren u, v und r so, dass\n\
1) t==r.mmul(u)==v.mmul(r) \n\
2) r orthogonal\n\
3) u und v symmetrisch positiv definit"
# u
( q, s, p )=Tensor( self.tr( ).mmul( self ) ).svd( )
s = s.map( math.sqrt )
U = Tensor( q.mmul( s ).mmul( p ) )
( q, s, p )=Tensor( self.mmul( self.tr( ) ) ).svd( )
V = Tensor( q.mmul( s ).mmul( p ) )
R = Tensor( self.mmul( U.inverse( ) ) )
return ( U, V, R )
#--------------------------------------------------------------------------
#==============================================================================
class Triangular(Square):
def eigs( self ): return self.diag()
def det( self ): return self.diag().prod()
#==============================================================================
class UpperTri(Triangular):
def _solve( self, b ):
'Solve an upper triangular matrix using backward substitution'
#print 'UpperTri._solve'
x = Vec([])
for i in range(self.rows-1, -1, -1):
# assert NPRE or self[i][i], 'Backsub requires non-zero elements on the diagonal'
x.insert(0, (b[i] - x.dot(self[i][i+1:])) / self[i][i] )
return x
#==============================================================================
class LowerTri(Triangular):
def _solve( self, b ):
'Solve a lower triangular matrix using forward substitution'
#print 'LowerTri._solve'
x = Vec([])
for i in range(self.rows):
# assert NPRE or self[i][i], 'Forward sub requires non-zero elements on the diagonal'
x.append( (b[i] - x.dot(self[i][:i])) / self[i][i] )
return x
#==============================================================================
def Mat( elems ):
'Factory function to create a new matrix.'
m, n = len(elems), len(elems[0])
if m != n: return Matrix(elems)
# m == n
if n <= 1: return Square(elems)
for i in range(1, len(elems)):
if not iszero( max(map(abs, elems[i][:i])) ):
break
else: return UpperTri(elems)
for i in range(0, len(elems)-1):
if not iszero( max(map(abs, elems[i][i+1:])) ):
return Square(elems)
return LowerTri(elems)
#------------------------------------------------------------------------------
def funToVec( tgtfun, low=-1, high=1, steps=40, EqualSpacing=0 ):
'''Compute x,y points from evaluating a target function over an interval (low to high)
at evenly spaces points or with Chebyshev abscissa spacing (default) '''
if EqualSpacing:
h = (0.0+high-low)/steps
xvec = [low+h/2.0+h*i for i in range(steps)]
else:
scale, base = (0.0+high-low)/2.0, (0.0+high+low)/2.0
xvec = [base+scale*math.cos(((2*steps-1-2*i)*math.pi)/(2*steps)) for i in range(steps)]
yvec = map(tgtfun, xvec)
return Mat( [xvec, yvec] )
#------------------------------------------------------------------------------
def funfit( arg, basisfuns ):
'Solves design matrix for approximating to basis functions'
(xvec, yvec) = arg
return Mat([ map(form,xvec) for form in basisfuns ]).tr().solve(Vec(yvec))
#------------------------------------------------------------------------------
def polyfit( arg, degree=2 ):
'Solves Vandermonde design matrix for approximating polynomial coefficients'
(xvec, yvec) = arg
return Mat([ [x**n for n in range(degree,-1,-1)] for x in xvec ]).solve(Vec(yvec))
#------------------------------------------------------------------------------
def ratfit( arg, degree=2 ):
'Solves design matrix for approximating rational polynomial coefficients (a*x**2 + b*x + c)/(d*x**2 + e*x + 1)'
(xvec, yvec) = arg
return Mat([[x**n for n in range(degree,-1,-1)]+[-y*x**n for n in range(degree,0,-1)] for x,y in zip(xvec,yvec)]).solve(Vec(yvec))
#------------------------------------------------------------------------------
def genmat(m, n, func):
if not n: n=m
return Mat([ [func(i,j) for i in range(n)] for j in range(m) ])
#------------------------------------------------------------------------------
def zeroes(m=1, n=None):
'Zero matrix with side length m-by-m or m-by-n.'
return genmat(m,n, lambda i,j: 0.)
#------------------------------------------------------------------------------
def eye(m=1, n=None):
'Identity matrix with side length m-by-m or m-by-n'
return genmat(m,n, lambda i,j: i==j and 1. or 0.)
#------------------------------------------------------------------------------
def hilb(m=1, n=None):
'Hilbert matrix with side length m-by-m or m-by-n. Elem[i][j]=1/(i+j+1)'
return genmat(m,n, lambda i,j: 1.0/(i+j+1.0))
#------------------------------------------------------------------------------
def rand(m=1, n=None):
'Random matrix with side length m-by-m or m-by-n'
return genmat(m,n, lambda i,j: random.random())
#------------------------------------------------------------------------------
