inverse Matrix problem, Newton's method

classic Classic list List threaded Threaded
3 messages Options
Reply | Threaded
Open this post in threaded view
|

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


Reply | Threaded
Open this post in threaded view
|

Re: inverse Matrix problem, Newton's method

Sebastian Schöps
mecht wrote
> 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)

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
problems?

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...

Sebastian



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


Reply | Threaded
Open this post in threaded view
|

Re: inverse Matrix problem, Newton's method

mecht
This post was updated on .
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.

Piotr

EDIT:
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 :<.

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