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

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

Назначение

Решение одномерного интегрального уравнения I рода с разностным ядром на прямой (- ∞, + ∞) методом регуляризации с применением алгоритма Быстрого Преобразования Фурье и (или) вычисление значений критериальных функций.

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

Приближенное решение одномepнoгo интегрального уравнения I рода типа свертки 

                   +∞
(1)    Ax  ≡   ∫  k(t - τ)  x(τ) dτ  =  y(t)  ,    -∞ < t < +∞ ,
                  -∞

с гладким ядром  k осуществляется методом однопараметричecкой регуляризации А.Н.Тихонова [1] путем сведения задачи к минимизации по  x сглаживающего функционала 

(2)     Mα [x, y]  =  || Ax - y ||2 + α Ω[x] ,
 где :
                                        +∞
       || Ax - y ||2  =  1/2π    ∫   | K(λ) X(λ) - Y(λ) |2 dλ  -  квадрат невязки,
                                       -∞
      α > 0 - параметр регуляризации,
                             +∞
       Ω[x]  =  1/2π  ∫  | λ |2p | X(λ) |2 dλ  -  стабилизирующий функционал
                            -∞
 порядка  p ≥ 0 (p не предполагается целым);

    K (λ), X (λ), Y (λ) - преобразования Фурье ядра  k, решения  x и приближенной правой части  y. 

 Функционал (2) минимизирует функция
                                +∞
(3)     xα(t)  =  1/2π  ∫  K*(λ) Y(λ) / ( | K(λ) |2 + α λ2p ) exp (iλt) dt ,
                              -∞
          -∞ < t < +∞ ,
 где  *  - знак комплексного сопряжения.

Считается, что  k (t) → 0,  x (t) → 0,  y (t) → 0 при  t → ±∞ и функции  k (t),  y (t) заданы на конечном отрезке [-Т/2, Т/2] в узлах равномерной сетки по  t : 

     ts  =  s δt ,   s =  -N/2, -N/2 + 1, ... , -1, 0, 1, ... , N/2 - 1 ,

где  N - четное число узлов,  δ (t) = T/N - шаг сетки.

Задача состоит в решении дискретного уравнения типа свертки с ядром  ks = k (ts) и приближенной правой частью  ys = = y (ts) методом регуляризации  p - ого порядка в спектральной области  λ. Численный алгоритм основан на применении Дискретного Преобразования Фурье временных рядов [2].

Приведем задачу к "каноническому" виду, вычислив ДПФ функций  ks и  ys : 

                         N/2-1
             Km   =    ∑      ks exp(-i λm ts) ,
                        S = -N/2
                         N/2-1
             Ym   =    ∑      ys exp(-i λm ts) ,
                       S = -N/2

где  λm = m δλ ,  m =  - N/2, - N/2 + 1 , ..., - 1, 0, 1, ..., N/2 - 1 ,  δλ = 2π/(N δt) ,  i = √-1.

Аналогично определяются  Xm,  Xmα  -  ДПФ рядов  xs = x (ts)  и  xsα =xα (ts).

Далее, аппроксимируя интегралы в pавенствах (2), (3) и в преобразованиях Фурье по формуле прямоугольников, получаем дискретные аналоги сглаживающего функционала  Mα [xs, ys] и минимизирующего его регуляризованного решения 

                                 N/2-1
(4)   xsα  = 1/ Nδt      ∑      KTmYm /( | Km |2 + α/δt2 λm2p ) exp(2π i m s / N) ,
                               m= -N/2
        s  =  -N/2, -N/2 + 1, ..., N/2 - 1 .

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

1) невязка    ρ(α)  =  || Axα - y ||  ,         

2) ноpма решения     γ(α)  =  Ω1/2 [ xα ]  ,    

3) значение    регуляризующего    функционала 
    φ(α)  =  [ ρ2(α) + αγ2(α) ]1/2  ,
4) функция "чувствительности"
   (для  выбора  квазиоптимального значения  α)
    τ(α)  =  α Ω1/2 [ dxα/dα ] ,
   где
                           +∞
    Ω[x]  =  1/2π  ∫   | λ | | X(λ) |2 dλ  ≡  Ω[x]
                         -∞
    при  p = 1/2.

Вычисление всех критериальных функций происходит в спектральной области  λ (через компоненты ДПФ). Это требует порядка  N операций и значительно облегчает задачу выбора параметра регуляризации. Восстановление искомого решения по формуле (4) целесообразно производить лишь для выбранных значений  λ.

Подпрограмма предусматривает работу в следующих режимах (в зависимости от значения параметра  L).

1.  Задача приводится к "каноническому" виду; этот режим позволяет избежать повторного вычисления  Km,  Ym, когда происходит повторное обращение к подпрограмме с тем же ядром и той же правой частью, например, при выборе  α  (L = 1).
2.  При условии, что задача приведена к "каноническому" виду, строится решение  xsα и вычисляются значения крикритериальных функций  ρ (α),  γ (α),  φ (α),  τ (α) при заданном значении  α (L = 2).
3.  При условии, что задача приведена к "каноническому" виду, вычисляются только значения критериальных функций  ρ (α),  γ (α),  φ (α),  τ (α) при при заданном значении  α без построения решения  xsα (L = 3).
4.  Вычисляется только ДПФ правой части  ys; этот режим необходим для приведения задачи к "каноническому" виду, когда повторно решается уравнение с тем же ядром, но с другой правой частью (L = 4).

Для оптимального вычисления прямого и обратного ДПФ применяся алгоритм быстрого преобразования Фурье [2], за счет которого реализованный в подпрограмме алгоритм численного решения задачи является эффективным в смысле быстродействия, экономии памяти ЭВМ и уменьшения ошибок округления, при этом время вычислений пропорционально  N log2N, а объем памяти пропорционален  N.

Полное описание изложенного алгоритма содержится в статье [3].

1.  А.Н.Тихонов. O регуляризации некоppектно поставленных задач, ДАН CCCP, 1963, т.153, N 1, 49 - 52.
2.  В.А.Морозов, Н.Н.Кирсанова, А.Ф.Сысоев. Комплекс алгоритмов быстрого преобразования Фурье дискретных рядов, сб. "Численный анализ на ФОРТРАНе", вып. 15. Изд - во МГУ, M., 1976.
3.  М.В.Арефьева. Быстрый регуляризирующий алгоритм решения интегральных уравнений Фредгольма I рода типа свертки, сб. "Численный анализ на ФОРТРАНе. Методы и алгоритмы". Изд - во МГУ, M., 1978.

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

    int ec00r_c(integer *n, real *dt, real *a, real *y, real *alp,
         real *p, integer *l, real *h, real *x, real *z)

   Параметры
n - количество заданных значений ядра и правой части уравнения,  n ≥ 4 - целая степень числа два (тип: целый);
dt - заданный шаг сетки по переменным  t,  τ,  dt > 0 (тип: вещественный);
a - вещественный вектоp длины  n, содержащий заданные значения ядра уравнения на сетке:  a (i) = kI-n/2-1,  i = 1, 2, ..., n;
y - вещественный вектоp длины  n, содержащий заданные значения правой части уравнения на сетке:  y (i) = yI- n/2- 1,  i = 1, 2, ..., n;
alp - заданное значение параметра регуляризации  α,  alp ≥ 0 (тип: вещественный);
p - заданное значение порядка регуляризации  p,  p ≥ 0 (тип: вещественный);
l - параметр, определяющий режим использования подпрограммы (тип: целый):
l = 1 - выполняется приведение задачи к "каноническому" виду,
l = 2 - строится решение и вычисляются значения критериальных функций,
l = 3 - вычисляются только значения критериальных функций;
l = 4 - выполняется приведение задачи с новой правой частью с тем же ядром к "каноническому" виду;
h - вещественный вектоp длины 4, содержащий вычисленные значения критериальных функций при заданном значении  alp:  h (1) = ρ(α),  h (2) = γ(α),  h (3) = φ(α),  h (4) = τ(α);
x - вещественный вектоp длины  n, содержащий вычисленные значения регуляризованного решения на сетке:  x (i) = xαI- n/2- 1,  i = 1, 2, ..., n;
z - вещественный вектоp длины  n, используемый в подпрограмме как рабочий. 
   Версии:  нет 
   Вызываемые подпрограммы
ftf1c_c - подпрограмма вычисления дискретного или обратного дискретного преобразования Фурье комплексного ряда длины, равной степени двух, методом быстрого преобразования Фурье. 
   Замечания по использованию
  1. 

B результате работы подпрограммы при  l = 1 в массиве  a содержится преобразование Фурье ядра  km, а в массиве  y - преобразование Фурье правой части  ym; на выходе подпрограммы при  l = 3 в массиве  x содержится преобразование Фурье регуляризованного решения  xmα.

C учетом симметрии вещественной части и антисимметрии мнимой части преобразований Фурье вещественных рядов относительно центра  m = n/2, хранятся только половины значений  km,  ym,  xmα,  m = 0, 1, ..., n/2.

Например, в массиве  a спектр располагается следующим образом: 
     a(i) = re kI-1 ,         i = 1, 2, ..., n/2 + 1 ,
     a(i) = im kI-n/2-1 ,   i = n/2 + 2, n/2 + 3, ..., n ;

при этом учитывая также, что  im km = 0,  m = 0, n/2.

  2. 

При первом обращении к подпрограмме (l = 1) необходимо задать параметры  n, dt, a, y, определяющие уpавнение. B дальнейшем (l = 2, 3) можно менять только параметры  alp и  p, характеризующие регуляризатор, сохраняя значения  n, dt, a, y. При  l = 4 значения параметров  n, dt, а также должны сохраняться.

  3.  Массивы  a, y, x используются при всех значениях  l, а рабочий массив  z - только в одном режиме при  l = 2.

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

Рассматривается решение интегрального уравнения типа свертки ( с ядром  k (t) = exp (- t2) и правой частью  y (t) = (π/2) 1/2 exp (- t2/2) ( точное решение  x (t) = exp (- t2) ) или y1 (t) = (π/3) 1/2 exp (- 2/3 t2) ( точное решение  x1 (t) = exp (- 2t2) ).

Ниже приводится фрагмент программы, вычисляющей решения и соответствующие значения критериальных функций при фиксированном параметре регуляризации  α = 10- 2 с использованием регуляризатора порядка  p = 1 и равномерной на отрезке [ - 1, 1 ] сетки из  n = 8 точек. 

int main(void)
{
    /* Initialized data */

    static int n = 8;
    static float dt = .25f;
    static float alp = .01f;
    static float p = 1.f;
    static float pi = 3.14159265359f;

    /* Builtin functions */
    double exp(double), sqrt(double);

    /* Local variables */
    extern int ec00r_c(int *, float *, float *, float *,
                       float *, float *, int *, float *, float *, float *);
    static float a[8], h__[4];
    static int i__, i, l;
    static float t, x[8], y[8], z__[8], y1[8], xt[8], xt1[8];
    int i__1;

    i__1 = n;
    for (i__ = 1; i__ <= i__1; ++i__) {
        t = dt * (-n / 2 + i__ - 1);
        a[i__ - 1] = (float)exp((float)(-t * t));
        xt[i__ - 1] = (float)exp((float)(-t * t));
        y[i__ - 1] = (float)sqrt((float)(pi / 2.f)) *
                (float)exp((float)(-t * t / 2.f));
        xt1[i__ - 1] = (float)exp((float)(t * -2.f * t));
/* l1: */
        y1[i__ - 1] = (float)sqrt((float)(pi / 3.f)) *
                 (float)exp((float)(-t * t * 2.f / 3.f));
    }
    for (i = 0; i <= 4; i+= 4) {
        printf("\n %16.7e %16.7e %16.7e %16.7e \n",
                  a[i], a[i+1], a[i+2], a[i+3]);
    }
    for (i = 0; i <= 4; i+= 4) {
        printf("\n %16.7e %16.7e %16.7e %16.7e \n",
                  xt[i], xt[i+1], xt[i+2], xt[i+3]);
    }
    for (i = 0; i <= 4; i+= 4) {
        printf("\n %16.7e %16.7e %16.7e %16.7e \n",
                  y[i], y[i+1], y[i+2], y[i+3]);
    }
    l = 1;
    ec00r_c(&n, &dt, a, y, &alp, &p, &l, h__, x, z__);
    l = 2;
    ec00r_c(&n, &dt, a, y, &alp, &p, &l, h__, x, z__);

        printf("\n %16.7e %16.7e %16.7e %16.7e \n",
                  h__[0], h__[1], h__[2], h__[3]);
    for (i = 0; i <= 4; i+= 4) {
        printf("\n %16.7e %16.7e %16.7e %16.7e \n",
                  x[i], x[i+1], x[i+2], x[i+3]);
    }
    for (i = 0; i <= 4; i+= 4) {
        printf("\n %16.7e %16.7e %16.7e %16.7e \n",
                  xt1[i], xt1[i+1], xt1[i+2], xt1[i+3]);
    }
    for (i = 0; i <= 4; i+= 4) {
        printf("\n %16.7e %16.7e %16.7e %16.7e \n",
                  y1[i], y1[i+1], y1[i+2], y1[i+3]);
    }
    l = 4;
    ec00r_c(&n, &dt, a, y1, &alp, &p, &l, h__, x, z__);
    l = 2;
    ec00r_c(&n, &dt, a, y1, &alp, &p, &l, h__, x, z__);

        printf("\n %16.7e %16.7e %16.7e %16.7e \n",
                  h__[0], h__[1], h__[2], h__[3]);
    for (i = 0; i <= 4; i+= 4) {
        printf("\n %16.7e %16.7e %16.7e %16.7e \n",
                  x[i], x[i+1], x[i+2], x[i+3]);
    }
    return 0;
} /* main */


Результаты:

       h__  = ( 0.13239,  1.08822,  0.17138,  0.33298 )
       x    = ( 0.378,  0.475,  0.713,  0.963,  1.072,  0.963,  0.713,  0.475 )



       l = 4
       ec00r_c(&n, &dt, a, y1, &alp, &p, &l, h__, x, z__);
       l = 2;
       ec00r_c(&n, &dt, a, y1, &alp, &p, &l, h__, x, z__);

Результаты:

       h__  = ( 0.13270,  1.11089,  0.17306,  0.33989 )
       x    = ( 0.211,  0.310,  0.554,  0.809,  0.919,  0.809,  0.554,  0.310 )