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?
