inverse Matrix problem, Newton's method

[mecht]

I'm writing some Octave code to solve a finite element-like problem using the
Newton's method.
My Jacobian matrix looks like this (spy(J)):
<img src="https://preview.ibb.co/hknY5A/spy.jpg"; alt="spy" border="0" />
<https://ibb.co/fCC0kA>  
Diagonal elements are negative sums of other elements in a row.
My problem is that i can not produce any result due to false(?) matrix
singularity. I tried left-hand matrix division, LU method, inv(A) but Octave
warns about matrix singularity due to machine precision.
Octave returns 0 or extremely low value (like 1.0E-100) when asked for
det(A)
At first glance, i thought the reason is that if I multiply large amount
diagonal values of 1.0E-1 order the product falls under machine precision.
BUT
I tried my code with minimal number of finite elements 2x2x2, which gives
Jacobian on 8x8 size, and the result was not that low.

Result example:

  11.20161  -12.53016    0.00000    0.00000    1.32856    0.00000    0.00000   
0.00000
  -12.53016    8.54449    0.00000    0.00000    0.00000    3.98567   
0.00000    0.00000
    0.00000    0.00000  -13.85872   12.53016    0.00000    0.00000   
1.32856    0.00000
    0.00000    0.00000   12.53016  -16.51584    0.00000    0.00000   
0.00000    3.98567
    1.32856    0.00000    0.00000    0.00000  -13.85872   12.53016   
0.00000    0.00000
    0.00000    3.98567    0.00000    0.00000   12.53016  -16.51584   
0.00000    0.00000
    0.00000    0.00000    1.32856    0.00000    0.00000    0.00000 
-13.85872   12.53016
    0.00000    0.00000    0.00000    3.98567    0.00000    0.00000  
12.53016  -16.51584

Octave says, that det(J):

ans = 0


circshifted diagonals (for positive term of Sarrus' scheme, negative term
looks to be 0):

   11.2016    8.5445  -13.8587  -16.5158  -13.8587  -16.5158  -13.8587 
-16.5158

    0.00000  -12.53016    0.00000   12.53016    0.00000   12.53016   
0.00000   12.53016

   0   0   0   0   0   0   0   0

   0   0   0   0   0   0   0   0

   1.3286   3.9857   1.3286   3.9857   1.3286   3.9857   1.3286   3.9857

   0   0   0   0   0   0   0   0

   0   0   0   0   0   0   0   0

  -12.53016    0.00000   12.53016    0.00000   12.53016    0.00000  
12.53016    0.00000

And prod(x) of each diagonal:

   1.1477e+09  -0.0000e+00   0.0000e+00   0.0000e+00   7.8619e+02  
0.0000e+00   0.0000e+00  -0.0000e+00

>>

Obviously my determinant in this case  should be on the order of 1.0E+9
then. Why does this happens?
I'm not experienced much in FEM or CFD code in Octave, but i really wanted
to recreate some simulation from certain publication. How to solve huge set
of equations then? Maybe i'm missing something on maths level here? Do
people even inverse matrices in their finite elements codes then? It is easy
to get
+-Inf or 0 diagonal product if your matrix is large enough. Can (or should)
I increase the precision of computations somehow?



--
Sent from: http://octave.1599824.n4.nabble.com/Octave-General-f1599825.html