Трехдиагональная линейная система уравнений может быть записана в следующем виде:
Матрица этой системы состоит из главной диагонали (с коэффициентами -bi) и примыкающих к ней сверху и снизу диагоналей (соответственно ci и ai).
Прямое (безытерационное) решение этой системы уравнений можно проводить методом ее факторизации (прогонка, или алгоритм Томаса), либо методом последовательного исключения сначала нечетных, затем делящихся на 2, но не на 4, и т.д. уравнений и соответствующего преобразования коэффициентов. Последний метод называется редукцией или бинарным исключением.
Алгоритм Томаса
Решение ищется в виде:
Алгоритм принципиально не может иметь ведущего элемента, и есть вероятность того, что он не сойдется даже для несингулярной матрицы. Для этого в программе, реализующей прогонку, необходимо контролировать значения fi. Еще два замечания: длины рекуррентных цепочек порядка N, поэтому при экспоненциальном накоплении погрешностей (что происходит при преобладании значений |beti|>1) прогонка практически обязательно разойдется. Второе: алгоритм сохраняет исходные значения коэффициентов. Доказано, что для устойчивой работы алгоритма Томаса достаточно диагонального преобладания, однако, во многих случаях он сходится и при отсутствии такового.
Редукция
Общее описание алгоритма таково.
Если провести исключение из имеющейся системы уравнений всех уравнений, находящихся на нечетных позициях, то от системы N исходных уравнений мы перейдем к системе из (N-1)/2+1 уравнений только для четных значений u, нечетные же значения будут вычисляться через четные. Повторяя эту процедуру для укороченной системы, мы в конце концов придем к единственному уравнению для узла, лежащего в позиции (N-1)/2 (только если N-1 - степень двойки). Затем, проводя обратную подстановку, определяются крайние, а затем и все прочие узлы.
Формулы для пересчета коэффициентов на каждом этапе исключения таковы:
Если число уравнений не является степенью двойки, алгоритм реализуется методом фиктивного восполнения области определения до степени 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];
}
}