I am stuck in making the pentadiagonal matrix A

a0=[1 (1+((tm+tp).*LBx(2:nx)))+(tm_y+tp_y).*LBy(2:ny) 1]; 
S=zeros(nx+1,ny+1);
s=diag(a0)+S; 
%%%%% Lower diagonal matrix b%%%

b=[0 -tm.*LBx(1:nx-1) 0];
b1=b(2:end);
XX=diag(b1,-1)+S;
%%%% upper diagonal matrix c%%%%
c=[0 -tp.*LBx(3:nx+1) 0];
w=c(1:end-1);
W=diag(w,1)+S;
%%%% double lower diagonal matrix%%%%
d=[0 -tm_y.*LBx(1:ny-1) 0];
dd=d(3:end);
D=diag(dd,-2)+S;
%%%%%%%% double upper diagonal matrix%%%%
e=[0 -tp_y.*LBx(3:ny+1) 0];
r=e(1:end-2);
E=diag(r,2)+S;
%%%%%%%%%Right hand side of the matrix%%%%%%%%%
f0=[0 b_n(2:nx)+(2.*(dt./(dx.^2))).*(S_iter(2:nx).*LBx(2:nx))+(2.*(dt./(dy.^2))).*(S_iter(2:ny).*LBy(2:ny))-(2.*(dt./(dx.^2))).*Q_iter_x(2:nx)-(2.*(dt./(dy.^2))).*Q_iter_y(2:ny)-(dt/(dx^2)).*(LBx(1:nx-1)).*(S_iter(1:nx-1))+(dt/dx^2).*Q_iter_x(1:nx-1)-(dt/(dx^2)).*(LBx(3:nx+1)).*(S_iter(3:nx+1))+(dt/dx^2).*Q_iter_x(3:nx+1)-(dt/(dx^2)).*(LBy(1:ny-1)).*(S_iter(1:ny-1))+(dt/dy^2).*Q_iter_y(1:ny-1)-(dt/(dy^2)).*(LBy(3:ny+1)).*(S_iter(3:ny+1))+(dt/dx^2).*Q_iter_y(3:ny+1) 0];

%%% matrix calculation%%%%
A=zeros((nx+1)*(ny+1),(nx+1)*(ny+1));
f=zeros(1,(nx+1)*(ny+1));
for ind=1:(nx+1)*(ny+1)
    
    %%%%%%%This line uses the quorem function to compute the quotient and remainder of ind divided by (nx+1)*(ny+1). The function returns two outputs: k for the quotient and j for the remainder. The sym function is used to convert the inputs to symbolic variables.
    %%'sym' is a function that converts a numerical value to a symbolic value.
    [k,j]=quorem(sym(ind),sym((nx+1).*(ny+1)));
    %%%This line checks whether j is equal to 0 or 1, or whether k is equal to 0 or ny. If any of these conditions are true, the code inside the if block will be executed.   
    if j==0||j==1||k==0||k==ny   %||j==nx why ??
        
        %%This line sets the ind-th diagonal element of the matrix A to 1. This effectively enforces a Dirichlet boundary condition at the node corresponding to index ind.
        
        A(ind,ind)=1;
        f(ind)=0;
    else
        s_vec = reshape(s, [], 1); % Reshape the matrix s into a column vector
        A = A+spdiags(s_vec, 0, (nx+1)*(ny+1), (nx+1)*(ny+1)); % Add s to the main diagonal of A
        XX_vec = reshape(XX, [], 1);
        A = A + spdiags(XX_vec, -1, (nx+1)*(ny+1), (nx+1)*(ny+1)); % Add XX to the subdiagonal of A
        W_vec = reshape(W, [], 1);
        A = A + spdiags(W_vec, 1, (nx+1)*(ny+1), (nx+1)*(ny+1)); % Add W to the superdiagonal of A
        D_vec = reshape(D, [], 1);
        A = A + spdiags(D_vec, -2, (nx+1)*(ny+1), (nx+1)*(ny+1)); % Add D to the sub-subdiagonal of A
        E_vec = reshape(E, [], 1);
        A = A + spdiags(E_vec, 2, (nx+1)*(ny+1), (nx+1)*(ny+1)); % Add E to the super-superdiagonal of A
        A(ind,ind) = a0(1);
        A(ind,ind-1) = a0(2);
        A(ind,ind+1) = a0(3);
        A(ind,ind-(nx+1)) = b(2);
        A(ind,ind+(nx+1)) = c(2);
        A(ind,ind-2*(nx+1)) = d(2);
        A(ind,ind+2*(nx+1)) = e(2);
        f(ind) = f0;
    end
end
end

[Dyuman Joshi]

Stuck how exactly? Are you not getting the desired output? If so, then provide the output for the given input.

Also, there are undefined variables in your code, we can't run your code.

Answer 1

[John D'Errico]

This sort of thing is a BAD idea:

A = A+spdiags(s_vec, 0, (nx+1)*(ny+1), (nx+1)*(ny+1)); % Add s to the main diagonal of A
XX_vec = reshape(XX, [], 1);
A = A + spdiags(XX_vec, -1, (nx+1)*(ny+1), (nx+1)*(ny+1)); % Add XX to the subdiagonal of A
W_vec = reshape(W, [], 1);
A = A + spdiags(W_vec, 1, (nx+1)*(ny+1), (nx+1)*(ny+1)); % Add W to the superdiagonal of A
D_vec = reshape(D, [], 1);
A = A + spdiags(D_vec, -2, (nx+1)*(ny+1), (nx+1)*(ny+1)); % Add D to the sub-subdiagonal of A
E_vec = reshape(E, [], 1);
A = A + spdiags(E_vec, 2, (nx+1)*(ny+1), (nx+1)*(ny+1)); % Add E to the super-superdiagonal of A

That is, you create a sparse matrix with ONE diagonal, then you add to it, ANOTHER spade matrix with ONE diagonal, etc. Then you do it 5 times.

DON'T DO THAT!!!!!!!!

DON'T DO THAT!!!!!!!!

DON'T DO THAT!!!!!!!!

DON'T DO THAT!!!!!!!!

DON'T DO THAT!!!!!!!!

Call spdiags ONCE, creating the entire matrix with all 5 diagonals in ONE call.

The problem is, when you create one diagonal at a time, is then you force MATLAB to constantly re-sort the elements of the array. It is far less efficient to generate the matrix that way.

Yes, people do it that way with FULL matrices using diag. That is ok, as long as the matrices are full. It is a bad idea with sparse matrices.

LEARN TO USE SPDIAGS, PROPERLY.