OPEM2                    Библиотека  "JINRLIB"                   
                                       
    Автор: Богданова Н.Б.
      (Institute of Nuclear Research and Nuclear Energy, BAS,
       Tzarigrasko chaussee 72, 1784 Sofia, Bulgaria)                Вы
    Язык: Фортран 77                                                 
    Платформа: Windows                                               посетитель.
   
       МЕТОД РАЗЛОЖЕНИЯ ПО ОРТОНОРМИРОВАННЫМ ПОЛИНОМАМ ОТ ДВУХ ПЕРЕМЕННЫХ
                 И ЕГО РЕАЛИЗАЦИЯ В ВИДЕ ПАКЕТА ПРОГРАММ OPEM2
  
    Пакет программ OPEM2 разрабатывался с 1981 для калибровочных задач
    измерительных систем в физике высоких энергий в сотрудничестве с 
    ЛИТ ОИЯИ, г.Дубна. 

    Последние версии пакета сделаны при поддержке Болгарского национального
    фонда научных исследований по физике - Phi 1001.

    1.Описание метода
   
    В OPEM2 используется метод ортогонализации полиномов, основанный на 
    трехчленном соотношении Хаусхолдера-Форсайта [1,2] (для одномерного 
    случая) и рекуррентном соотношении Вейсфельда [3] (для многомерного 
    случая). Семейство полиномов {PL} задается рекуррентным соотношением:
                       (1.1)
    где Pk - базовый полином. K и I вычисляются в процессе работы программы.
    {PL} определяются на конечном дискретном подмножестве D пространства ,
    D={q1,q2,...,qM}. Каждая точка qj=qj(x1j1,x2j2,...,xnjn) представляет собой
    n-координатный вектор (в данном случае n=2). Веса wj=1/S2j, ассоциированные
    с этими точками, зависят от стандартных отклонений в каждой точке.
    Точки, веса и значения экспериментальной функции {qj,wj,fjexp}j=1M
    в этих точках должны быть заданы. Рекурентные коэффициенты {L}, {L-1}
    и нормализующий фактор {cL} вычисляются как скалярное произведение
    от заданных значений. Данное описание есть развитие работ [3,4,7].
    В OPEM2 аппроксимируемая функция разлагается по ортонормированным
    полиномам с использованием специальных ортонормированных коэффициентов ak.
    Оптимальное значение степени полиномов Lr вычисляется, исходя из двух
    критериев.
                                       (1.2)
    Тогда ортонормированные коэффициенты {ak} в (1.2) легко вычисляются
    из заданных значений:
                                         (1.3)
    Представим двумерные полиномы в виде пирамиды в порядке возрастания
    степени полиномов: 
              
    В каждой ячейке пирамиды сверху записан номер полинома, а степени j1 и j2
    первого и второго аргумента x1 и x2 приведены ниже.

    Пример вычисления индексов K и I в (1.1) для полинома P13 (k=2, j1=2, j2=2):
                       (1.4)
    (P8 - базовый полином: K=8, I=4).

    В работе [7] для двумерного случая представлены новые теоретические
    исследования - для случая дифференцирования и интегрирования в уравнение
    (1.1) добавлен четвертый член. Более детальное описание математического
    метода для многомерных полиномов, включая дифференцирование и
    интегрирование, будет опубликовано позднее в подходящем журнале
    в соответствии с настоящим описанием и с некоторыми дополнениями.

    2.Описание пакета OPEM2 (Фортран 77)

    OPEM2 состоит из главной программы CALIB и пяти подпрограмм: 
    TSTWO, ORTTWO, PRETWO (основные), DEGREE, NUMBER (вспомогательные).

    В главной программе CALIB задаются:
       M     - количество точек (узлов);
       X и Y - массивы координат точек;
       F     - значения аппроксимируемой функции в узлах;
       W     - значения весовых коэффициентов.

    Основные подпрограммы содержат COMMON-блок /LINKS/, куда записываются
    коэффициенты рекуррентного соотношения и вспомогательные переменные.
    Размеры массивов в /LINKS/ зависят от суммы индексов j1 и j2, и когда эта
    сумма равна 16, размеры массивов определяются как Lmax=(16+1)(16+2)/2=153.
    Они могут изменяться в зависимости от задачи.

    Подпрограмма TSTWO(X,Y,W,M,F,A,NDEG) организует вызов других подпрограмм. 
    Сначала TSTWO вызывает PRETWO для подготовки коэффициентов рекуррентных 
    соотношений с предельной степенью NDEG, подсчитывает ортонормированные
    коэффициенты A, вызывая подпрограмму ORTTWO, определяет аппроксимируемые
    значения и отклонения в точках, выбирает оптимальную степень Lr и 
    минимум ()2, подсчитывает ошибки и печатает результаты.

       X и Y - одномерные массивы координат точек (узлов);
       W     - значения весовых коэффициентов в каждом узле;
       M     - количество точек;
       F     - заданная функция;
       A     - ортонормированные коэффициенты;
       NDEG  - максимальная предельная степень полиномов.
  
    Подпрограмма ORTTWO(NUMBMX,X,Y,POLY) в POLY-массиве возвращает конкретные
    значения полиномов от 1 до NUMBMX при заданных значениях x и y. 

       NUMBMX - значение Lr - оптимальная степень полиномов;
       X и Y  - координаты точек, где вычисляются полиномы;
       POLY   - одномерный массив результатов.
     
    Подпрограмма PRETWO(M,MAXD,X,Y,W,POLY) подготавливает коэффициенты 
    рекуррентного соотношения двумерных полиномов, ортонормированных на 
    дискретном точечном множестве.

       M     - количество точек;
       MAXD  - максимальная предельная степень, для которой вычисляются 
               коэффициенты рекуррентного соотношения (как NDEG в TSTWO);
       X и Y - одномерные массивы, содержащие координаты точек;
       W     - значения весовых коэффициентов;
       POLY  - одномерный массив результатов.

    PRETWO делает последовательный рекурсивный вызов подпрограммы ORTTWO,
    начиная с NUMBMX = L-1. Следующий вызов ORTTWO выполняется со значением
    L, увеличенным на 1, и так далее до тех пор, пока L не достигнет 
    значения Lmax.
   
    Подпрограмма DEGREE(N,J,JX,JY) считает для данного порядкового номера N 
    предельную степень его главного члена J, степень x-переменных JX и 
    степень y-переменных JY.

    Функция NUMBER(JX,JY) наоборот, по заданным JX и JY возвращает порядковый 
    номер полинома в соответствии с [3].

    3.Применение
    
    В работах [5,6,7] описаны применения пакета OPEM2 для получения 
    калибровочных коэффициентов между идеальной и измеренной сетками в 
    оптических измерительных приборах физики высоких энергий. Последнее
    применение пакета - для эксперимента H1 в DESY (Hamburg) - нахождение
    координат максимума поверхности, полученной суммированием гауссианов.
    Аналитическая поверхность построена на сетке, заданной более чем 400
    точками. Работа сделана в коллаборации с ЛИТ ОИЯИ (г.Дубна).

    4.Литература

    1. A.Householder, Principles of numerical analysis
       (Mc Graw-Hill, New York,1953), p.221.
    2. G.Forsythe, J.Soc.Ind.Math. 5 (1957) 74.
    3. W.Weisfeld, Numer.Math. 1 (1959) 38.
    4. V.Gadjokov, Commun.JINR E11-80-782, Dubna (1980).
    5. N.Bogdanova, V.Gadjokov, Comp.Phys.Commun.,24 (1981) p.225-229.
    6. N.Bogdanova, J.Comp.and Appl.Mathematics, 14 (1986) p.345-351.
    7. N.Bogdanova, T.Kupenova, Comp.Phys.Commun.,73 (1992) p.170-178.
 
    Архив программы OPEM2  с исходными текстами и результатами.