Текст подпрограммы и версий
als1r_c.zip 		Тексты тестовых примеров
tals1r_c.zip

Подпрограмма:  als1r_c

Назначение

Вычисление устойчивого приближения к нормальному псевдорешению переопределенной системы линейных алгебраических уравнений методом регуляризации А.Н.Тихонова с выбором параметра регуляризации по невязке.

Математическое описание

Решение системы линейных aлгeбpaических уравнений: 

(1)                            A * U  =  F

методом регуляризации А.Н.Тихонова [1] сводится к минимизации параметрического функционала: 

(2)                           || A * U - F || 2  +  ALFA * || U || 2 ,

где ALFA > 0 - параметр регуляризации. Значение параметра регуляризации ALFA определяется из условия: 

(3)                           RO(ALFA)  =  E * || F || ,

где:
RO(ALFA) - невязка на регуляризованном решении U[ALFA](X),
E                 - заданный относительный уровень невязки,
|| F ||            - евклидова норма правой части уравнения (1).

Задача минимизации (2) по U эквивалентна решению системы линейных алгебраических уравнений: 

(4)      (TRANSP(A) * A + ALFA * I ) * U[ALFA]  =  TRANSP(A)*F ,

при этом параметр ALFA выбирается из условия: 

(5)                   RO(ALFA)  =  || A * U[ALFA] - F ||  =  E*|| F || ,

где U[ALFA] - решение (4) ,  I - единичная матрица.

При решении системы (4) используется сведение задачи к более простой "канонической" задаче с двухдиагональной матрицей A согласно методу В.В.Воеводина [2].

Для решения уравнения (5) используется метод касательных Ньютона в соответствии с вычислительной схемой, предложенной В.А.Морозовым [3].

Подпрограмма вычисляет семь характеристик полученного регуляризированного решения, которые могут быть использованы для оценки его качества, а именно, в подпрограмме вычисляются приближенные:

1.  невязка  RO (ALFA) = || A*U [ALFA] - F || 2.  функция "чувствительности"  TAU (ALFA) = ALFA*|| DU [ALFA] / D [ALFA] || 3.  норма решения  GAMMA(ALFA) = || U [ALFA] || 4.  значение регуляризирующего функционала
FI(ALFA) = [RO (ALFA)] 2 + ALFA* [GAMMA (ALFA)] 2 5.  найденное значение ALFA как решение уравнения (5) 6.  достигнутый относительный уровень невязки 7.  количество итераций при решении уравнения (5)

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

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

1)  решение интегрального уравнения для фиксированных правой части и относительного уровня невязки (это требует порядка MN 2 операций). 2)  решение того же уравнения, но с другим относительным уровнем невязки (это требует порядка N 2 операций). 3)  решение интегрального уравнения с теми же ядром и стабилизирующим функционалом, но с другой правой частью и относительным уровнем невязки (это требует порядка MN операций).
 
1.  Тихонов А.Н., О регуляризации некорректно поставленных задач, ДАН СССР, 1963, т.153, N 1. 2.  Воеводин В.В. О методе регуляризации, ЖВМ и МФ, 1969, т.9, N3. 3.  Морозов В.А. О принципе невязки при решении несовместных уравнений методом регуляризации А.Н.Тихонова, ЖВМ и МФ, 1973, т.13, N 5.

Использование 

    int als1r_c (real *a, real *u, real *f, integer *m, integer *n,
             real *eps, real *h, real *w, integer *l)

   Параметры
a - вещественный двумерный массив размера m на n, содержащий заданную матрицу системы линейных уравнений;
u - вещественный вектор длины n, в результате работы подпрограммы содержит значения найденного решения;
f - вещественный вектор длины m, содержащий вектор правых частей системы линейных уравнений;
m - число строк матрицы A (тип: целый);
n - число столбцов матрицы A (тип: целый);
eps - заданный относительный уровень невязки (тип: вещественный);
h - вещественный вектор длины 50. В результате работы подпрограммы содержит вычисленные характеристики решения:
h(1) - невязка,
h(2) - функция "чувствительности",
h(3) - норма регуляризованного решения,
h(4) - значение регуляризирующего функционала,
h(5) - значение параметра регуляризации,
h(6) - достигнутый относительный уровень невязки,
h(7) - число итераций;
w - двумерный вещественный рабочий массив размера n на 5;
l - параметр, определяющий режим использования подпрограммы (тип: целый):
l= 1 - решение системы первый раз (режим i);
l= 2 - решение системы с новым значением относительного уровня невязки eps (режим ii);
l= 3 - решение системы с новыми правой частью f и относительным уровнем невязки eps (режим iii). 
   Версии:  нет 
   Вызываемые подпрограммы:  нет 
   Замечания по использованию
  1. 

Подпрограмма als1r_c обращается к вспомогательным подпрограммам с именами: als1r1_c, als1r2_c, als1r3_c, als1r4_c.

  2. 

При первом обращении к подпрограмме (l = 1) необходимо задать параметры, определяющие систему: a, f, m, n, eps, l.
В дальнейшем можно менять только значения eps (при l = 2) и f, eps (при l = 3).
Значения остальных перечисленных выше параметров (и рабочего массива w) должны сохраняться, т.к. они используются в программе при всех значениях l.

  3. 

Если требуемое значение eps не может быть достигнуто (при нахождении параметра регуляризации ALFA методом Ньютона значение ALFA достигает наименьшего положительного числа, представимого на ЭВМ, или при уменьшении ALFA невязка начинает возрастать), то вычисляется наилучшее возможное решение (с минимальной на вычисленных значениях ALFA невязкой или при наименьшем возможном значении ALFA > 0).

  4.  Значение относительного уровня невязки eps должно быть неотрицательно: eps ≥ 0. Если eps ≥ 1, то начальное приближение u = 0 является решением задачи.

Пример использования

Рассмотрим решение системы линейных алгебраических уравнений
(точное решение - прямая  x + y = 2,  z = 0): 

                 x  +      y  +      z  =  2 
             2*x  +  2*y  +      z  =  4 
             3*x  +  3*y  +  2*z  =  6

Ниже приводится фрагмент программы, вычисляющей решения и соответствующие значения критериальных функций при заданном значении относительного уровня невязки eps. 

int main(void)
{
    /* Initialized data */
    static float a[9] /* was [3][3] */ = { 1.f,2.f,3.f,1.f,2.f,3.f,1.f,1.f,
                                           2.f };
    static float f[3] = { 2.f,4.f,6.f };
    static int m = 3;
    static int n = 3;
    static float eps = 1e-4f;

    /* System generated locals */
    int i__1;

    /* Local variables */
    extern int als1r_c(float *, float *, float *, int *, int *, float *,
                       float *, float *, int *);
    static float h__[50];
    static int i__;
    static float u[3], w[15] /* was [3][5] */;

    als1r_c(a, u, f, &m, &n, &eps, h__, w, &c__1);

    printf("\n eps = %11.3e \n", eps);
    printf("\n u = %11.3e %11.3e %11.3e \n", u[0], u[1], u[2]);
    for (i__ = 1; i__ <= 7; ++i__) {
         printf("\n h = %11.3e \n", h__[i__-1]);
    }
    return 0;
} /* main */


Результаты: 

      eps:           1.000e-04 
      peшeниe:   9.993e-01   9.997e-01   1.587e-03
      h__:             7.483e-04   1.770e-03   1.413e + 00   1.061e-03 
                         5.307e-04   1.000e-04   2.000e + 00