On Fri, 2015-06-19 at 08:49 +0000, Carrico, Paul wrote: (top posting fixed) > -----Message d'origine----- > De : users [mailto:[email protected]] De la part de Serge Steer > Envoyé : jeudi 18 juin 2015 10:08 À : [email protected] Objet : Re: > [Scilab-users] eigs calculation > > Le 17/06/2015 22:18, Carrico, Paul a écrit : > > Dear All > > > > Thanks for the answers. > > > > To give more information's on what I'm doing (That's quite new I confess), > > I'm performing a (basic) finite element calculation with beams using > > sparse matrix (stiffness matrix K and mass matrix M). > > [u,v] = > > eigs(K((ddl_elem+1):$,(ddl_elem+1):$),M((ddl_elem+1):$,(ddl_elem+1):$) > > ,n,"SM"); > > > > Nota: ddl means dof > > > > I'm calculated first the natural frequencies using (K - omega^2.M).x=0 ... > > the pulse (or circular frequencies) are no more and no less than the > > eigenvalues of the above system (u = omega^2). > > > > Just a "mechanical" remark: since the beam is clamped in one side (and free > > on the tip), it is absolutely normal that you find twice the same natural > > frequency (1rst mode in one direction, the second one in a new direction at > > 90°) .... I've been really surprised to find it, but happy at the same time > > ... > > > > The origin of my question was: since it obvious that the results are > > strongly sensitive to the "units" (i.e. the numbers), I'm wondering if > > there is a way to control the accuracy of the eigenvalues calculation using > > eigs keywords ... > There is no way to improve accurary with an option. In general the numerical > algorithms try to return the best solution. > But it should be possible to improve accuracy by a well suited normalisation > (unit change for example!) > > The condition number of u gives a measure of the numerical difficulty to > solve the problem > > > Note that as stated the eigenvalue computation may be a ill conditionned > problem in particular when there are clusters of eigenvalues > > Please find below a little script which illustrate a typicall case > A=zeros(10,10)+diag(ones(1,9),1);A(10,1)=%eps; > s=spec(A); > clf;plot(real(s),imag(s),'+'); > > Serge Steer > > In any way, thanks for the debate > > > > Paul Dear > > I've been reading all the advices that (all of them) speak about the need to normalize the data in order to avoid the ill-conditioned issues (and to be the closest as possible to the value 1.0) ; the example shared by Serge Steer is quite interesting and highlights the ill-conditioning issue. > > Nevertheless I'm a bit confuse on how to do it in order to rebuild the results afterward: > [u,v] = eigs(K((ddl_elem+1):$,(ddl_elem+1):$),M((ddl_elem+1): $,(ddl_elem+1):$),n,"SM"); > > pulse_ =sqrt(diag(u)) //pulse - circular frequency > freq_ = pulse_/(2*%pi); //natural frequency > > Naively I thought that the following would have fix the problem: > - using max(K) and max(M) to normalize the sparse matrix, > - removing the smallest values (i.e. bellow a threshold value) > > ... but the matrix becomes singular (why ?). > > So is there any additional advice ? > > Paul > "Normalization" refers to jiggering the numbers around in a way that does not change the problem in a strictly mathematical sense, but which makes it more tractable.
d = eigs(A, B) computes the solutions to A * v = lambda * B * v So any nonsingular square matrix N won't change the problem if it's multiplied in: N * A * v = N * lambda * B * v Because lambda is a scalar (well, I hope I'm getting that right) you can change the problem to A_ = N * A, B_ = N * B, and d = eigs(A_, B_) will, theoretically, get the same answers as with the original matrices, but possibly with better numerical conditioning. This is what you're doing when you change units. If you don't mind the meaning of your eigenvectors changing, you can do a similarity transform. Start with N * A * N^(-1) * N * v = N * lambda * B * N^(-1) * N * v Now set A_ = N * A * N^(-1) B_ = N * B * N^(-1) v_ = N * v [d, v_] = eigs(A_, B_); v = N^(-1) * v_; will, again, theoretically give you the same numbers as before, but it may be better conditioned numerically. Actually _finding_ N, or giving you advise on how to do so, is beyond my powers -- but maybe this will set you on a better road. -- Tim Wescott www.wescottdesign.com Control & Communications systems, circuit & software design. Phone: 503.631.7815 Cell: 503.349.8432 _______________________________________________ users mailing list [email protected] http://lists.scilab.org/mailman/listinfo/users
