I'm writing some Octave code to solve a finite element-like problem using the
My Jacobian matrix looks like this (spy(J)):
<img src="https://preview.ibb.co/hknY5A/spy.jpg" alt="spy" border="0" />
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
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.
I tried my code with minimal number of finite elements 2x2x2, which gives
Jacobian on 8x8 size, and the result was not that low.
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
+-Inf or 0 diagonal product if your matrix is large enough. Can (or should)
I increase the precision of computations somehow?
> 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
Are you sure that your boundary conditions are correct? In particular, may
there a Dirichlet condition be missing? A missing boundary condition is the
most common reason for a FEM matrix to be singular. Does it work for linear
In either case, do not use "inv", "rank" or "det" to solve, check the
regularity or rank of your matrix. Those algorithms are not reliable for
large problems. Just use Gauss or lu but make sure to use a sparse matrix
decomposition. Indeed, the most convenient choice is using the backslash
("x=A\b"). It will automatically take the best approach. If it fails then
almost surely you have made a mistake...
Thank you for your answer.
Boundary conditions might be okay, but today i noticed that my whole concept
is wrong. If diagonal is the negative sum of other elements you'll have zero
column if you add up them to one (while doing gaussian elimination), so it
seems my matrix is actually truely singular (i haven't noticed it at first,
because it got distorted with page formatting). I will have to reformulate
my whole problem once again.
So actually my code fixes one value so matrix is A=jacobian_mat(2:end,2:end) and matrix is no singular (my boundary condition ). It just does not converge at any conditions and the nature of my problem does not belong to this forum :<.