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

Written on 13 Мая 2007.

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

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

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

/* 
  Алгоритм Томаса, решающий трехдиагональную линейную систему. 
	Описание:
   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);
}

Редукция

Если провести исключение из имеющейся системы уравнений всех уравнений, находящихся на нечетных позициях, то от системы 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.

/* 
  Алгоритм редукции для решения трехдиагональных систем.
	Описание:
   void Reduce(double *u,double *a,double *b,double *c,double *f,int n);
	Параметры:
   u - разрешение (размера n), на выходе;
   a,b,c,f - массивы коэффициентов (как в алгоритмах, описанных выше, размера n)
   n - размер разрешения и массива коэффициентов;
   Примечание:
   Алгоритм реализован без проверки на производительность. 
*/

/* Метод редукции */
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;
  /* Вычисление четверти, центральной и "правой" точки чисел */
  i=n; full=1;
  while(i>>=1) full<<=1;
  if(n>full+1) full<<=1;	
  half=full>>1; quarter=half>>1;
  /* Первый проход: пересчет коэффициентов 
     Коэффициенты на нечетных позициях остаются нетронутыми! */
  for(i=1;i<half;i<<=1) {
    /* левая точка*/
    c1=c[0]/b[i]; a[0]=0.; c[0]=c[i]*c1; b[0]-=a[i]*c1; f[0]+=f[i]*c1;
    /* другие точки */
    for(j=i+i;j<n;j+=i+i) {
      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;
      }
      /* общий случай */
      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;
      }
    }
  }
  /* Вычисление центральных коэффициентов b и f*/
  if(full>=n) {		/* случай, когда n!=2^order+1 */
    a1=a[half]/b[0]; b[half]-=c[0]*a1; f[half]+=f[0]*a1;
  }
  else {			/* когда они равны */
    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;
  }
  /* результат в центральной и левой точках */
  u[half]=f[half]/b[half]; u[0]=(c[0]*u[half]+f[0])/b[0];
  /* результат в правой точке, если n==2^order+1 */
  if(full<n) u[full]=(a[full]*u[half]+f[full])/b[full];
  /* Обратный проход: результаты во всех остальных точках */
  for(i=quarter;i>=1;i>>=1) for(j=i;j<n;j+=i+i) {
    if(j+i>=n) u[j]=(f[j]+a[j]*u[j-i])/b[j];
    else u[j]=(f[j]+c[j]*u[j+i]+a[j]*u[j-i])/b[j];
  }
}

форум firemonkey

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

ОГЛАВЛЕНИЕ

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

Редукция