Prev: -- Packing unit circles in squares: new results
Next: On the defense of cranks, Re: for x,y > 7 twins(x+y) <= twins(x)+ twins(y)
From: Peter Spellucci on 25 Apr 2010 11:59 In article <4bd2c77b$0$30906$5fc30a8(a)news.tiscali.it>, matt <none(a)none.com> writes: >Dear all, > > I would like to be able to estimate _numerically_ the eigenvalues (ev) >of a matrix M(x) as x changes, and display them in a plot. Let me use >Matlab notation to further clarify what I'm trying to achieve. > >Let us suppose that I have M = [x-1, 3; x^2, x+4]. I would give some >values to x, e.g x = -3:.01:3, and for each of these values I would >update Mx_ev(:,i) = eig([x(i)-1, 3; x(i)^2, x(i)+4]), being Mx_ev a >matrix of 2x601 storing the evolution of the ev of M(x) as x changes. >Unfortunately, nothing guarantees that the "corresponding" ev always >appear in the same position at different values of x: they may well >exchange their positions for some x, showing crossing on the plot at >those x. > >An obvious workaround to this is explicit sorting. But I would like to >avoid it, because requires a significant amount of computation. And by >the way, how am I supposed to check the occurrency of discontinuties in >the ev evolutions, without knowing their rate of rise (or fall) in >advance? Although checking for a variation of 2-3% from the previous >value worked correctly for all the cases I tried, I suppose there must >be a better (more accurated) way to checking for discontinuities. > >So the question is: do you know any algorithm that, starting from very >good estimations, Mx_ev(:,i), computes the ev of M(x(i+1)), namely >Mx_ev(:,i+1), with the same order in which the estimates are given? > >I tried the Shifted-QR algorithm with deflation, but without success. As >an example, consider the following matrix (A) standing for M(x(i+1)). > >A = > > 0.0015 - 0.1645i -0.0054 + 0.0012i > 0.0032 0.0626 + 0.0122i > >We have the following estimate (s) of its second ev, i.e. we know that >the second ev of M(x(i)), namely Mx_ev(2,i), is > >s = > > 0.0016 - 0.1684i > >What I was expecting is that, after some iteration, the Shifted-QR >algorithm > > >> n = size(A,1); I = eye(n,n); > >> [Q,R] = qr(A - s*I); A = R*Q + s*I, s = A(n,n), > >returned 0.0015 - 0.1646i, i.e. the continue evolution at x(i+1) of the >second ev of M(x(i)), namely Mx_ev(2,i+1), in A(2,2). Instead, I got > >A = > > 0.0015 - 0.1646i -0.0019 + 0.0024i > 0 0.0626 + 0.0123i > >which shows clearly that the initial estimate 0.0016 - 0.1684i converged >to 0.0626 + 0.0123i (the other ev) instead of 0.0015 - 0.1646i first of all the good message: eigenvalues are Hoelder continuous functions of the matrix elements (hence your x) of index 1/n at least. Normally as long as they are all pairwise different, even smooth. I wonder what you mean by "sorting" for complex eigenvalues. this makes no sense, the more as your curves you obviously think of (you want a n curves in the complex plane, each one representing the development of one eigenvalue as function of x, right?) may cross themselves. Hence the advice is to use small steps and extrapolate the curves using an appropriate model for it, e.g. polynomial interpolation of low degree (quadratic, cubic) of back values for "ordering" the new set. (also a double eigenvalue can develop Lipschitz continuous) In the case of a switch from real ->real double -> complex conjugate this does no longer work, which branches to fit together is a question of convention. hth peter |