|
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 с исходными текстами и результатами.
|