How to fast solve a system with block tridiagonal matrix

Hi everyone, I am trying to solve a system of pdes and I am a new one at Matlab. I used Euler implicit discretization method and it leads to solving a system of linear equations with block tridiagonal matrix in every iterations (respect to time discretization). It converges actually and runs normally in the small dimension of matrix but when I increase the dimension of matrix (100x100), it becomes extremely consuming time (over a week). So could you please advice me how to optimize the code. Here it is the piece of my code, that I should optimize

    %calculate matrix elements
    dd = (delta_z*delta_z/delta_t)*(2.*r_next - 3.*r_backup + delta_t*beta_tmp);
    aa(2:J) = r_backup(2:J)./dd(2:J);
    bb(2:J) = -(2*(2.*r_backup(2:J)- r_backup(1:J-1)) + (delta_z*delta_z/delta_t).*r_backup(2:J))./dd(2:J);
    cc(2:J) = (3*r_backup(2:J)- 2*r_backup(1:J-1))./dd(2:J);
    bb(J) = bb(J) + cc(J);
    %
    % creat matrice 
     A = diag(bb(2:J))+diag(aa(3:J),-1)+ diag(cc(2:(J-1)),1); % A has dimenssion (J-1)
     % vector y in the right side of equation Ac_next=y
     y(1)=c_backup(2)-aa(2)*c_backup(1);
     y(2:J-1)=c_backup(3:J);
     %
     %solving the system of equation
     c_next(2:J)=A\y';

Here c_next(2:J) are variables. The matrix elements can be solved in every iteration. The two steps are time consuming is the step to define matrix elements (but I think there is no way to optimize it any more), rest is only the step to solve the equation Ax=y. So if you have any idea, please advice me.

Thank you very much and I am looking forward to hearing from you. Cheers,

Answer 1

[Alan Weiss]

I wonder whether you have specified your matrix elements as sparse. If not, then I think that you should try using sparse matrices. See the documentation, such as Computational Advantages of Sparse Matrices.

Alan Weiss

MATLAB mathematical toolbox documentation

[Huyen Nguyen]

Thank you Alan, I am trying with using sparse matrix

Answer 2

[John D'Errico]

Alan is certainly correct, in that your matrix is not a sparse one. Were it sparse, then the solve (using backslash) would be quite fast.

So create a sparse matrix. A very good way to do so is to use my blktridiag , as found on the file exchange for download. You give it a list of blocks, and the matrix is created (and this will be fast too) and it will be in sparse form.

[Huyen Nguyen]

Thank you John, I was trying with your function, now solving systems Ac_next=y does not take long time like before but running the function creating sparse A really takes time (40% of the total time of program) :(. Maybe because my program requires lots of iterations :(