[ncl-talk] Errors with inverting a matrix

[Dennis Shea]

Most numerical analysis experts would agree with what was written at:
*dont-invert-that-matrix*
<https://www.johndcook.com/blog/2010/01/19/dont-invert-that-matrix/>

In fact, if the *determinant*
<http://www.ncl.ucar.edu/Document/Functions/Built-in/determinant.shtml> of
a matrix is 0.0, it is not possible to invert the matrix.
The matrix you are using has a determinant value of -1.02827e-41. Possibly,
only numerical rounding issues prevented the determinant value from being
exactly 0.0. I guess the matrix would be characterized as ill-conditioned.

Matlab likely uses *much more sophisticated *approaches than using the
LAPACK's LU factorization method used by *inverse_matrix *
<http://www.ncl.ucar.edu/Document/Functions/Built-in/inverse_matrix.shtml>.That
is what Matlab does ... and VERY well!

If you need the matrix inverse, perhaps use Matlab and write out the
correct matrix inversion to a file and read that into NCL?

===

%> ncl invert.ncl_talk.ncl


On Fri, Oct 25, 2019 at 7:31 PM Prashanth Bhalachandran via ncl-talk <
ncl-talk at ucar.edu> wrote:

> Dear NCL'ers,
> I am having a problem with inverting the matrix using NCL that does not
> make any sense.
>
> I double-checked both the routines: inverse_matrix() as well as
> solve_linsys() and both are functioning well for sample matrices that I am
> trying out. I even compared their answers to Matlab and they were perfectly
> right. However, when I use any of the two routines for the below matrix A:
>
> write_matrix(A,"9f11.7", False)
>
>
> A =
>
>   0.4000000 -0.0000000 -0.1000000 -0.0000000 -0.1000000  0.3077683
> -0.0000000  0.0726542 -0.0000000
>
>  -0.0000000  0.1500000 -0.0000000 -0.1000000 -0.0000000 -0.0000000  0.1902113
> -0.0000000  0.0363271
>
>  -0.1000000 -0.0000000  0.1500000  0.0000000 -0.1000000 -0.1175570
> -0.0000000  0.1538842 -0.0000000
>
>  -0.0000000 -0.1000000  0.0000000  0.1500000 -0.0000000 -0.0000000
> -0.1538842 -0.0000000  0.1175570
>
>  -0.1000000 -0.0000000 -0.1000000 -0.0000000  0.1500000 -0.0363271
> -0.0000000 -0.1902113 -0.0000000
>
>   0.3077683 -0.0000000 -0.1175570 -0.0000000 -0.0363271  0.2500000  0.0000000
> -0.0000000 -0.0000000
>
>  -0.0000000  0.1902113 -0.0000000 -0.1538842 -0.0000000  0.0000000
> 0.2500000  0.0000000  0.0000000
>
>   0.0726542 -0.0000000  0.1538842 -0.0000000 -0.1902113 -0.0000000
> 0.0000000  0.2500000 -0.0000000
>
>  -0.0000000  0.0363271 -0.0000000  0.1175570 -0.0000000 -0.0000000  0.0000000
> -0.0000000  0.2500000
>
>
> The inverse that NCL computes is totally wrong. That is: One, it does not
> satisfy the identity, A Ainv = Identity Matrix. Two, for the same values,
> it is very different from what Matlab produces and the one that Matlab
> produces satisfies the identity matrix test. Can someone please help me
> understand as to what is wrong with the above matrix and how to circumvent
> this problem?
>
> I read in previous NCL help that I could use solve_linsys(A, Identity
> matrix) and that should give me the inverse. However, that does not work
> for this particular A as well. It does work for other sample As that I
> tried.
>
> Please offer any advice if you have run into something like this before.
>
> Warm Regards,
> Prashanth
> _______________________________________________
> ncl-talk mailing list
> ncl-talk at ucar.edu
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20191027/a9473f51/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: invert.ncl_talk.ncl
Type: application/octet-stream
Size: 3319 bytes
Desc: not available
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20191027/a9473f51/attachment.obj>