I am stuck in making the pentadiagonal matrix A
6 views (last 30 days)
Show older comments
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
1 Comment
Dyuman Joshi
on 9 Apr 2023
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.
Answers (1)
John D'Errico
on 18 May 2023
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.
0 Comments
See Also
Categories
Find more on Operating on Diagonal Matrices in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!