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

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

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

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);
}