Diagonal matrices with spdiags

6 views (last 30 days)
Matthew Hunt
Matthew Hunt on 25 Oct 2019
Commented: Matthew Hunt on 25 Oct 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?

Accepted Answer

Matt J
Matt J on 25 Oct 2019
Edited: Matt J on 25 Oct 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);
  5 Comments
Matthew Hunt
Matthew Hunt on 25 Oct 2019
So I tried out your code, and it changes the solution of my equation which is not something I expected. It sped up the code by a factor of 2 which isn't bad. However, it completely changed the solution of the equation. So I'm not too sure what's going on.

Sign in to comment.

More Answers (0)

Tags

Community Treasure Hunt

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

Start Hunting!