[SciPy-user] do complex numbers default to double precision?

Ryan Krauss ryanfedora at comcast.net
Thu Jun 30 10:17:06 CDT 2005


An HTML attachment was scrubbed...
URL: http://www.scipy.net/pipermail/scipy-user/attachments/20050630/e5ebafd7/attachment-0001.htm
-------------- next part --------------
from scipy import cos, cosh, sin, sinh, array, r_, c_
import scipy
from scipy.linalg import det
import copy
#from rwkmisc import null, norm2

class TMMElement:
    def __init__(self,elemtype,params):
        self.elemtype=self.parsetype(elemtype)
        self.params=params
        self.maxsize=8

    def parsetype(self, elemtype):
        elemout=-1
        if isinstance(elemtype,str):
            teststr=elemtype.lower()
            if teststr=='beam':
                elemout=1
            elif teststr.find('rigid')==0:
                elemout=2
        else:
            elemout=elemtype
        return elemout


    def GetMat(self, s):
        if self.elemtype==1:
            #this is a beam element
            L=self.params['L']
            mu=self.params['mu']
            EI1=self.params['EI1']
            EI2=self.params['EI2']
            mat1=bendmat(s,L,mu,EI1)
            mat2=bendmat(s,L,mu,EI2)
            zmat=scipy.zeros(scipy.shape(mat1))
            bigmat1=c_[mat1,zmat]
            bigmat2=c_[zmat,mat2]
            return r_[bigmat1,bigmat2]
            

def bendmat(s,L,mu,EI):
    beta=pow((-1*pow(s,2)*pow(L,4)*mu/EI),0.25)
    a=pow(L,2)/EI
    c0=(cosh(beta)+cos(beta))/2.0
    c1=(sinh(beta)+sin(beta))/(2.0*beta)
    c2=(cosh(beta)-cos(beta))/(2.0*pow(beta,2));
    c3=(sinh(beta)-sin(beta))/(2*pow(beta,3));
    matout=array([[c0, L*c1, a*c2, a*L*c3], [pow(beta,4)*c3/L, c0, a*c1/L, a*c2],[pow(beta,4)*c2/a, pow(beta,4)*L*c3/a, c0, L*c1],[pow(beta,4)*c1/(a*L), pow(beta,4)*c2/a, pow(beta,4)*c3/L, c0]])
    return matout


class TMMSystem:
    def __init__(self,elemlist,bcend=[2,3,6,7],bcbase=[0,1,4,5]):
        self.elemlist=elemlist#a list of TMMElement instances
        self.bcend=bcend#a list of the boundary condition indices that are 0 at the end of the chain (the free end in the case of a cantilever configuration) - the defaults are for a 1D cantilever where the states are [-w, phi, M, V] and M and V are 0 at the free end
        self.bcbase=bcbase#a list of the boundary configuration indices that are 0 as the base of the chain (the clamped beginning of a cantilever configuration) - the defaults are for a 1D cantilever where -w and phi are 0 at the base
        self.maxsize=elemlist[0].maxsize
        templist=range(self.maxsize)
        self.bccols=[]
        for ent in templist:
            if ent not in self.bcbase:
                self.bccols.append(ent)#in evaluating the eigvalues, the columns of the system U matrix that multiply the non-zero elements of the base state vector are the ones whose sub-determinant must go to 0



    def FindEig(self, guess,mytol=1e-15,maxfun=1e4,maxiter=1e3):
#        eigout=scipy.optimize.fmin(self.EigError,guess,xtol=mytol,ftol=mytol,maxfun=maxfun,maxiter=maxiter)
        if scipy.shape(guess)==():
            guess=scipy.array([scipy.real(guess),scipy.imag(guess)])
        eigout=scipy.optimize.fmin(self.EigError,guess,xtol=mytol,ftol=mytol,maxfun=maxfun,maxiter=maxiter)
#        eigout=scipy.optimize.fmin_bfgs(self.EigError,guess,gtol=1e-20)#,epsilon=1e-16)#,norm=2)
        return eigout


    def EigError(self, value):
        submat=self.FindSubmat(value)
        mydet=det(submat)
#        print('mydet='+'%1.8e'+' '+'%1.8e'+'j') % (scipy.real(mydet),scipy.imag(mydet))
#        mye=scipy.sqrt(norm2(mydet))
#        print('mye='+str(mye))
        mye=abs(mydet)
#        mye=pow(mye,0.1)
#        return norm2(mydet)
#        absval=scipy.absolute(mydet)
#        print('abs(mydet)='+str(absval))
        return mye

    def FindU(self,value):
        if scipy.shape(value):
            value=value[0]+value[1]*1j
        templist=copy.deepcopy(self.elemlist)
        curelem=templist.pop(0)
        U=curelem.GetMat(value)
        for curelem in templist:
            tempU=curelem.GetMat(value)
            U=tempU*U
        return U

    def FindModeShape(self, eigval):
        submat=self.FindSubmat(eigval)
        ns=null(submat)
        basevect=scipy.zeros((self.maxsize,1))#this will be the vector of boundary conditions at the base of the cantilever
        basevect=basevect*1j#force the vector to be cast into complex
#        print('ns[0,0]='+str(ns[0,0]))
#        print('type(ns[0,0])='+str(type(ns[0,0])))
#        print('shape(ns)='+str(scipy.shape(ns)))
        nsvect=ns[:,0]#*1e5#assume distinct eigenvalues/vectors and only 1D null space
        for curind,curns in zip(self.bccols,nsvect):
            basevect[curind,0]=curns#insert non-zero values into the boundary condition vector at the base of the beam
        endvect=scipy.matrixmultiply(self.FindU(eigval),basevect)#this is the vector of boundary conditions at the free end of the cantilever
        endvect=endvect/(scipy.linalg.norm(endvect))
        return endvect, basevect

    def FindSubmat(self,value):
#        print('s='+str(value))
#        print('shape(s)='+str(scipy.shape(value)))
#        if scipy.shape(value):
#            value=value[0]
        U=self.FindU(value)
        
        first=1
        for curri in self.bcend:
            currow=[]
            for curci in self.bccols:
#                print('curri='+str(curri))
#                print('curci='+str(curci))
#                print('shape(U)='+str(scipy.shape(U)))
                currow.append(U[curri,curci])
            if first:
                submat=array([currow])
                first=0
            else:
                submat=r_[submat,[currow]]
        return submat


def null(A, eps=1e-10):
    u, s, vh = scipy.linalg.svd(A)
    null_mask = (s <= eps)
    null_space = scipy.compress(null_mask, vh, axis=0) 
    return scipy.transpose(scipy.conj(null_space)) 


def norm2(x):
    """compute |x|^2 = x*conjugate(x)"""
    if iscomplexobj(x):
#        t1=time.time()
#        mat1=x.real**2 + x.imag**2
#        t2=time.time()
        mat2=multiply(x.real,x.real) + multiply(x.imag,x.imag)
#        t3=time.time()
#        print('---------------------------')
#        print('pow time='+str(t2-t1))
#        print('multiply time='+str(t3-t2))
#        print('pow time/multiply time='+str((t2-t1)/(t3-t2)))
#        print('shape(x)='+str(shape(x)))
#        print('x.typecode='+str(x.typecode()))
#        print('mat1.typecode='+str(mat1.typecode()))
#        print('mat2.typecode='+str(mat2.typecode()))
#        if len(shape(x))==1:
#            print('type(x[0])='+str(type(x[0])))
#            print('type(mat1[0])='+str(type(mat1[0])))
#            print('type(mat2[0])='+str(type(mat2[0])))
#        else:
#            print('type(x[0,0])='+str(type(x[0,0])))
#            print('type(mat1[0,0])='+str(type(mat1[0,0])))
#            print('type(mat2[0,0])='+str(type(mat2[0,0])))
        return mat2
    else:
        return multiply(x,x)
-------------- next part --------------
(dp0
S'EI1'
p1
F339137.17963813758
sS'm11'
p2
F0.026659562476630289
sS'EI2'
p3
F508705.76945720636
sS'EA'
p4
F145735941.95585597
sS'L'
p5
F4.6482000000000001
sS'mu'
p6
F5.7281487889477534
sS'GJ'
p7
F255774.40014479251
s.
-------------- next part --------------
import pickle
import TMM, scipy
f1=file("myparams.txt",'r')
myparams=pickle.load(f1)
f1.close()
mybeam=TMM.TMMElement('beam',myparams)#create a transfer matrix beam element
mysys=TMM.TMMSystem([mybeam])#create a system that contains only the beam - default boundary conditions are for a cantilever configuration
eig1=mysys.FindEig(40j)
subU=mysys.FindSubmat(eig1)#this method returns the system submatrix that should have a null space
print('det='+str(scipy.linalg.det(subU)))
ns=TMM.null(subU)#find the null space vector (4x1)
print('ns='+str(ns))
test=scipy.matrixmultiply(subU,ns)
print('test='+str(test))
freevect,basevect=mysys.FindModeShape(eig1)#this method combines the 4x1 null space vector with boundary conditions that are 0 to from the 8x1 basevect which is the vector of boundary conditions at the base of the beam.  freevect is the boundary condition vector at the free end of the beam.  It is found be freevect=U*basevect, where U is the transfer matrix for the system.
print('basevect='+str(basevect))
print('freevect='+str(freevect))
U=mysys.FindU(eig1)
freecheck=scipy.matrixmultiply(U,basevect)
print('freecheck='+str(freecheck))
freecheck=freecheck/(scipy.linalg.norm(freecheck))
print('freecheck (normalized)='+str(freecheck))


More information about the SciPy-user mailing list