How can I deal with a badly conditionned matrix when computing strain energy?

I am trying to compute strain energy with a program based on a Rayleigh-Ritz method. The base I used is not orthogonal and therefore, my stiffness matrix, CK, is badly conditionned (rcond = 8.02e-21).

It's not that much of a problem to solve the system and obtain the complex displacement vector Xmn, but it is a bigger problem to computer the strain energy:

 SE = 1/2*Xmn'*CK*Xmn;

If CK is real and symetric, SE is complex when it should be real.

If I replace CK by a normaly conditionned real symetric matrix, then SE is real.

Does anybody have any trick help me to compute SE such that it is real instead of complex?

Answer 1

[Andrew Newell]

The specifics about strain are beside the point - it is a rounding problem. Take any badly-conditioned matrix and some complex vector and you will get a similar result:

>> A = diag([1e-10 1 1e10]);

>>  rcond(A)

ans =

   1.0000e-20
>>  chi = rand(3,1) + sqrt(-1)*rand(3,1);

>>  chi'*A*chi

ans =

   1.0985e+10 + 4.7684e-07i

Even though the imaginary part looks small enough to ignore, it is a symptom of a more serious problem: the real part is inaccurate too!

>> B = sym(A); v = sym(chi); double(v'*B*v)

ans =

   1.2952e+10

So if you need accuracy, do the calculation using the Symbolic Toolbox. However, the badly conditioned matrix may be a sign that your problem itself is ill-defined and needs to be revised.