[Scilab-users] eigs calculation

Thu Jun 18 10:38:02 CEST 2015

Le 06/17/2015 10:18 PM, 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 ...
>
> In any way, thanks for the debate

Hi Paul,

I won't go into the ill-conditionned matrix inversion, I don't know much 
about it.
What I say below might be out of topic but I might also be useful.
Here is a rule: whenever you use floating point arithmetic: stay close 
to one!
That is, apply whatever normalization makes sense for your problem so 
that the numbers your are crunching are around 1 (your matrix elements 
in your case).
Never, never use "units" in your actual calculation, just apply the 
proper renormalization to the result of your number crunching.
The rational behind this rule is that the distance to the next float (or 
double) you can represent scales with the float value (see attached plot).
Matlab has a function eps to measure this distance: 
http://fr.mathworks.com/help/matlab/ref/eps.html (no equivalent in scilab)
You should have a look here for more fun with floats: 
http://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html

Cheers,

Antoine

>
> Paul
>
> -----Message d'origine-----
> De : users [mailto:users-bounces at lists.scilab.org] De la part de roger.cormier at ncf.ca
> Envoyé : 17 June 2015 19:50
> À : International users mailing list for Scilab.
> Objet : Re: [Scilab-users] eigs calculation
>
> Paul:
>
> Did you consider checking the condition number of the matrix?
>
> What I sometimes do, for problems that are ill conditioned by their very nature, is normalise the entries in an attempt to be as far away from singularity (i.e. to keep the rows and columns as orthogonal as possible) and then do a check on the condition number to see if there is an improvement. The lower the condition number the better. Once your calculations are done you can convert back to the un-normalised entries, does improve a bit when inverting matrices of ill-conditioned problems such as target-motion analysis.
>
> I'm willing to bet that the condition number of the second matrix is a bit lower than of the first one. Also, I'm noticing that there are  four double poles, possibly six in total (resonnances 5 and 6, and 11 and 12), in the second calculation. You knowing your system, that could provide you with some hints as to which calculation to trust.
>
> Let me know how you make out.
>
> Regards,
>
>
> Roger.
>   
>
> ___________________________
> Dr. Roger Cormier, P.Eng.
>
>
> Le mer. 17 juin 2015 à 11:55, tim at wescottdesign.com a écrit :
>
>> On 2015-06-17 06:50, Carrico, Paul wrote:
>>> Dear All,
>>> I'm performing a (mechanical) calculation using the eigs and I've been
>>> noticing that the results are strongly sensitive on the unit system
>>> I'm using; I can understand that high numbers can lead to some
>>> numerical "issues" .
>>> Is there a way to increase the accuracy ?
>>> Paul
>>> PS: the 2 types of results
>>> _NB_:
>>> 1 (MPa) = 1E6 (Pa)
>>> 1 (mm) = 1E-3 (m)
>>> 1 (Kg/m^3) = 1E12 (T/mm^3)
>>> [u,v] =
>>> eigs(K((ddl_elem+1):$,(ddl_elem+1):$),M((ddl_elem+1):$,(ddl_elem+1):$),n,"SM");
>>> a) calculation 1 in Pa, m, Kg/m^3
>>> Natural frequency calculation:
>>> - Resonance 1 : 497.956 Hz
>>> - Resonance 2 : 3120.64 Hz
>>> - Resonance 3 : 5277.8 Hz
>>> - Resonance 4 : 6948.69 Hz
>>> - Resonance 5 : 8737.88 Hz
>>> - Resonance 6 : 15832.1 Hz
>>> - Resonance 7 : 17122.8 Hz
>>> - Resonance 8 : 20847.8 Hz
>>> - Resonance 9 : 26382.5 Hz
>>> - Resonance 10 : 28305.1 Hz
>>> - Resonance 11 : 34752 Hz
>>> - Resonance 12 : 36926.4 Hz
>>> b) Calculation in MPa, mm, T/mm^3 ..
>>> Natural frequency calculation:
>>> - Resonance 1 : 497.955 Hz
>>> - Resonance 2 : 497.956 Hz
>>> - Resonance 3 : 3120.59 Hz
>>> - Resonance 4 : 3120.64 Hz
>>> - Resonance 5 : 6948.69 Hz
>>> - Resonance 6 : 7463.93 Hz
>>> - Resonance 7 : 8737.56 Hz
>>> - Resonance 8 : 8737.88 Hz
>>> - Resonance 9 : 17121.6 Hz
>>> - Resonance 10 : 17122.8 Hz
>>> - Resonance 11 : 20847.8 Hz
>>> - Resonance 12 : 22390 Hz
>> Hi Paul:
>>
>> I can't tell you about the innards of Scilab specifically, but eigenvalue calculation in general can be very sensitive to numerical issues.  If you're entering the data by hand or otherwise truncating the source data your entire difference in results may just be from rounding error in your source data.
>>
>> If you're starting from one set of source data and multiplying by conversion constants, then you can try changing the tolerance (if it's not in the function then there's a global one, called, I think, %TOL).
>>
>> There are ways to make the matrices more numerically stable.  I am absolutely positively not an expert on this, but I think that the more you can make your matrix into something with a band of non-zero numbers around the main diagonal and zeros elsewhere, the more stable the problem will be.
>>
>> If mechanical systems are like control systems, then the numerical stability of the matrix that describes the system dynamics is just a reflection of the real sensitivity of the real system to manufacturing variations -- it may be that, in a group of several physical units all assembled to the same specification, you'll find that much variation in the real world!
>>
>> _______________________________________________
>> users mailing list
>> users at lists.scilab.org
>> https://urldefense.proofpoint.com/v2/url?u=http-3A__lists.scilab.org_mailman_listinfo_users&d=AwIF-g&c=0hKVUfnuoBozYN8UvxPA-w&r=4TCz--8bXfJhZZvIxJAemAJyz7Vfx78XvgYu3LN7eLo&m=X5rGSmMzj_mrDPhlxODiYEfNco7jW0285o_vP15-sDI&s=AUkkoXcGhr1-zgfahFwMmtOJuuCDaaCeKtS2-qHM8NQ&e=
> _______________________________________________
> users mailing list
> users at lists.scilab.org
> https://urldefense.proofpoint.com/v2/url?u=http-3A__lists.scilab.org_mailman_listinfo_users&d=AwIF-g&c=0hKVUfnuoBozYN8UvxPA-w&r=4TCz--8bXfJhZZvIxJAemAJyz7Vfx78XvgYu3LN7eLo&m=X5rGSmMzj_mrDPhlxODiYEfNco7jW0285o_vP15-sDI&s=AUkkoXcGhr1-zgfahFwMmtOJuuCDaaCeKtS2-qHM8NQ&e=
>
> EXPORT CONTROL :
> Cet email ne contient pas de données techniques
> This email does not contain technical data
> _______________________________________________
> users mailing list
> users at lists.scilab.org
> http://lists.scilab.org/mailman/listinfo/users
>

-- 
+++++++++++++++++++++++++++++++++++++++++++++++++++++++

  Antoine Monmayrant LAAS - CNRS
  7 avenue du Colonel Roche
  BP 54200
  31031 TOULOUSE Cedex 4
  FRANCE

  Tel:+33 5 61 33 64 59

  email : antoine.monmayrant at laas.fr
  permanent email : antoine.monmayrant at polytechnique.org

+++++++++++++++++++++++++++++++++++++++++++++++++++++++

-------------- next part --------------
A non-text attachment was scrubbed...
Name: untitled.png
Type: image/png
Size: 16440 bytes
Desc: not available
URL: <https://lists.scilab.org/pipermail/users/attachments/20150618/52ae7fe3/attachment.png>