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

Подпрограмма:  EV01R (модуль EV01R_p)

Назначение

Решение интегрального уравнения Вольтерра со стационарным ядром методом рягуляризации первого порядка без выбора параметра регуляризации.

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

Данная подпрограмма реализует алгоритм, изложенный в [1].

Численное решение интегрального уравнения Вольтерра 1 рода вида 

                T
(1)           ∫   k (t - x)  u (x)  dx  =  f (t)
              t
                t
           (   ∫   k (t - x)  u (x)  dx  =  f (t) )
              0

         t ∈ [ 0, T ]

на равномерной сетке  ti = i H,  i = 1, 2, ..., n,  H = Т/n, сводится к задаче минимизации 

(2)           || K u - f || E - min u ,

где K - верхняя (нижняя) треугольная теплицева матрица размера  n на n, полученная при дискретизации (1):   Ki j  =  H * k (H ( i - j )).

Используя метод регуляризации А.Н.Тихонова [2], приходим к следующей задаче 

(3)           || K u - f || E2  +  α || D u || E2  -  min u ,

где α > 0 - заданный параметр регуляризации, D - дискретный аналог оператора дифференцирования   d /dx.

Задача (3) эквивалентна задаче 

                ||     |    K      |           |  f   |     ||  
(4)           ||     |             |  u   -   |      |     ||      -   min u ,                 
                ||     | α1/2 D |           |  0  |     || E

для решения которой при помощи метода вращений [3] определим ортогональную матрицу Q, такую, что 

              |     K    |            |  Kα  |                             |   f    |            |   f1    |
    QT     | α1/2 D |     =     |   0    |  ,               QT     |   0   |     =     |   f2    |

где Kα - верхняя (нижняя) треугольная матрица. Тогда решение задачи (4) сводится к решению системы линейных алгебраических уравнений с треугольной матрицей: 

(5)           Kα u  =  f1 ,

котоpое осуществляется методом обратной подстановки.

1.  L.Elden. An efficient algorithm for the regularization of ill - conditioned least squares problems with triangular toeplitz matrix. - REPORT  LITH - MAT - R - 1981 - 25, Sweden, 1981.
2.  А.Н.Тихонов, В.Я.Арсенин. Методы решения некорректных задач. - М., Наука, 1979.
3.  В.В.Воеводин. Численные методы алгебры. - M., Hаука, 1966.

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

procedure EV01R(var A :Array of Real; var F :Array of Real; N :Integer;
                L :Integer; ALFA :Real; H :Real; var B :Array of Real;
                var C :Array of Real; var U :Array of Real;
                var IERR :Integer);


   Параметры
A - вещественный вектоp длины N, содержащий первую стpоку (первый столбец) теплицевой матрицы K;
F - вещественный вектоp длины 2N, содержащий в первых N ячейках компоненты заданного вектоpа правой части; последующие N ячеек используются как рабочие;
N - заданная размерность вектоpов F и U и число узлов сетки разбиения (тип: целый);
L - заданный параметр - указатель типа матрицы K:
 

L = 0, если матрица верхняя треугольная,

L = 1, если матрица нижняя треугольная (тип: целый);

ALFA - заданный параметр регуляризации (тип: вещественный);
H - заданный шаг разбиения сетки (тип: вещественный);
B - вещественный вектоp длины N (N + 1)/2, используемый как рабочий;
C - вещественный вектоp длины N, используемый как рабочий;
U - вещественный вектоp длины N, содержащий на выходе из подпрограммы вычисленное решение интегрального уравнения в узлах заданной сетки;
IERR - целая переменная, сообщающая об ошибках, обнаруженных в ходе вычислений;
IERR= 0 - если ошибок не было обнаружено;
IERR=65 - в случае, если при заданном значении ALFA матрица системы вырождается. 
   Версии:  нет 
   Вызываемые подпрограммы
UTEV10 - подпрограмма выдачи диагностических сообщений при работе подпрограммы EV01R. 
   Замечания по использованию:  нет

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

 Решается интегральное уравнение
                      t
                     ∫   K (t - s)  f (s)  ds  =  g (t)   ,      0 ≤ t ≤ 1
                   0
с ядром
                              ∞
      K(t)  =  0.5 π   ∑    (- 1)n (2nH )  exp (-(2n + 1)2 π2 t /8) .
                             n=0

Точное значение функции  f (s) задано таблицей. Используя формулу трапеции, вычисляем приближенное значение правой части  g (t), затем по  g (t) восстанавливаем значение функции  f (s) = u с помощью программы EV01R.

Параметр регуляризации  α = 10 - 8 

Unit TEV01R_p;
interface
uses
SysUtils, Math, { Delphi }
LStruct, Lfunc, UtRes_p, EV01R_p;

function TEV01R: String;

implementation

function TEV01R: String;
var
N,N1,I,J,L,IERR :Integer;
ALFA,H,PI,S,RАВ :Real;
R :Array [0..49] of Real;
G :Array [0..99] of Real;
U :Array [0..49] of Real;
WORK :Array [0..1324] of Real;
K :Array [0..49] of Real;
const
F :Array [0..49] of Real = ( 0.025,0.215,0.405,0.595,0.785,0.975,1.0,0.92,
0.71,0.5,0.29,0.08,0.0,0.08,0.29,0.5,0.71,0.92,
1.0,0.975,0.785,0.595,0.405,0.215,0.025,0.0,0.0,
0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0 );
label
_12,_11,_13,_15,_14;
begin
Result := '';
N := 50;
ALFA := 1.E-8;
H := 2.E-2;
PI := 3.14159265;
N1 := N-1;
for I:=1 to N do
 begin
  G[N+I-1] := 0;
  S := Exp(-(PI*PI*I*H/8));
  for J:=1 to 15 do
   begin
    RАВ := Exp(-IntPower(2*J+1,2)*PI*PI*I*H/8);
_12:
    S := S+Power(-1,J)*(2*J+1)*RAB;
   end;
_11:
  R[I-1] := 0.5*PI*S;
 end;
K[0] := R[0]*H/2;
K[N-1] := R[N-1]*H/2;
for I:=2 to N1 do
 begin
_13:
  K[I-1] := R[I-1]*H;
 end;
for I:=1 to N do
 begin
  S := 0;
  for J:=1 to I do
   begin
_15:
    S := S+F[J-1]*K[I-J+0];
   end;
_14:
  G[I-1] := S;
 end;
L := 1;
EV01R(R,G,N,L,ALFA,H,WORK[N+0],WORK,U,IERR);
Result := Result + Format('%s',[' ALFA=']);
Result := Result + Format('%20.16f',[ALFA]) + #$0D#$0A;
Result := Result + Format('%s',[' U=']);
Result := Result + #$0D#$0A;
I := 1;
WHILE ( I<=N ) do
 begin
  Result := Result + Format('%20.16f ',[U[I-1]]) + #$0D#$0A;
  inc(I,5);
 end;
Result := Result + #$0D#$0A;
UtRes('TEV01R',Result);  { вывод результатов в файл TEV01R.res }
exit;
end;

end.

Результаты:

          ALFA  =  1.00000000-08

               |  4.18939426-02   9.47262731-01   2.69209987-01    4.94652786-01
      U  =  |  8.07735236-01   9.34217208-04   8.94488047-03   -2.22705150-03
               | -3.81489539-04   5.87342988-03