Hi Markus, Thanks for your patch. Would you mind sending us a pull request against the official repository?
https://github.com/sympy/sympy Then we can easily review it and push it in. Thanks, Ondrej On Sat, Sep 17, 2011 at 7:41 AM, Markus Müller <markus.mueller....@googlemail.com> wrote: > Hallo > Since I needed not only the jordan normal form but also the > transformation I had a look at the code. > (How nice to have open source CAS, Thanks for that) > I tried to combine the computation of P and J and therefore had to > look very closely. > So I realized the bug > It can be expressed very easy: > In sympy 0.71 only the geometric and the algebraic multiplicity > are used to determine the blocksizes. > > This is expressed by the interface of function: > self._jordan_split(multiplicity, geometrical) > in the matrix class. > However this is not possible in general since there exist different > Jordan matrices > (with different blocksizes ) that have the same eigenvectors with the > same geometrical multiplicity as an counter example shows which I > added in the comment for the patch. > > I hope that helps (took me a week nearly inspite of my plan to > implement it in 4h...) > > Best regards > Markus > > So here comes the patch with a lot of comments how it works and some > testing code around it > I inherited from matrix and overloaded the jordan_cells method: > > > #!/usr/bin/python > from sympy.interactive.printing import init_printing > init_printing(use_unicode=True, wrap_line=False, no_global=True) > from sympy.matrices import * > def jordan_cell(eigenval, n): > n = int(n) > out = zeros(n) > for i in range(n-1): > out[i, i] = eigenval > out[i, i+1] = 1 > out[n-1, n-1] = eigenval > return out > > class Matrix2(Matrix): > def jordan_cells(self, calc_transformation = True): > n=self.rows > Jcells = [] > Pcols_new=[] > # to every eingenvalue may belong i blocks with size s(i) > # and a chain of generalized eigenvectors > # which will be determined by the following computations > # for every eigenvalue we will add a dictionary > # containing for all blocks the blockssizes and the attached > chain vectors > # that will eventually be used to form the transformation P > JordanBlockStructures={} > _eigenvects=self.eigenvects() > for eigenval, multiplicity, vects in _eigenvects: > geometrical = len(vects) > if geometrical == multiplicity: > Jcell = diag( *([eigenval] * multiplicity)) > Jcells.append(Jcell) > Pcols+=vects > > elif geometrical==0: > raise MatrixError("Matrix has the eigen vector with > geometrical multiplicity equal zero.") > else: > > # Up to now we know nothing about the sizes of the blocks of our > Jordan matrix > # Note that knowledge of algebraic and geometrical multiplicity > # will >> n o t << be sufficient to determin this structure > # The blocksize s could be defined as the minimal k where > # (self-lI)^k=0 (Since (self-lI) is a nilpotent matrix such k is > bound to exist. > # The worst case would be that k= (multiplicity-geometrical+1) > but the blocks could be smaller. > # Consider for instance the following matrix > # [2 1 0 0] > # [0 2 1 0] > # [0 0 2 0] > # [0 0 0 2] > # which coincides with it own Jordan canonical form. > # It has only one eigenvalue l=2 of (algebraic) multiplicity=4. > # It has two eigenvectors. one belonging to the last row > (blocksize 1) > # and one being the last part of a jordan chain of length 3 > (blocksize of the first block) > # Note again that it is not not possible to obtain this from the > algebraic and geometrical > # multiplicity alone. This only gives us an upper limit for the > dimension of one of > # the subspaces (blocksize of according jordan block) given by > # max=(multiplicity-geometrical+1) which is reached for our > matrix > # but not for > # [2 1 0 0] > # [0 2 0 0] > # [0 0 2 1] > # [0 0 0 2] > # although multiplicity =4 and geometrical=2 are the same for > this matrix. > # Actually we have to increase k until (A-lI)^k=0 > Z=zeros(self.rows) > I=eye(self.rows) > l=eigenval > M=(self-l*I) > # We will store the matrices (self-l*I)^k for further > computations > # for convinience only we store Ms[0]=(sefl-lI)^0=I > # so the index is the same as the power for all further Ms > entries > # We also store the vectors that span these kernels (Ns[0] =[] ) > # and also their dimensions a_s > # this is mainly done for debugging since the number of blocks of > a given size > # can be computed from the a_s, in order to check our result > which is obtained simpler > # by counting the number of jordanchains for a given s > # a_0 is dim(Kernel(Ms[0]) =dim (Kernel(I))=0 since I is regular > smax=0 > Ms=[I] > Ns=[[]] > a=[0] > l_JordanChains={} > l_JordanBlockNumbers={} > ChainVectors=[] > while Ms[-1] != Z: > M_new=Ms[-1]*M > Ns_new=M_new.nullspace() > a_new=len(Ns_new) > Ms.append(M_new) > Ns.append(Ns_new) > a.append(a_new) > smax+=1 > # We now have Ms[-1]=((self-l*I)**s)=Z=0 > # We now know the size of the biggest jordan block > # associatet with l to be s > # now let us proceed with the computation of the associate part > of the transformation matrix P > # We already know the kernel (=nullspace) K_l of (self-lI) which > consists of the > # eigenvectors belonging to eigenvalue l > # The dimension of this space is the geometric multiplicity of > eigenvalue l. > # For every eigenvector ev out of K_l, there exists a subspace > that is > # spanned by the jordan chain of ev. The dimension of this > subspace is > # represented by the length s of the jordan block. > # The chain itself is given by {e_0,..,e_s-1} where: > # e_k+1 =(A-lI)e_k (*) > # and e_s-1=ev > # So it would be possible to start with the already known ev and > work backwards until one > # reaches e_0. Unfortunately this can not be done by simply > solving system (*) since its matrix > # is singular (by definition of the eigenspaces) > # This approach would force us a choose in every step the degree > of freedom undetermined > # by (*). This is difficult to implement with computer algebra > systems and also quite unefficient > # We therefore reformulate the problem in terms of nullspaces. > # To do so we start from the other end and choose e0's out of > # E=Kernel(self-lI)^s / Kernel(self-lI)^(s-1) > # Note that Kernel(self-lI)^s =Kernel(Z)=V (the whole vector > space) > # So in the firs step s=smax this restriction turns out to > actually restrict nothing at all > # and the only remaining condition is to chose vectors in > Kernel(self-lI)^(s-1) > # Subsequently we compute e_1=(self-lI)e_0, e_2=(self-lI)*e_1 and > so on. > # The subspace E can have a dimension larger than one > # That means that we have more than one Jordanblocks of size s > for the eigenvalue l > # and as many jordanchains (This is the case in the second > example) > # In this case we start as many jordan chains and have as many > blocks of size s in the jcf > # We now have all the jordanblocks of size s but there might be > others attached to the same > # eigenvalue that are smaller. > # So we will do the same procedure also for s-1 and so on until 1 > the lowest possible order > # where the jordanchain is of lenght 1 and just represented by > the eigenvector > for s in reversed(xrange(1,smax+1)): > #we want the vectors in Kernel((self-lI)^s) (**) > S=Ms[s] > # but without those in Kernel(self-lI)^s-1 so we will add > these as additional equations > # to the sytem formed by S (S will no longer be quadratic but > this does not harm > # since S is rank deficiant) > excludeVectors=Ns[s-1] > for k in range(0,a[s-1]): > S=S.col_join((excludeVectors[k]).transpose()) > # we also want to exclude the vectors in the chains for the > bigger blogs > # that we have already computed (if there are any) > # (That is why we start wiht the biggest s) > ######## implementation remark:############# > # Doing so for >> a l l << already computed chain vectors > # we actually exclude some vectors twice because they are > already excluded > # by the condition (**) > # This happens if there are more than one blocks attached to > the same eigenvalue >>and<< > # the current blocksize is smaller than the block whose chain > vectors we exclude > # If the current block has size s_i and the next bigger block > has size s_i-1 then > # the first s_i-s_i-1 chainvectors of the bigger block are > allready excluded by (**) > # The unnecassary adding of these equations could be avoided > if the algorithm would > # take into account the lengths of the already computed chains > which are already stored > # and add only the last s items. > # However the following loop would be a good deal more nested > to do so. > # > # A much more elegant alternative approach might be to drop > condition (**) altogether > # because it is added implicitly by excluding the chainvectors > # > # Since adding a linear dependent equation does not change the > result > # it can harm only in terms of efficiency. > # So to be sure i let it there for the moment > # > l=len(ChainVectors) > if l>0: > for k in range(0,l): > old=ChainVectors[k].transpose() > S=S.col_join(old) > > e0s=S.nullspace() > #determine the number of chain leaders which equals the number > of blocks with that size > n_e0=len(e0s) > s_chains=[] > s_cells=[] > for i in range(0,n_e0): > chain=[e0s[i]] > for k in range(1,s): > v=M*chain[k-1] > chain.append(v) > > chain.reverse() #we want the chain leader appear as the last > of > the block > ChainVectors+=chain > s_chains.append(chain) > l_JordanChains[s]=s_chains > JordanBlockStructures[eigenval]=l_JordanChains > P=zeros(n) > for l in reversed(sorted(JordanBlockStructures.keys())): #start > with the biggest eigenvalue > l_JordanChains=JordanBlockStructures[l] > for s in reversed(sorted((l_JordanChains).keys())): #start with the > biggest block > s_chains=l_JordanChains[s] > block= jordan_cell(eigenval,s) > number_of_s_chains=len(s_chains) > for i in range(0,number_of_s_chains): > Jcells.append(block) > chainVectors=s_chains[i] > lc=len(chainVectors) > assert(lc==s) > for j in range(0,lc): > generalizedEigenVector=chainVectors[j] > Pcols_new.append(generalizedEigenVector) > > J = diag(*Jcells) > P=zeros(n) > for j in range(0,n): > P[:,j]=Pcols_new[j] > assert(J==P.inv()*self*P) > return (P,Jcells) > > #A=Matrix2(5,5, > # [5,1,0,0,0, > # 0,5,1,0,0, > # 0,0,5,0,0, > # 0,0,0,5,1, > # 0,0,0,0,5] > # ) > A=Matrix2(7,7, > [5,1,0,0,0,0,0, > 0,5,1,0,0,0,0, > 0,0,5,1,0,0,0, > 0,0,0,5,0,0,0, > 0,0,0,0,5,0,0, > 0,0,0,0,0,5,1, > 0,0,0,0,0,0,5] > ) > #A=Matrix2(5,5, > # [25, -16, 30, -44, -12, > # 13, -7, 18, -26, -6, > # -18, 12, -21, 36, 12, > # -9, 6, -12, 21, 6, > # 11, -8, 15, -22, -3] > #) > > > print(A) > (P,cells)=A.jordan_cells() > print "P=",P > print "P^-1 A P=" > print P.inv()*A*P > > > -- > You received this message because you are subscribed to the Google Groups > "sympy-patches" group. > To post to this group, send email to sympy-patches@googlegroups.com. > To unsubscribe from this group, send email to > sympy-patches+unsubscr...@googlegroups.com. > For more options, visit this group at > http://groups.google.com/group/sympy-patches?hl=en. > > -- You received this message because you are subscribed to the Google Groups "sympy-patches" group. To post to this group, send email to sympy-patches@googlegroups.com. To unsubscribe from this group, send email to sympy-patches+unsubscr...@googlegroups.com. For more options, visit this group at http://groups.google.com/group/sympy-patches?hl=en.
