Diagonal matrices with spdiags

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?

Answer 1

[Matt J]

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);

[Matthew Hunt]

What was I doing wrong?

[Matt J]

I've marked the lines in the code that changed. Your B_plus didn't pre-pad a_plus with anything.

[Matthew Hunt]

Pre-pad?

[Matt J]

Notice how I specified the upper diagonal with a_plus pre-padded by a zero:

[0;a_plus(:)]

and conversely for a_minus.

The behavior is different depending on whether have you have more rows than columns or vice versa

https://www.mathworks.com/help/matlab/ref/spdiags.html#mw_6de431d3-8887-4b16-b947-86a610ddea8c

[Matthew Hunt]

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.