Pentadiagonal matrix to solve for Pressure variable

Hello everyone,

I want to solve my pressure equation implicitly by pentadiagonal matrix method. Here is the following equation.

---------------> (1)

Eq. (1) is the linear system of equaion with 5 pressure unknowns. For 3 x 4 = 12 cells, considering and as 1, the matrix can be determined as shown below, representing Pentadiagonal matrix.

Thus, in order to represent this pentagonal matrix in MATLAB, for-loop should be used.

% for internal cells excluding boundaries

for i=2:Ny-1

for j=2:Nx-1

U(i,j) = 1;

V(i,j) = -4;

W(i,j) = 1;

X(i,j) = 1;

Y(i,j) = 1;

end

end

% for left boundary

//

% for right boundary

//

% for top boundary

//

% for right boundary

//

In case of tri-diagonal matrix (1D problem), the algorithm is based on forward elimination and backward substitution to calculate the unknown variable, as shown below. Link: https://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm

%TriDiagonal_Matrix_Algorithm(N%, A#(), B#(), C#(), D#(), X#())

for i = 2:N

W = A(i) / B(i - 1)

B(i) = B(i) - W * C(i - 1)

D(i) = D(i) - W * D(i - 1)

end

for i = N-1:-1:1

X(i) = (D(i) - C(i) * X(i + 1)) / B(i)

end

for i = 1:N

X(N) = D(N) / B(N)

end

But, I'm not sure about the algorithm of pentadiagonal matrix. I have searched through forum links and documents, but I couldn't able to get the clarity of thoughts. Can someone guide me to find a way and kindly judge my above implementation.

Thank you

Answer 1

[John D'Errico]

WHY? Don't write code to do what was already written well by world class experts in linear algebra. You are using MATLAB. So use sparse matrices to represent a sparse matrix. A banded matrix is indeed sparse, and very nicely so.

help spdiags

 SPDIAGS Sparse matrix formed from diagonals.
    SPDIAGS, which generalizes the function "diag", deals with three
    matrices, in various combinations, as both input and output.
 
    [B,d] = SPDIAGS(A) extracts all nonzero diagonals from the m-by-n
    matrix A. B is a min(m,n)-by-p matrix whose columns are the p nonzero
    diagonals of A. d is a vector of length p whose integer components
    specify the diagonals in A.
 
    B = SPDIAGS(A,d) extracts the diagonals specified by d.
 
    A = SPDIAGS(B,d,A) replaces the diagonals of A specified by d with
        the columns of B. The output is sparse.
 
    A = SPDIAGS(B,d,m,n) creates an m-by-n sparse matrix from the
        columns of B and places them along the diagonals specified by d.
 
    Roughly, A, B and d are related by
        for k = 1:p
            B(:,k) = diag(A,d(k))
        end
 
    Example: These commands generate a sparse tridiagonal representation of
    the classic second difference operator on n points.
        e = ones(n,1);
        A = spdiags([e -2*e e], -1:1, n, n)
 
    Some elements of B, corresponding to positions "outside" of A, are not
    actually used. They are not referenced when B is an input and are set
    to zero when B is an output. See the documentation for an illustration
    of this behavior.
 
    See also DIAG, SPEYE.

    Documentation for spdiags
       doc spdiags

    Other uses of spdiags

       codistributed/spdiags    gpuArray/spdiags

Then use backslash for the solve.

B = repmat([1 1 -4 1 1],10,1);
M = spdiags(B,[-4 -1 0 1 4],10,10);
spy(M)

full(M)

ans = 10×10

-4 1 0 0 1 0 0 0 0 0 1 -4 1 0 0 1 0 0 0 0 0 1 -4 1 0 0 1 0 0 0 0 0 1 -4 1 0 0 1 0 0 1 0 0 1 -4 1 0 0 1 0 0 1 0 0 1 -4 1 0 0 1 0 0 1 0 0 1 -4 1 0 0 0 0 0 1 0 0 1 -4 1 0 0 0 0 0 1 0 0 1 -4 1 0 0 0 0 0 1 0 0 1 -4

Leave the matrix as a sparse matrix of course.

S = rand(10,1); % I don't have your data.
P = M\S;

[Kumaresh Kumaresh]

Thank you John for your answer. It is helpful.

However, if I need to analyse my problem with increase in grid numbers, it should be hard in creating the sparse matrix every time. In that case, sufficient algorithm based on for-loop will be helpful, right ? I can find that algorithm for TDMA (tridiagonal matrix), but not for Pentadiagonal matrix.

[Torsten]

However, if I need to analyse my problem with increase in grid numbers, it should be hard in creating the sparse matrix every time.

Is this hard ?

B = repmat([1 1 -4 1 1],10,1);
M = spdiags(B,[-4 -1 0 1 4],10,10);
S = rand(10,1); % I don't have your data.
P = M\S;

In case your systems are too big for sparse direct methods like above, use an iterative scheme. It should be convergent since your matrix is weakly diagonal-dominant.

[John D'Errico]

You have a viable solution to solve the problem. It will be very quick to solve such a problem, even for huge matrices. A quick test on my computer, for matrices of size 1e6 has the solve take 1/4 of a second, which includes building the matrix itself.

Why spend a fair amount of time to find the algorithm, then write and debug the code, to then need to maintain it? Good programming policy says not to waste programmer time on something you have that works, and works well. I'm not even at all sure your code would be faster than the simple sparse matrix solution. Lengthy loops can be slow.

[Kumaresh Kumaresh]

Thank you both of you.

I prefer using for loop for "Pentadiagonal matrix" rather than "repmat" and "spdiags" to understand my problem better.

clc; clear all;

dx = 1; dy = 1;

%U P(i-1,j) + V P(i,j) + W P(i+1,j) + X P(i,j-1) + Y P(i, j+1) = S

U = 1/dx^2;

V = -(2/dx^2)-(2/dy^2);

W = 1/dx^2;

X = 1/dy^2;

Y = 1/dy^2;

Nx = 3; Ny = 4;

A = zeros(Nx*Ny, Nx*Ny);

b = rand(Nx*Ny, 1);

%Constructing the matrix

for i = 1:(Nx*Ny)

A(i,i) = V;

for i = 1:(Nx*Ny-1)

A(i,i+1) = W;

A(i+1,i) = U;

end

for i = 1:(Nx*Ny-Ny)

A(i,i+Ny) = X;

A(i+Ny,i) = Y;

end

end

% for loop // Left, right, top and bottom boundaries

P = A\b;

[Torsten]

Your matrix A from above is wrong in many ways.

I suggest you try to solve

d^2u/dx^2 + d^2u/dy^2 = f

on a rectangle with Dirichlet boundary conditions for u to test and correct it.

[John D'Errico]

First, you apparently want to solve your problem efficiently. There is no reason why you would be asking to use a pentadiagonal solver, if it was not that.

But then you say you want to use a loop to construct your matrices, because you want to "understand" your problem better. Understanding how to use the right software to solve your problem means that you should use the proper tools.

If you REALLY, REALLY want to fully understand every aspect about your problem, you should write the code to multiply, add, subtract and divide numbers, from scratch. You should write EVERYTHING from scratch. Why are you using MATLAB anyway? Don't even use C. Go directly to assembly code. Perform those adds and multiplies using registers. Hey, I wrote code on machines where I needed to toggle in the instructions using dip switches. I learned nothing more than I hated writing code that way.

My point is, this is silly. You gain no more understanding of a problem by doing things the hard way. And you will gain no speed. Anyway, you already know how to construct the matrix you have decided to construct, doing it the inefficient way here with loops. So there is no understanding to be gained from that.

We are trying to teach you how to solve your problem, using MATLAB. If you want to understand something, then learn to use MATLAB.

[Kumaresh Kumaresh]

Hello Torsten,

I've checked my A matrix (created for Pentadiagonal matrix). In my case, top and left have Dirichlet boundaries, whereas bottom and right have Neumann boundaries. I have got my pressure value, that looks correct. If posssible, could you please ellaborate, why my A matrix is wrong in many ways ?

Thank you Torsten once again.

Thank you John for your concern over my question. I don't understand your exaggerated response. I believe everyone here doesn't know everything. To answer one of your point, my answers and question are never deviated from the actual topic. You suggested me to use "repmat" and "spdiags" for pentadiagonal matrix, whereas I preferred "for-loop" for the same pentadiagonal matrix. That's it.