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

Ryan Krauss ryanfedora at comcast.net
Thu Jun 30 09:32:58 CDT 2005


Thanks for your offer Nils.  Attached are two files.  TMM.py defines the 
transfer matrix classes.  initial_TMM_scipy_v2 calls these classes.  If 
you run initial_TMM_scipy_v2 it should work.  I run it from Ipython.  It 
generates a plot of the determinant to tell me where to start the 
optimization.  Because there is no damping in this model yet, the 
eigenvalue should actually be purely imaginary in this case, but I need 
the algorithm to be able to handle damping and therefore complex values.

I will look at the code and see if I can get you an easier to follow 
example if you would like.  But feel free to start with this.

Thanks,

Ryan


Nils Wagner wrote:

> Ryan Krauss wrote:
>
>> Thanks to Nils and Robert for their quick responses.  This is 
>> definitely one thing to love about SciPy.  I posted my null space 
>> question yesterday and Robert responded in 20 minutes and Nils in 
>> 30.  I posted my optimization question this morning and Nils reponded 
>> in 8 minutes!  It is almost like chatting with tech support.  You 
>> make SciPy a pleasure to use!  Thanks.
>>
>> Coincidentally, these two problems are linked and I don't know if my 
>> numerical error problems are from fmin or the null space/svd stuff.  
>> Once I find the input that drives my matrix to have a null space, I 
>> find the vector that corresponds to the null space (assuming the 
>> sub-matrix is only rank deficient by 1).  Then I combined the null 
>> space vector with zeros that correspond with the boundary condition 
>> on one end of the problem.  So, a 4x1 null space vector would give me 
>> an 8x1 full vector.  I then take the 8x8 full matrix which has a 4x4 
>> submatrix whose determinant is roughly 0 and multiply it by the 8x1 
>> vector of the boudnary conditions I just solved for.  This should 
>> then give me an 8x1 vector of the boundary conditions on the other 
>> end.  The problem is that there are some elements of this second 
>> vector that should be 0 because of the boundary conditions and they 
>> are actually of order 1e-10, if the vector is normalized so that its 
>> magnitude is 1.  This physically means that I have a cantilever beam 
>> with a free end that has just a little bit of force and moment at the 
>> free tip.
>>
>> Ryan
>
>
> Ryan,
>
> Please can you send me your matrix (matrices). I will try solve your 
> problem asap.
>
> A good starting point for your problem is an article by Ram
> Transcendental eigenvalue problem and its applications
> AIAA Journal Vol.40 (2002) pp. 1402-1407
>
> Nils
>
>>
>> Ryan Krauss wrote:
>>
>>> The matrix is currently 4x4 but will grow to probably 6x6.  It is 
>>> definitely nonlinear.  The matrix contains sinh, cosh, sin, and 
>>> cos.  I am using the transfer matrix method to analyze structures.  
>>> When you say two-parameter, do you mean the real and imaginary part 
>>> of the independent variable?  I guess you are right that I don't 
>>> necessarily need to use the determinant.  In order to satisfy the 
>>> boundary conditions of the problem this 4x4 or 6x6 matrix (which is 
>>> really a submatrix of an 8x8 or 12x12) must have a null space.  So, 
>>> what would be the better thing to look for?  An eignevalue that 
>>> approaches 0?
>>>
>>> Ryan
>>>
>>> Nils Wagner wrote:
>>>
>>>> Ryan Krauss wrote:
>>>>
>>>>> I have a matrix that is a function of a complex valued input.  I 
>>>>> am trying to find that value of that input that drives the 
>>>>> determinant of the matrix to zero.  I am searching for this value 
>>>>> using fmin. The error I am trying to minimize is the 
>>>>> abs(det(complex matrix)).
>>>>
>>>>
>>>>
>>>>
>>>>
>>>> It's not a good idea to use the determinant directly since det(A) 
>>>> is a rapidly varying function. As far as I understand your problem,
>>>> you are interested in the solution of a two-parameter nonlinear 
>>>> eigenvalue problem. Is that correct ? How about the size of your 
>>>> complex matrix A ?
>>>>
>>>> Nils
>>>>
>>>>> I don't seem to be able to drive this error lower that roughly 
>>>>> 9e-17, regardless of the values for ftol and xtol I use.
>>>>>
>>>>> Am I hitting some internal limitation?  Are complex values by 
>>>>> default single or double precision?
>>>>>
>>>>> Thanks,
>>>>>
>>>>> Ryan
>>>>>
>>>>> _______________________________________________
>>>>> SciPy-user mailing list
>>>>> SciPy-user at scipy.net
>>>>> http://www.scipy.net/mailman/listinfo/scipy-user
>>>>
>>>>
>>>>
>>>>
>>>>
>>>>
>>>>
>>>>
>>>> _______________________________________________
>>>> SciPy-user mailing list
>>>> SciPy-user at scipy.net
>>>> http://www.scipy.net/mailman/listinfo/scipy-user
>>>>
>>>
>>> _______________________________________________
>>> SciPy-user mailing list
>>> SciPy-user at scipy.net
>>> http://www.scipy.net/mailman/listinfo/scipy-user
>>>
>>
>> _______________________________________________
>> SciPy-user mailing list
>> SciPy-user at scipy.net
>> http://www.scipy.net/mailman/listinfo/scipy-user
>
>
>
>
>
>
> _______________________________________________
> SciPy-user mailing list
> SciPy-user at scipy.net
> http://www.scipy.net/mailman/listinfo/scipy-user
>

-------------- next part --------------
import TMM
from scipy import pi, sqrt, arange, vectorize
import pylab, scipy
import pdb

def main():
    mybeam=SAMIIBeam()
    print('mybeam.elemtype='+str(mybeam.elemtype))
    print('mybeam.params='+str(mybeam.params))
    mysys=TMM.TMMSystem([mybeam])
    svect=arange(1,500,1)
    svect=svect*1j
    vec_eigerror=vectorize(mysys.EigError)
    evect=vec_eigerror(svect)
#    pdb.set_trace()
#    evect=[]
#    for curs in svect:
#        e=mysys.EigError(curs)
#        evect.append(e)
    pylab.figure(1)
    pylab.cla()
    pylab.plot(scipy.imag(svect),evect)
    pylab.ylim((0,100))
    pylab.show()
    mytol=1e-30
#    if 0:
    eig1=mysys.FindEig(40j,mytol)
    e1=mysys.EigError(eig1)
    print('eig1='+str(eig1))
    print('error1='+str(e1))
#    if 0:
    eig2=mysys.FindEig(50j,mytol)
    print('eig2='+str(eig2))
    eig3=scipy.optimize.fmin(mysys.EigError,250j,xtol=mytol,ftol=mytol)
    print('eig3='+str(eig3))
    eig4=scipy.optimize.fmin(mysys.EigError,300j,xtol=mytol,ftol=mytol)
    print('eig4='+str(eig4))
#    pdb.set_trace()
    ms1=mysys.FindModeShape(eig1)
    print('ms1='+str(ms1))
    ms2=mysys.FindModeShape(eig2)
    print('ms2='+str(ms2))
    return mysys,eig1

def SAMIIBeam():
    #rho=2700;#kg/m^3
    ro=0.1412875/2.0;#meters ==> 5.5625in outer diameter
    ri=0.1314188/2.0;#meters (Lynnane's theis p. 175)
    t=ro-ri;
    L=4.6482;#meters
    density=2710.0;#kg/m^3
    E=6.8948e10;#N/m^2
    A=pi*(pow(ro,2)-pow(ri,2));
    mu=A*density;
    I=pi/4.0*(pow(ro,4)-pow(ri,4));
    J=2.0*I
    w1=pow(1.875,2)*sqrt((E*I)/(mu*pow(L,4)))
    f1=w1/(2.0*pi)
    print('I='+str(I))
    print('mu='+str(mu))
    print('EI1='+str(E*I))
    print('1.5*EI1='+str(1.2*E*I))
#    print('EA='+str(E*A))
    Ao=pi*pow(ro,2)
    Ai=pi*pow(ri,2)
    muo=density*Ao
    mui=density*Ai
    m22=(muo*pow(ro,2)-mui*pow(ri,2))/4.0
    m33=m22
    m11=2.0*m22
    G=2.6e10;#Pa
    Ks=G*A
    TK=J*G
    myparams={'EI1':E*I,'EI2':1.5*E*I,'L':L,'mu':mu,'EA':E*A,'GJ':G*J,'m11':m11}
    return TMM.TMMElement('beam',myparams)

if __name__=="__main__":
    mysys,eig1=main()
-------------- 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))
        basevect=basevect*1j
#        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
        endvect=scipy.matrixmultiply(self.FindU(eigval),basevect)
        endvect=endvect/(scipy.linalg.norm(endvect))
        return endvect

    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)


More information about the SciPy-user mailing list