| IDL inverse matrix problem?? [message #90065] |
Thu, 22 January 2015 11:24  |
amin farhang
Messages: 39
Registered: November 2010
|
Member |
|
|
Dear All,
I'm trying to inverse a big matrix with heavy elements. According to http://www.idlcoyote.com/math_tips/sky_is_falling.html discussion by David, I try to produce all elements in double precision. But the INVERT and LA_INVERT return wrong matrices. For example:
First define a double precision matrix:
IDL> m1 = 9.5d15 + 0.5d18*randomu(1,1,5,/double)
IDL> m2 = -1d0 + 2d0*randomu(2,5,1,/double)
IDL> Am = m1 ## m2
IDL> print,Am
-6.8599934e+16 -8.3590379e+16 -1.8742199e+17 1.5697952e+17 1.1607144e+17
-1.2228319e+17 -1.4900449e+17 -3.3409009e+17 2.7982470e+17 2.0690379e+17
-1.4972730e+17 -1.8244569e+17 -4.0907020e+17 3.4262599e+17 2.5333938e+17
-1.0612015e+17 -1.2930950e+17 -2.8993102e+17 2.4283828e+17 1.7955584e+17
-1.1696076e+17 -1.4251900e+17 -3.1954867e+17 2.6764521e+17 1.9789821e+17
IDL> max_value = double(max(abs(Am)))
IDL> Am2 = la_invert(Am/max_value)/(max_value)
% LA_INVERT: Singular matrix encountered, STATUS=4.
% Execution halted at: $MAIN$
As you see la_invert encounter error. Since the matrix determinant is big.
Now again examine the following definition:
IDL> Am = [[-6.8599934d+16,-8.3590379d+16,-1.8742199d+17,1.5697952d+17 , 1.1607144d+17],$
[-1.2228319d+17,-1.4900449d+17,-3.3409009d+17,2.7982470d+17, 2.0690379d+17],$
[-1.4972730d+17,-1.8244569d+17,-4.0907020d+17,3.4262599d+17, 2.5333938d+17],$
[-1.0612015d+17,-1.2930950d+17,-2.8993102d+17,2.4283828d+17, 1.7955584d+17],$
[-1.1696076d+17,-1.4251900d+17,-3.1954867d+17,2.6764521d+17, 1.9789821d+17]]
IDL> max_value = double(max(abs(Am)))
IDL> Am2 = la_invert(Am/max_value)/(max_value)
As you see this time IDL return the inverse BUT this inverse is WRONG!!!
Let me test the IDL inversion:
IDL> Bm = Am ## Am2
IDL> print,Bm
0.99999999 -1.3038516e-08 -5.5879354e-09 -4.4703484e-08 4.4703484e-08
-7.4505806e-09 1.0000000 1.8626451e-09 -4.4703484e-08 5.9604645e-08 -7.4505806e-09 -7.4505806e-09 0.99999999 -8.9406967e-08 1.1920929e-07
-7.4505806e-09 -1.1175871e-08 -7.4505806e-09 0.99999994 1.7881393e-07
-7.4505806e-09 -1.1175871e-08 -4.6566129e-09 -1.1920929e-07 1.0000001
As you see the returned matrix is not a identity matrix, therefore the inversion is wrong.
The correct inversion is:
1.0e-08 *
0.0166 -0.0010 -0.0008 0.0277 -0.0228
0.0292 -0.3699 0.0482 0.0017 -0.0105
-0.0221 0.0376 0.0224 -0.0485 0.0154
-0.0100 -0.0207 0.0110 -0.0025 -0.0005
-0.0103 0.0520 -0.0125 0.0092 0.0035
Am I doing something wrong?
Best regards,
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
