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

На этой странице представлены два альтернативных друг другу алгоритма решения линейных систем с трехдиагональной матрицей. Подобные системы возникают при реализации решений параболических или эллиптических уравнений в одномерном случае либо в многомерном при покоординатном расщеплении решения.
Первый из этих алгоритмов использует факторизацию системы, называется прогонкой (или алгоритмом Томаса) и входит во все классические учебники по численным методам. Второй алгоритм (редукция) использует последовательное бинарное исключение уравнений из системы; он менее известен и более сложен, однако обладает очень высокой устойчивостью и меньшим накоплением погрешности в процессе счета.

Трехдиагональная линейная система уравнений может быть записана в следующем виде:

a_i*u_i-1-b_i*u_i+c_i*u_i+1+f_i=0, где 0<=i<N, a₀=c_N-1=0, u - искомое решение.

Матрица этой системы состоит из главной диагонали (с коэффициентами -b_i) и примыкающих к ней сверху и снизу диагоналей (соответственно c_i и a_i).

Прямое (безытерационное) решение этой системы уравнений можно проводить методом ее факторизации (прогонка, или алгоритм Томаса), либо методом последовательного исключения сначала нечетных, затем делящихся на 2, но не на 4, и т.д. уравнений и соответствующего преобразования коэффициентов. Последний метод называется редукцией или бинарным исключением.

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

Решение ищется в виде:

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) прогонка практически обязательно разойдется. Второе: алгоритм сохраняет исходные значения коэффициентов. Доказано, что для устойчивой работы алгоритма Томаса достаточно диагонального преобладания, однако, во многих случаях он сходится и при отсутствии такового.

Редукция

Общее описание алгоритма таково.
Если провести исключение из имеющейся системы уравнений всех уравнений, находящихся на нечетных позициях, то от системы N исходных уравнений мы перейдем к системе из (N-1)/2+1 уравнений только для четных значений u, нечетные же значения будут вычисляться через четные. Повторяя эту процедуру для укороченной системы, мы в конце концов придем к единственному уравнению для узла, лежащего в позиции (N-1)/2 (только если N-1 - степень двойки). Затем, проводя обратную подстановку, определяются крайние, а затем и все прочие узлы.

Формулы для пересчета коэффициентов на каждом этапе исключения таковы:

a_i = a_i-1*A,
c_i = c_i+1*C,
b_i = b_i - (a_i+1*C + c_i-1*A),
f_i = f_i + f_i-1*A + f_i+1*C;
где A=a_i/b_i-1, C=c_i/b_i+1 -- естественно, до пересчета. Формулы для обратной подстановки: u_i = (a_i*u_i+1 + c_i*u_i-1 + f_i)/b_i.

Если число уравнений не является степенью двойки, алгоритм реализуется методом фиктивного восполнения области определения до степени 2. При этом никаких реальных действий, кроме нескольких целочисленных проверок, не производится, а только коэффициенты, индекс которых выпадает из области определения [0,N-1] включительно, считаются равными нулю и в преобразованиях не участвуют (подробности в программе).

Особая устойчивость этого алгоритма по сравнению с алгоритмом Томаса связана с тем, что длина рекуррентных цепочек здесь не превышает log2(N)+1, и в связи с этим даже при экспоненциальном росте погрешности исходный результат будет ограниченным. Кроме того, в алгоритме Томаса в рекуррентных цепочках производятся операции между коэффициентами соседних узлов, при редукции на этапах больших начального - между далеко отстоящими, что уменьшает вероятность образования погрешности при вычитании двух близких чисел. Проверено, что алгоритм редукции успешно работает во многих случаях, когда алгоритм Томаса (прогонка) расходится.

Алгоритм не требует дополнительной памяти, но сохраняет только половину из имеющихся в начале его работы коэффициентов: те, что находятся на нечетных позициях (независимо от четности N). Скорость работы почти не зависит от того, является ли N-1 целочисленной степенью 2.

Ниже приводятся программы для прогонки и редукции.


/* 
  Thomas algorithm solving tridiagonal linear system. 
	Description:
   int Sweep(double *u,double *a,double *b,double *c,double *f,
             int n,double *bet,double *gam);
	Parameters:
   u - solution (size n), on output;
   a,b,c,f - arrays coefficients (as in the algorithm's description above, size n)
   n - the solution and coefficients arrays size;
   bet, gam - additional arrays (size n+1).
	   Returns: 
   0 - if the algorithm fails,
   1 - if OK.

  Reduction algorithm to solve tridiagonal system.
	Description:
   void Reduce(double *u,double *a,double *b,double *c,double *f,int n);
	Parameters:
   u - solution (size n), on output;
   a,b,c,f - arrays coefficients (as in the algorithm's description above, size n)
   n - the solution and coefficients arrays size;
   Note:
   The algorithm is implenmented without singularity checkup for better 
   performance. To apply it, one must include checkup for division by 
   zero any time the division takes place.

*/

/* Thomas algorithm */
int Sweep(double *u,double *a,double *b,double *c,double *f,int n,
          double *bet,double *gam) {
  register int n;
  double fi;
  /* calculating the additional arrays (factorization of the system) */
  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;
  }
  /* back substitution */
  u[n-1]=gam[n];
  for(i=n-1;i>0;i--) u[i-1]=bet[i]*u[i]+gam[i];
  return(1);
}

/* Reduction techniques */
void Reduce(double *u,double *a,double *b,double *c,double *f,int n) {
  register int i,j;
  int quarter,half,full;
  double a1,c1;
  /* Calculating quarter, central and "right" point numbers */
  i=n; full=1;
  while(i>>=1) full<<=1;
  if(n>full+1) full<<=1;	
  half=full>>1; quarter=half>>1;
  /* Forward pass: coefficients recalculation 
     Note: coefficients at odd positions remain intact! */
  for(i=1;i<half;i<<=1) {
    /* leftmost point */
    c1=c[0]/b[i]; a[0]=0.; c[0]=c[i]*c1; b[0]-=a[i]*c1; f[0]+=f[i]*c1;
    /* other points */
    for(j=i+i;j<n;j+=i+i) {
      /* points placed too far to have influence from the right */
      if(j+i>=n) {
        a1=a[j]/b[j-i]; c[j]=0.; a[j]=a[j-i]*a1; b[j]-=c[j-i]*a1; f[j]+=f[j-i]*a1;
      }
      /* the general case */
      else {
        c1=c[j]/b[j+i]; a1=a[j]/b[j-i];
        c[j]=c[j+i]*c1; a[j]=a[j-i]*a1;
        b[j]-=a[j+i]*c1+c[j-i]*a1; f[j]+=f[j+i]*c1+f[j-i]*a1;
      }
    }
  }
  /* Central coefficients b and f calculation */
  if(full>=n) {		/* case when n!=2^order+1 */
    a1=a[half]/b[0]; b[half]-=c[0]*a1; f[half]+=f[0]*a1;
  }
  else {			/* when it's equal */
    c1=c[half]/b[full]; a1=a[half]/b[0];
    b[half]-=a[full]*c1+c[0]*a1; f[half]+=f[full]*c1+f[0]*a1;
  }
  /* results in the central and leftmost points */
  u[half]=f[half]/b[half]; u[0]=(c[0]*u[half]+f[0])/b[0];
  /* result in the rightmost point if n==2^order+1 */
  if(full<n) u[full]=(a[full]*u[half]+f[full])/b[full];
  /* Backward pass: results in all the other points */
  for(i=quarter;i>=1;i>>=1) for(j=i;j<n;j+=i+i) {
    /* too far to have right influence */
    if(j+i>=n) u[j]=(f[j]+a[j]*u[j-i])/b[j];
    /* the general case */
    else u[j]=(f[j]+c[j]*u[j+i]+a[j]*u[j-i])/b[j];
  }
}

Скачать текст программы tridiag.c

Линейные уравнения
Индекс