[Scilab-users] eigs calculation

[Tim Wescott]

On Fri, 2015-06-19 at 08:49 +0000, Carrico, Paul wrote:
(top posting fixed)
> -----Message d'origine-----
> De : users [mailto:users-bounces at lists.scilab.org] De la part de Serge Steer Envoyé : jeudi 18 juin 2015 10:08 À : users at lists.scilab.org 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