www.pudn.com > NumericalComputing.rar > tridisolve.m, change:2003-02-18,size:831b


function x = tridisolve(a,b,c,d) 
%   TRIDISOLVE  Solve tridiagonal system of equations. 
%     x = TRIDISOLVE(a,b,c,d) solves the system of linear equations 
%     b(1)*x(1) + c(1)*x(2) = d(1), 
%     a(j-1)*x(j-1) + b(j)*x(j) + c(j)*x(j+1) = d(j), j = 2:n-1, 
%     a(n-1)*x(n-1) + b(n)*x(n) = d(n). 
% 
%   The algorithm does not use pivoting, so the results might 
%   be inaccurate if abs(b) is much smaller than abs(a)+abs(c). 
%   More robust, but slower, alternatives with pivoting are: 
%     x = T\d where T = diag(a,-1) + diag(b,0) + diag(c,1) 
%     x = S\d where S = spdiags([[a; 0] b [0; c]],[-1 0 1],n,n) 
 
x = d; 
n = length(x); 
 
for j = 1:n-1 
   mu = a(j)/b(j); 
   b(j+1) = b(j+1) - mu*c(j); 
   x(j+1) = x(j+1) - mu*x(j); 
end 
 
x(n) = x(n)/b(n); 
for j = n-1:-1:1 
   x(j) = (x(j)-c(j)*x(j+1))/b(j); 
end