Tridiagonal matrix (thomas algorithm)

hi.

i want to solve a second order, homogeneous differential equation by using tridiagonal matrix. can u help me?

example: http://img846.imageshack.us/i/diffeq.jpg/

so, i want only a general matlab code to solve like this equation.

because i am using "finite difference method"

Answer 1

[John D'Errico]

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.

[Mehmet]

i am very new on matlab and didnt understood anything from your answer. But thank for your answer.i dont know anything about backslash.

i need only general example code (a basic application) like my differential equation.this example should solve second order, homogeneous diff eq. i think it was been solved by using "thomas algorithm" but if there another solution you can give me as .m file. i need this solution immediately.

thank for all John

[John D'Errico]

I gave you the solution as m-code. If you don't understand it, then it is high time to start reading the tutorials. Read the help for slash, or for linsolve, and for spdiags.

[Royi Avital]

Does MATLAB identifies it as tridiagonal matrix and solve it accordingly (Using optimized solver for tridiagonal matrix)?

[John D'Errico]

No. MATLAB does not care that it is explicitly a tridiagonal matrix. However, because it IS a tridiagonal sparse matrix, AND because the sparse solver is efficient on sparse matrices, MATLAB effectively does use an extremely efficient scheme to solve the problem.

Can you do better? Perhaps. Sometimes since the time I originally answered this question, the decomposition utility was introduced. As well, there are tools for bandwith computation.

N = 100000;
T = spdiags(rand(N,3),[-1:1],N,N) + 2*speye(N);
[Lband,Uband] = bandwidth(T)
Lband =
     1
Uband =
     1

So MATLAB can determine the matrix is banded, and in this case, is tridiagonal.

rhs = rand(N,1);
timeit(@() T\rhs)
ans =
    0.0024706

So the simple backslash is fast as hell. Do you really care if it uses the Thomas algorithm? Can we compute an LU decomposition?

timeit(@() lu(T))
ans =
     0.029801

Actually, it is way slower to compute the LU. The decomposition function is faster. Although we need to be careful. When I created T, I chose it very carefully, so that T is very well conditioned, and there should be no desire to do any pivoting. If the linear solve needed to do some pivoting here, those tools will potentially slow down.

timeit(@() decomposition(T,'banded'))
ans =
    0.0027719

If you can write MATLAB code for the Thomas algorithm solution to run on a 100K tridiagonal matrix in 0.0024 second or less, then well done.

The point is, simple backslash is indeed fast. Does it explicitly use the Thomas algorithm? Again, why would you care? In fact, the Thomas algorithm can be unstable, IF pivoting were necessary, whereas backslash will be far better. (Do I really need to prove this, and show how the Thomas algorithm can be problematic?)

[Royi Avital]

John, Thanks for this great comment. So you say the general strategy to solve Sparse Matrices is so efficient it will do well with Tridiagonal without being explicitly aware it is Tridiagonal Matirx. Yet, what about the cases the matrix isn't Sparse? Is it fast when one transfer it into Sparse form (Namely taking into account the time converting the matrix into sparse form)?

[John D'Errico]

Are you asking, if someone for some incredibly silly reason generated a full matrix, that was actually tridiagonal, but huge?

May I ask why in the name of god and little green apples would you ever create a large, full tridiagonal matrix in the first place? If you will create a huge trigiagonal matrix, then create it as sparse to start. Anything you will do with it, including store it, will be faster, more efficient, etc. Just creating that matrix as a full matrix in the first place is inefficient.

What you don't want to do is something silly like create a full huge matrix, that should actually be sparse. Then convert it to sparse form just to solve a system, and then discard the sparse form. Yes, it takes some time to convert between full and huge sparse forms. For example:

N = 10000;
tic,Afull = diag(3*ones(N,1),0) + diag(ones(N-1,1),-1) + diag(ones(N-1,1),1);toc

Elapsed time is 2.750122 seconds.

So it took 2.75 seconds to create that 10K sparse tridiagonal matrix.

Only about 0.2 seconds to convert to sparse form.

timeit(@() sparse(Afull))
ans =
      0.20962

Now, computing the solution of the sparse system?

Asparse = sparse(Afull);
rhs = rand(N,1);
timeit(@() Asparse\rhs)
ans =
   0.00025803

So a quarter of a millisecond to solve it, in sparse form.

By way of comparison, using the looped Thomas algorithm code posted by Shantanu Vachhani as an answer here? After I removed the input statements, and put it into a function form, that code was 10 times slower than the sparse backslash.

As for creating the tridiagonal matrix directly as a sparse matrix?

timeit(@() spdiags(ones(N,3),-1:1,N,N) + 2*speye(N))
ans =
    0.0036851

It took close to 750 times as long to create it as full as it did to create the sparse matrix directly.

The point is, get used to using sparse form when you have a sparse matrix, and use the methods in MATLAB to work with your matrices. Don't create the matrix as a full one, then use methods that should be reserved for student homework assignments.

If you had to multiply two matrices in MATLAB, would you write it as A*B, or would you write the triply nested for loop form that we learned to write in the dark days of old Fortran? Use the tools provided in MATLAB. Learn to use them well, and your code will gain greatly.

[Guojun Chen]

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.

Answer 2

[Shantanu Vachhani]

%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

[John D'Errico]

For someone who wants to use this code, remember that it will be very much slower than just using backslash. Looped student code is rarely efficient code. At best, it was an answer to a homework assignment. But no more than that.

For example, on a quick test with a 10k by 10k tridiagonal matrix, this looped code was roughly 10 times lower than just using backslash properly.

[Mohammad Gohardoust]

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

Answer 3

[Troy]

Solve by using Thomas Method

[Walter Roberson]

I do not understand how your Answer will help the original poster?