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