Hi,

ich habe eine Frage zum SIP. Bei Wikipedia wird der Algorithmus erklärt (http://en.wikipedia.org/wiki/Stone_method). Leider geht für mich aus der Beschreibung nicht hervor, wie die Matrix N bestimmt wird.

Ax = (M-N)x = (LU-N)x = b

what did you say!

Hi,

I don't understand the SIP (http://en.wikipedia.org/wiki/Stone_method). One reason is, that there is no explanation what does N mean.

Ax = (M-N)x = (LU-N)x = b

If A and b are known where is the reason to go on with an iterative process? The L-U-decompensation leads to the exact result....

"If A and b are known where is the reason to go on with an iterative process? The L-U-decompensation leads to the exact result...."

Because A is generally a very large (sparse) matrix whose LU factorization is dense => very large amount of memory to store ; e.g. if you have 1 million grid points then in double precision you would reguire around 7450Gb of RAM just to store the complete LU factoriztion. For a Laplacian type stencil (7 nonzero entries per row in A) then the ILU will only require around 54Mb (and significantly less if the matrix entries are constant)*. Also, in gerneral, the number of floating point operations will be lower for a suitably accelerated iterative scheme with a good initial guess compared to the full LU (assuming you could store it).

(*) SIP will actually require the full 54Mb even when the entries of A are constant unlike ILU due to extra interpolations performed in SIP.

NOTE: SIP only works for matrices (A) obtained via a discretization of a PDE unlike ILU.

"Because A is generally a very large (sparse) matrix whose LU factorization is dense => very large amount of memory to store ; e.g. if you have 1 million grid points then in double precision you would reguire around 7450Gb of RAM just to store the complete LU factoriztion. For a Laplacian type stencil (7 nonzero entries per row in A) then the ILU will only require around 54Mb (and significantly less if the matrix entries are constant)*. Also, in gerneral, the number of floating point operations will be lower for a suitably accelerated iterative scheme with a good initial guess compared to the full LU (assuming you could store it)."

I know that A is a sparse matrix. I thought the L-U algorithm could be used in this way:

A*u = r and A=L*U

L*q = r --> U*u=q (forward/backward substitution)

This would lead to the exact solution u if the right hand side q and A are known. I cannot see why iterations would be necessary...

How did you calculate the required mem for 1.000.000 grid points?

Thx.

As you have already written: Ax = (M-N)x = (LU-N)x = b

so A=(LU-N) and not A=LU ,

LU is not the complete LU decomposition, meaning it has only nonzero elements where A has nonzero elements. It is sparse, while the full LU decomposition wouldn't be so sparse. The other elements are stored in the matrix N which of course could be calculated by N=LU-A. This decomposition is used in iterative methods: Mx = Nx + b can be solved iteratively by: M x_i+1 = N x_i + b or x_i+1 = M^(-1) N x_i + M^(-1)b This is useful If M can easily be inverted, and M^(-1)N has specific properties (abs of largest eigenvalue <1).

"How did you calculate the required mem for 1.000.000 grid points?"

It's what would be required to store a full 10^6 x 10^6 matrix which you would have to consider doing if you formed the full LU factorization.

As an exercise write down the finite difference equations for the Laplace equation on a square with 10^2 and 25^2 points and form the LU factorization of the matrix and see how many zeroes there are and compare this to the original matrix which has just 5 bands.

Could you give an explicit example for N. I read several papers about the L-U-decompensation and I got mixed up with it. I got the impression the word L-U-decompensation is interpreted in different ways.

[2.0 1.0 0.0 0.0 0.0 0.0 0.0]*****[1.0 0.0 0.0 0.0 0.0 0.0 0.0]*****[2.0 1.0 0.0 0.0 0.0 0.0 0.0] [1.0 2.0 1.0 0.0 0.0 0.0 0.0]***** [0.5 1.0 0.0 0.0 0.0 0.0 0.0]***** [0.0 1.5 1.0 0.0 0.0 0.0 0.0] [0.0 1.0 2.0 1.0 0.0 0.0 0.0]*****[0.0 0.667 1 0.0 0.0 0.0 0.0]*****[0.0 0.0 1.3 1.0 0.0 0.0 0.0] [0.0 0.0 1.0 2.0 1.0 0.0 0.0]**=**[0.0 0.0 0.749 1 0.0 0.0 0.0]*****[0.0 0.0 0.0 1.25 1 0.0 0.0] [0.0 0.0 0.0 1.0 2.0 1.0 0.0]*****[0.0 0.0 0.0 0.8 1.0 0.0 0.0]*****[0.0 0.0 0.0 0.0 1.2 1.0 0.0] [0.0 0.0 0.0 0.0 1.0 2.0 1.0]*****[0.0 0.0 0.0 0.0 0.833 1 0.0]*****[0.0 0.0 0.0 0.0 0.0 1.166 1] [0.0 0.0 0.0 0.0 0.0 1.0 2.0]*****[0.0 0.0 0.0 0.0 0.0 0.857 1]*****[0.0 0.0 0.0 0.0 0.0 0.0 1.14] *****A*****=**********L*********************U One paper calls this a L-U-decompensation....which would not be right with your definition!

I think I understood, why it makes sense to solve the system of equations iteratively. Using SIPSOL for non-linear equation systems means that A is not constant because the coefficients include mass-fluxes (rho*u)!

The LU in SIP is not the LU decomposition it is the Incomplete LU (or ILU see the book by Peric & Ferziger for an explanation). The full LU has N=0 but because ILU is not the complete factorization there is an error (N = LU-A is nonzero).