Решение линейных систем с трехдиагональной матрицей
Страница 2. Алгоритм Томаса

Алгоритм Томаса

Поиск решения производится в виде:

u_i-1=bet_i*u_i+gam_i

В этом случае для bet и gam получаются рекуррентные формулы:

bet_i+1=c_i*fi, gam_i+1=(f_i+a_i*gam_i)*fi,     fi=1/(b_i-a_i*bet_i);     bet₀=gam₀=0.

По этим формулам вспомогательные массивы вычисляются для i от 1 до N, затем находится u_N-1=gam_N, а затем при известных вспомогательных массивах вычисляются все значения u.

Алгоритм принципиально не может иметь ведущего элемента, и есть вероятность того, что он не сойдется даже для несингулярной матрицы. Для этого в программе, реализующей прогонку, необходимо контролировать значения fi. Еще два замечания: длины рекуррентных цепочек порядка N, поэтому при экспоненциальном накоплении погрешностей (что происходит при преобладании значений |bet_i|>1) прогонка практически обязательно разойдется. Второе: алгоритм сохраняет исходные значения коэффициентов. Доказано, что для устойчивой работы алгоритма Томаса достаточно диагонального преобладания, однако, во многих случаях он сходится и при отсутствии такового.

/* 
 Алгоритм Томаса, решающий трехдиагональную линейную систему. 
Описание:
 int Sweep(double *u,double *a,double *b,double *c,double *f,
 int n,double *bet,double *gam);
Параметры:
 u - разрешение (размера n), на выходе;
 a,b,c,f - массивы коэффициентов (как в алгоритмах, описанных выше, размера n)
 n - размер разрешения и массива коэффициентов;
 bet, gam - дополнительные массивы (размера n+1).
 Результат: 
 0 - сбой алгоритма,
 1 - положительный результат.
*/

/* Алгоритм Томаса */
int Sweep(double *u,double *a,double *b,double *c,double *f,int n,
 double *bet,double *gam) {
 register int n;
 double fi;
 /* Вычисление дополнительных массивов (умножение системы) */
 bet[0]=gam[0]=0.;
 for(i=0;i<n;i++) {
 fi=b[i];
 if(i>0) fi-=a[i]*bet[i];
 if(fi==0.) return(0);
 fi=1./fi; 
 if(i<n-1) bet[i+1]=c[i]*fi; 
 gam[i+1]=(f[i]+a[i]*gam[i])*fi;
 }
 /* обратная замена */
 u[n-1]=gam[n];
 for(i=n-1;i>0;i--) u[i-1]=bet[i]*u[i]+gam[i];
 return(1);
}