Diagonal matrices with spdiags
11 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Matthew Hunt
am 25 Okt. 2019
Kommentiert: Matthew Hunt
am 25 Okt. 2019
I'm working on a numerical solution to an equation and as part of this I have to solve a matrix solution. The system of equations in a tridiagonal matrix I have been informed that there is a routine called spdiags which allows me access to specialised solution/inversing routines which should speed up my code.
The code I use is:
s=0.12;
N_r=30;
r=linspace(0,1,N_r)';
dr=r(2);
r_plus=r+0.5*dr;
r_minus=r-0.5*dr;
a_plus=s*r_plus(1:end-1).^2;
a_minus=s*r_minus(1:end-1).^2;
a=-(r.^2+s*(r_plus.^2+r_minus.^2));
A=diag(a_plus,1)+diag(a)+diag(a_minus,-1);
A(1,1)=-1;A(1,2)=1;
A(N_r,N_r-1)=s*(r_plus(N_r)^2+r_minus(N_r)^2);
This provides the matrix that I want. I can run the code and it's pretty fast but I want to see that if I define the A matrix as a spdiags matrix:
B_plus=s*r_plus.^2;
B_minus=s*r_minus.^2;
B=spdiags([B_minus a B_plus],-1:1,N_r,N_r);
B(1,1)=-1;B(1,2)=1;
B(N_r,N_r-1)=s*(r_plus(N_r)^2+r_minus(N_r)^2);
Now hopefully, these should yield the same matrix, but they don't. What am I doing wrong?
0 Kommentare
Akzeptierte Antwort
Matt J
am 25 Okt. 2019
Bearbeitet: Matt J
am 25 Okt. 2019
s=0.12;
N_r=30;
r=linspace(0,1,N_r)';
dr=r(2);
r_plus=r+0.5*dr;
r_minus=r-0.5*dr;
a_plus=s*r_plus(1:end-1).^2;
a_minus=s*r_minus(1:end-1).^2;
a=-(r.^2+s*(r_plus.^2+r_minus.^2));
D=[[a_minus(:);0], a(:), [0;a_plus(:)]]; %<---changed
B=spdiags(D,-1:1,N_r,N_r); %<--changed
B(1,1)=-1;B(1,2)=1;
B(N_r,N_r-1)=s*(r_plus(N_r)^2+r_minus(N_r)^2);
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu Creating and Concatenating Matrices finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!