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

Подпрограмма:  ALS1R (модуль ALS1R_p)

Назначение

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

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

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

(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.

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

procedure ALS1R(var A :Array of Real; var U :Array of Real;
                var F :Array of Real; M :Integer; N :Integer;
                EPS :Real; var H :Array of Real; var W :Array of Real;
                L :Integer);


   Параметры
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 обращается к вспомогательным подпрограммам с именами: ALS1R1, ALS1R2, ALS1R3, ALS1R4.

  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. 

Unit TALS1R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, ALS1R_p;

function TALS1R: String;

implementation

function TALS1R: String;
var
I :Integer;
U :Array [0..2] of Real;
H :Array [0..49] of Real;
W :Array [0..14] of Real;
const
A :Array [0..8] of Real = ( 1.0,2.0,3.0,1.0,2.0,3.0,1.0,1.0,2.0 );
F :Array [0..2] of Real = ( 2.0,4.0,6.0 );
M :Integer = 3;
N :Integer = 3;
EPS :Real = 0.0001;
begin
Result := '';  { результат функции }
ALS1R(A,U,F,M,N,EPS,H,W,1);
Result := Result + Format('%s',['               ТЕСТ ПОДПРОГРАММЫ ALS1R']) + #$0D#$0A; 
Result := Result + Format('%s',[' СИСТЕМА УРАВНЕНИЙ:' + #$0D#$0A]) + #$0D#$0A; 
Result := Result + Format('%s',['    X +   Y  +   Z =  2' + #$0D#$0A +
 '  2*X + 2*Y  +   Z =  4' + #$0D#$0A +
 '  3*X + 3*Y  + 2*Z =  6' + #$0D#$0A]) + #$0D#$0A; 
Result := Result + Format('%s',['     EPS: ']);
Result := Result + Format('%11.3f ',[EPS]) + #$0D#$0A;
Result := Result + #$0D#$0A;
for I:=1 to N do
 begin
  Result := Result + Format(' РЕШЕНИЕ: %11.3f ',[U[I-1]]) + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
Result := Result + #$0D#$0A;
for I:=1 to 7 do
 begin
  Result := Result + Format('       H: %11.3f ',[H[I-1]]) + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
UtRes('TALS1R',Result);  { вывод результатов в файл TALS1R.res }
exit;
end;

end.

Результаты: 

      EPS:           1.000E-04 
      Решение:   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