Tridiagonal matrix (thomas algorithm)

415 views (last 30 days)
Mehmet
Mehmet on 11 Mar 2011
Commented: Guojun Chen on 21 Aug 2021
hi.
i want to solve a second order, homogeneous differential equation by using tridiagonal matrix. can u help me?
so, i want only a general matlab code to solve like this equation.
because i am using "finite difference method"

Accepted Answer

John D'Errico
John D'Errico on 11 Mar 2011
Why not just build it as a sparse matrix using spdiags, then solve using backslash? It will be quite fast for a tridiagonal matrix, and you won't need to write any solver at all.
For example, I won't bother to do more than create a random tridiagonal matrix, rather than building one directly from your equation, but the time is all that matters.
n = 100000;
A = spdiags(rand(n,3),-1:1,n,n);
b = rand(n,1);
tic,x = A\b;toc
Elapsed time is 0.023090 seconds.
So to solve a system with 1e5 unknowns took me .023 seconds, with virtually NO time invested to write anything.
  7 Comments
Guojun Chen
Guojun Chen on 21 Aug 2021
Well, I am a person who used Thmos algorithm a lot. Most of the time I still use the time-series routine instead of the sparse backslash. The reason is, when my ODE system is highly nonlinear, I need to get the three diagonals first as three vectors. Then after I have them, I am comparing the speed of: (1) create a sparse matrix based on them + use sparse backslash; (2) solve with time-series Thomas algorithm. (2) is usually faster than (1) for me.

Sign in to comment.

More Answers (1)

Shantanu Vachhani
Shantanu Vachhani on 24 Dec 2015
Edited: Walter Roberson on 24 Dec 2015
%of the form AX=B
n=input('enter the order for the matrix');
for(i=1:n)
for(j=1:n)
a(i,j)=input('enter the element of coefficient matrix');
end
end
for i=1:n
r(i)=input('enter the RHS');
end
R(1)=0;
P=zeros(1,n);
Q=zeros(1,n-1);
R=zeros(1,n);
Y=zeros(1,n-1);
for i=1:n
P(i)=a(i,i);
end
for i=1:n-1
Q(i)=a(i,i+1);
end
for i=1:n-1
R(i+1)=a(i+1,i);
end
Y(1)=Q(1)/P(1);
for i=2:n-1
Y(i)=Q(i)/(P(i)-R(i)*Y(i-1));
end
W(1)=r(1)/P(1);
for i=2:n
W(i)=(r(i)-R(i)*W(i-1))/(P(i)-R(i)*Y(i-1));
end
x(n)=W(n);
for i=n-1:-1:1
x(i)=W(i)-Y(i)*x(i+1);
end
  2 Comments
Mohammad Gohardoust
Mohammad Gohardoust on 1 Mar 2019
Thanks John for your complete answers in this page. In the case of tridiagonal matrix, I have tried what you have suggested and also tested the Thomas algorithm I have implemented. The results were comparable and even a bit to the favor of Thomas algorithm.
function h = Thomas(ld,md,ud,a)
% Solves linear algebraic equation where the coefficient matrix is
% tridiagonal. ld, md and ud stands for lower-, main, and upper-
% diagonal respectively. a is the answer matrix and h is the solution.
N = length(md) ;
w = zeros(N, 1) ; g = zeros(N, 1) ;
w(1) = ud(1)/md(1) ; g(1) = a(1)/md(1) ;
if isrow(ud)
ud = ud' ;
end
if isrow(ld)
ld = ld' ;
end
ud = [ud; 0] ; ld = [0; ld] ;
for i=2:N
w(i) = ud(i)/(md(i)-ld(i)*w(i-1)) ;
g(i) = (a(i)-ld(i)*g(i-1))/(md(i)-ld(i)*w(i-1)) ;
end
h = zeros(N, 1) ;
h(N) = g(N) ;
for i=N-1:-1:1
h(i) = -w(i)*h(i+1)+g(i) ;
end
end

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!