|
БИБЛИОТЕКА ПРОГРАММ JINRLIBOPEM2 - метод разложения по ортонормированным полиномам от двух переменных и его реализация в виде пакета программАвтор: Богданова Н.Б. Institute of Nuclear Research and Nuclear Energy, Sofia, Bulgaria |
|
Язык: Фортран 77, Windows Пакет программ OPEM2 разрабатывался с 1981 для калибровочных задач измерительных систем в физике высоких энергий в сотрудничестве с ЛИТ ОИЯИ, г.Дубна. Последние версии пакета сделаны при поддержке Болгарского национального фонда научных исследований по физике - Phi 1001.
1.Описание метода В OPEM2 используется метод ортогонализации полиномов, основанный на трехчленном соотношении Хаусхолдера-Форсайта [1,2] (для одномерного случая) и рекуррентном соотношении Вейсфельда [3] (для многомерного случая). Семейство полиномов {PL} задается рекуррентным соотношением:
где 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 вычисляется, исходя из двух критериев.
Тогда ортонормированные коэффициенты {ak} в (1.2) легко вычисляются из заданных значений:
Представим двумерные полиномы в виде пирамиды в порядке возрастания степени полиномов: В каждой ячейке пирамиды сверху записан номер полинома, а степени j1 и j2 первого и второго аргумента x1 и x2 приведены ниже. Пример вычисления индексов K и I в (1.1) для полинома P13 (k=2, j1=2, j2=2):
(P8 - базовый полином: K=8, I=4). В работе [7] для двумерного случая представлены новые теоретические исследования - для случая дифференцирования и интегрирования в уравнение (1.1) добавлен четвертый член. Более детальное описание математического метода для многомерных полиномов, включая дифференцирование и интегрирование, будет опубликовано позднее в подходящем журнале в соответствии с настоящим описанием и с некоторыми дополнениями.
2.Описание пакета OPEM2 (Фортран 77) OPEM2 состоит из главной программы CALIB и пяти подпрограмм: TSTWO, ORTTWO, PRETWO (основные), DEGREE, NUMBER (вспомогательные). В главной программе CALIB задаются:
Основные подпрограммы содержат 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, подсчитывает ошибки и печатает результаты.
Подпрограмма ORTTWO(NUMBMX,X,Y,POLY) в POLY-массиве возвращает конкретные значения полиномов от 1 до NUMBMX при заданных значениях x и y.
Подпрограмма PRETWO(M,MAXD,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.Литература:
Архив программы OPEM2 с исходными текстами и результатами. |