Решение линейных систем с трехдиагональной матрицей
Страница 3. Редукция

Редукция

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

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

Если число уравнений не является степенью двойки, алгоритм реализуется методом фиктивного восполнения области определения до степени 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];
 }
}