Tridiagonal matrix (thomas algorithm)
显示 更早的评论
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"
采纳的回答
John D'Errico
2011-3-11
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 个评论
Mehmet
2011-3-11
John D'Errico
2011-3-11
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
2018-5-2
Does MATLAB identifies it as tridiagonal matrix and solve it accordingly (Using optimized solver for tridiagonal matrix)?
John D'Errico
2018-5-2
编辑：John D'Errico
2018-5-2
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
2018-5-12
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
2018-5-12
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
2021-8-21
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.
更多回答（2 个）
Shantanu Vachhani
2015-12-24
编辑：Walter Roberson
2015-12-24
%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 个评论
John D'Errico
2018-5-12
编辑：John D'Errico
2018-5-12
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
2019-3-1
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
Troy
2024-10-26
0 个投票
Solve by using Thomas Method
1 个评论
Walter Roberson
2024-10-27
I do not understand how your Answer will help the original poster?
类别
在 帮助中心 和 File Exchange 中查找有关 Mathematics 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
