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?
Best Answer