# Diagonal matrices with spdiags

6 views (last 30 days)
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?
##### 0 CommentsShowHide -1 older comments

Sign in to comment.

### Accepted Answer

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 CommentsShowHide 4 older comments
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.

### Community Treasure Hunt

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

Start Hunting!