I am stuck in making the pentadiagonal matrix A

9 Ansichten (letzte 30 Tage)
saad raza
saad raza am 8 Apr. 2023
Beantwortet: John D'Errico am 18 Mai 2023
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 Kommentar
Dyuman Joshi
Dyuman Joshi am 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.

Melden Sie sich an, um zu kommentieren.

Antworten (1)

John D'Errico
John D'Errico am 18 Mai 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.

Kategorien

Mehr zu Operating on Diagonal Matrices finden Sie in Help Center und File Exchange

Tags

Produkte

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by