Получение профилей и контурных картин поверхностей вращения способом триангуляции с использованием зеркально отражeнного излучения, алгоритмы и расчeтные формулы Текст научной статьи по специальности «Математика»

Математическое моделирование

УДК 531.717.81 (088.8) (520)

О. А. Заякин

ПОЛУЧЕНИЕ ПРОФИЛЕЙ И КОНТУРНЫХ КАРТИН ПОВЕРХНОСТЕЙ ВРАЩЕНИЯ СПОСОБОМ ТРИАНГУЛЯЦИИ С ИСПОЛЬЗОВАНИЕМ ЗЕРКАЛЬНО ОТРАЖЁННОГО ИЗЛУЧЕНИЯ, АЛГОРИТМЫ И РАСЧЁТНЫЕ ФОРМУЛЫ

Описан и обоснован вывод расчётных формул и приведены алгоритмы решения задачи получения профилей и контурных картин поверхностей вращения способом триангуляции с использованием зеркально отражённого излучения. Отмечена перспективность данного способа для задачи активного контроля отклонений формы обработанных поверхностей.

В технике часто применяются детали с зеркально отражающими поверхностями вращения, имеющими высокое качество обработки. Такие поверхности требуют высокого технологического уровня контроля, а также бережного к себе отношения. Поэтому методы контроля к ним стремятся применять чувствительные и неразрушающие. В анализе микрогеометрии названного класса поверхностей преобладают оптические методы, поскольку они позволяют достичь большего пространственного разрешения и поддаются тарированию. Их применение способствует решению таких актуальных задач, как повышение производительности, экономичности и достоверности операций контроля, а также внедрение массового и активного контроля.

Из оптических методов контроля указанного класса поверхностей предпочтительней использовать те методы, в которых в качестве источника информации используется зеркальная, а не диффузная, компонента отражённого света. Это связано с тем, что контролируемая поверхность имеет следы обработки с сильно выраженной направленностью. Они оказывают сильное влияние на индикатрису рассеяния диффузной компоненты и её зависимость от ориентации контролируемой поверхности относительно зондирующего луча. В процессе сканирования контролируемой поверхности сложной формы это может приводить к нестабильности мощности излучения, попадающего в приёмную оптическую систему, а, следовательно, и уровня входного сигнала. Особенно это сказывается при картировании неровностей субмикронного диапазона высот, характерного для зеркально отражающих поверхностей, поскольку в таких задачах оптические профилометры приближаются к пределу своих возможностей и, соответственно, требования к качеству входного сигнала ужесточаются.

В настоящее время в названной предметной области используют устройства на основе интерференционных методов [1]. В основном используют поляризационные интерферометры [2, 3]. Однако до настоящего времени они не смогли вытеснить традиционно применяемые в данной области контактные профилометры. Прогресс идёт в направлении развития комбинированных устройств, в которых в качестве первичного преобразователя используется контактный щуп, а для дальнейшего преобразования сигнала используют оптические методы [4-6]. Очевидно, это является компромиссным решением.

Одна из основных причин, сдерживающих распространение оптических методов, применяемых сейчас в названной предметной области, состоит в том, что контролируемая поверхность вращения сложной формы вносит недопустимые искажения во входной, оптический сигнал. Исправление этих искажений в устройстве, способном работать с поверхностями вращения разнообразных форм, требует применения дорогостоящих элементов адаптивной оптики.

Метод триангуляции с использованием зеркально отражённого излучения может быть приемлемой альтернативой указанному подходу. Он позволяет обойтись без адаптивной оптики. Реализующая его оптическая схема более простая и дешёвая, чем интерферометр.

В работах [7, 8] впервые была обоснована применимость данного метода для решения задачи контроля формы рабочих поверхностей внутренних колец шарикоподшипников.

Проблема тарирования, сдерживавшая применение названной разновидности метода триангуляции, была успешно решена в работах, выполненных при участии автора данной статьи [9-11], и в работах других учёных [12, 13], выполненных независимо друг от друга.

В данной работе подробно описан и обоснован вывод расчётных формул и приведены алгоритмы решения задачи получения профилей и контурных картин поверхностей вращения способом триангуляции с использованием зеркально отражённого излучения.

Ограничимся рассмотрением таких случаев, когда контролируемая поверхность имеет слабый рельеф, то есть её отклонения от круглости незначительны по сравнению со средним радиусом профиля.

Рассматриваемые в данной работе контролируемые изделия имеют вид однополостного тела вращения с криволинейной наружной поверхностью вида

где р, р, г — цилиндрические координаты, Я(г) — функция номинального профиля, Н(р, 2) — функция отклонения от номинальной формы.

Рассмотрим модель оптического сканирования в приближении геометрической оптики при условии, что падающий и отражённый пучки света являются тонкими [14, с. 168]. При этом важно, что, во-первых, поперечная ширина тонкого пучка не учитывается в расчётах и, во-вторых, прохождение тонкого пучка в пространстве описывается положением его центрального луча.

На рис. 1 изображена тестируемая поверхность в зоне оптического контроля. Параллельно оси OZ со смещением на величину й в положительном направлении оси ОУ проведена прямая АВ. В направлении данной прямой производится сканирование поверхности вдоль оси вращения источником излучения. Луч а, падающий на поверхность, перпендикулярен плоскости УOZ и зеркально отражается относительно вектора нормали N.

Отражённый луч Ь регистрируется матричным фотоприёмником, расположенным так, чтобы столбцы матрицы были параллельны оси С этой осью совмещается серединный столбец матрицы. Для поиска луча Ь, отражённого от поверхностей с различными наклонами, перед началом измерений фотоприёмник имеет возможность перемещаться по дуге радиуса Ь в диапазоне углов [ф1, ц2].

Информативными параметрами являются угол ^, который отсчитывается от положительного направления оси ОХ и который соответствует энергетическому центру отражённого пучка на матричном фотоприёмнике, и координата $ этого энергетического центра на фотоприёмнике.

Сканирование всей поверхности осуществляется дискретным перемещением источника света по прямой АВ в точках с координатами = %к-1 + Дг и шаговым разворотом контролируемого

объекта вокруг оси OZ. Алгоритм сбора данных приведён на рис. 2.

Получим формулы для определения координат освещённой точки контролируемой поверхности на каждом шаге сканирования.

В уравнении падающего луча а заданы две из трёх координат — это у и г. Третью, неизвестную координату, можно найти из условий зеркального отражения луча света. Как известно, они описываются векторными соотношениями:

Р(р, р, z) = р - Я (г) - Н(р, z) = 0,

г

Рис. 1. Оптико-механическая схема координатных измерений: 1 — источник излучения; 2 — контролируемая поверхность; 3 — матричный фотоприёмник; АВ, С — направления

сканирования поверхности

-а х N N хЬ (-а ■ N (Ь ■ N

\а\ \Ь\ ’ 1а1 \Ь\ '

Эти соотношения определяют принадлежность входящих в них векторов одной плоскости и равенство углов падения и отражения. Записывая данные соотношения в развёрнутом виде и дополняя их уравнением падающего луча, получаем математическую модель оптической схемы в виде системы из четырех уравнений:

Ру ■ Ъ2 Р2 ■ Ьу — О,

Рх-Ъя-Ря-Ьх

х 2 2 х _ р, —о,

у!ЪХ+ъ2у+ъ2

Рх ■ Ъу _ Ру ■ Ъх

^Ъ2х + ъу + Ъу

_ Ру — О,

(1)

Рх ■ Ъх + Ру ■ Ъу + Р2 ■ Ъ2

УЪу+Ъу+Ъу

_ Рх — О,

где Ьх, Ьу, Ьг — проекции луча Ь на соответствующие координатные оси; Рх, Ру, Рг — частные производные функции поверхности в освещённой точке.

Для упрощения уравнений (1) переходим от декартовых координат {х, у, z}, к системе цилиндрических координат {рс, рс, гс}. Располагаем её так, чтобы направления осей г и гс совпадали, а отсчёт углов рс начинался от положительного направления оси ОХ в сторону положительного направления оси ОУ. Уравнения связи между обеими системами имеют вид:

х — Рс ■ 0ОБРс, у — Рс ■ БІПРс,

В новом базисе теперь две координаты из трёх оказываются неизвестными — это рс и рс. Однако в оптической схеме они связаны между собой следующим соотношением

й — Рс ■ БІП Рс

(2)

Решая уравнение (1) в новом базисе, какую-нибудь одну из них можно выразить через другую. В результате опять остаётся только одна неизвестная переменная, в качестве которой естественно выбрать радиус рс.

Следующие уравнения представляют собой связь между производными функции поверхности, записанными в двух выбранных нами системах координат оптической схемы:

Рх — 0ОБ Рс ---------

дрс

дР БІП рс

Рс

Ру — Біп Рс + ~-------

дрс

дР ООБ Рс

Рс

дР

д2с '

В (3) предполагается, что функцию Р можно представить в виде:

Р = р -/{р> г) =0,

где / — функция отклонения от формы поверхности вращения. Тогда дР/др = 1.

(3)

(4)

Подставив (3) в (1), можно видеть, что математическая модель приобретает вид системы нелинейных дифференциальных уравнений, содержащей кроме одной неизвестной координаты рс ещё две неизвестные величины — производные от функции поверхности дР/дрс и дР/дгс.

Попытаемся найти неизвестную величину рс, рассматривая уравнения (1) как систему четырёх алгебраических уравнений с тремя неизвестными.

Используя тот факт, что уравнения (1) в новом базисе являются линейными относительно неизвестных дР/дрс и дР/дгс, выразим из этой системы указанные производные в явном виде. Для этого, оказывается, достаточно взять любые два уравнения из четырёх, имеющихся в системе (1), так как, если рассматривать эту систему как линейную относительно дР/дрс и дР/дхс, то матрица её коэффициентов имеет ранг, равный двум. Это означает, что два оставшихся уравнения, из которых путём подстановки исключены неизвестные дР/дрс и дР/дгс, представляют собой тождества.

Таким образом, решить эту систему относительно рс не удаётся. Но можно получить из неё дифференциальное уравнение, связывающее рс и его первую производную по какой-либо другой координате рс или гс. Покажем это.

Выбрав, для определённости, первое и четвёртое уравнения из системы (1), получаем выражения для производных следующего вида:

дР дР

-— = -рс ■ Ф {рс, гс, ря, гя), — = -Г (рс, гс, ря, гк), (5)

дрс дгс

где ря — угол поворота контролируемого объекта, гя — смещение источника излучения. Перечисленные переменные задаются независимо, известным способом. Их отсчёт ведём в той же выбранной нами системе цилиндрических координат оптической схемы.

Функции в правых частях уравнений (5) имеют следующий вид:

Ф [Рс, 2с, Ря, 2я) — г (Рс, 2с, Ря, 2я) —

[сОБ (у _ Рс) _ А _ В ■ ООБ рс] • (б1пу _ А ■ БІП Рс) + tgу в ■ БІП Рс [б1п (у _ Рс) + В ■ БІПрс] ■ (б1пу _ А ■ БІПРс) +tgу в ■ ООБРс (А ■ ООБ Рс _ ООБу _ В) ■ tgв

(6)

[б1п (у _ Рс) + В ■ БІП рс] ■ (б1п у _ А ■ БІП Рс) +tgу в ■ ООБ Рс’

где

А — —, БІпРс — —, tgв — -——, В ^/Ї^Ау_2^А^СО$(у_РсУ+їдУ0.

Т п~ Г V ^ }

Отметим, что В ■ Ь — \Ъ|, а величины у, - и в — есть функции переменных Ря и 2я, задаваемых дискретными отсчётами в процессе сканирования.

Используя, что в рассматриваемой оптической схеме 2с — 2я, исключаем переменную 2с из правой части рассматриваемых дифференциальных уравнений.

Из (5) получаем дифференциальные уравнения относительно неизвестного радиуса контролируемой поверхности в её освещённой точке:

дР = Рс ■ Ф (Рс, Ря, 2я) , ^ — г (Рс, Ря, 2я) , (7)

в которых правые части уравнений представляем как производные неявной функции радиуса Рс, зависящей от переменных Рс и 2с и определяемой уравнением поверхности (4).

Запишем уравнения (7) в системе координат {р, р, 2}, жёстко связанной с контролируемой поверхностью. Уравнения связи между системами координат имеют следующий вид:

Р — Рс, Р — Рс _РА, 2 — 2с, (8)

где ра — угловая координата начала отсчёта (р — 0). Она изменяется с поворотом поверхности при сканировании в соответствии с формулой:

рА — ря _ рСОшЬ (9)

где рСОш1; — постоянная величина угла. Её мы вводим для удобства расчётов. Она определяет начальное положение контролируемого объекта в оптической схеме. Для угла ря принимаем,

что он изменяется в процессе сканирования от 0 до 2п, или до —2п радиан, в зависимости от направления сканирования. Все формулы в данной работе записаны, предполагая сканирование в положительном направлении (см. правую часть рис. 1).

Из (8) следует, что ^ ^, др = дрс. Тогда уравнения (7) переходят в следующие уравнения:

др

— = р • Ф (р, р«, Zr) , (10)

др = Г (р, р«, Zr) . (11)

Функция вида р(р«, zr) является частным решением любого из них в отдельности. Из неё,

используя (2) и (8), можно получить профиль радиуса р(р, z) в системе координат контролиру-

емой поверхности, а из сетки таких профилей — и изображение всей поверхности.

Так, решая уравнение (10) при каком-либо постоянном zr, получаем функцию р(ря, zr), от которой, используя (2), (8) и (9), переходим к р(р, z) при заданном z = zr. Граничное условие

р(ря = 0, Zr) = р(ря = 2п, Zr)

при малых отклонениях р от своей средней величины вдоль восстанавливаемого профиля поверхности позволяет получить однозначное решение при замене нелинейного дифференциального уравнения (10) упрощённым дифференциальным уравнением, полученным путём его линеаризации в окрестности средних интегральных, по этому профилю, величин у и в.

Для уравнения (11) граничного условия нет, но его можно получить, решив уравнение (10) при названных выше условиях.

Из (10) получаем дифференциальное уравнение следующего вида:

др = Фapprox (р, V,в) = Ф («О, Уо,во) + dФ ^ • (р — R0) +

+ diRVoA}• , — уо) + dФ(я£о,во) • (в —во), (12)

где «о, Vo и во — средние интегральные величины р, у и в по радиальному профилю.

Подставим в первое из уравнений (6) вместо у и в значения ,о и во и приравняем нулю его правую часть. Тогда р в этом уравнении будет равно искомому «о.

Уравнение (12) переходит в исходное уравнение (10) в предельном случае, а именно, когда максимальная высота рельефа с ограниченным по частоте спектром пространственных гармоник стремится к нулю. Нетрудно показать, что и решение приближенного (линеаризованного) уравнения (12) стремится в этом пределе к точному восстановлению функции профиля контролируемого объекта.

Мы рассматриваем линеаризацию уравнения (10), а не (11), поскольку этот путь выглядит более предпочтительным. В самом деле, при картировании поверхности сеткой радиальных профилей условие малости относительных отклонений от круглости в пределах отдельного профиля выполняется для любой поверхности вращения. В то же время при сканировании по прямой AB это верно только для цилиндра.

При сканировании контролируемой поверхности имеют место вариации переменной рс вокруг своего среднего значения, которые прямо связаны с радиусной координатой поверхности в освещённой точке (см. (2)). Но при условиях, когда допустима замена исходного уравнения (10) на приближенное линеаризованное (12), этими вариациями можно пренебречь. Тогда

р = —ря + const,

и можно считать, что производная в левой части уравнения (12) берётся по независимой координате, то есть, координату р можно считать в них независимой переменной.

При этом допущении уравнение (12) решается как обыкновенное дифференциальное уравнение. Получим его аналитическое решение. Запишем (12) в следующем виде:

1 (dH ^

--------a • H = b • е + с • y, (13)

«о Ы р I

где a = Яо • dФ(яо,во), b = dФ(Уво), с = dФ(яо,во), е = у — Уо, Y = в — во, у = у(р, z), в = в(р, z).

Граничное условие для (13) имеет вид:

H (р = о, z) = H (р = 2п, z) . (14)

Для каждого значения z имеется своё уравнение (13), отличающееся значением коэффици-

ентов a, b, с.

Применим для его решения хорошо известный метод Бернулли [15, с. 29-31]. Решение имеет вид:

H (р) = Яо • | С + J [b • е (р1) + с • y (р1)] • exp [—a • р — р)] dр-i |,

где

C =

п2л

I b■ ep) + С■ jp) ■ exp [-a■ p -p)] dpi

Jo

exp (-2n ■ a) - i

С целью аппроксимации, либо фильтрации измерительной информации, функцию

W (p) = b ■ e(p) + с ■ j[p) (15)

(см. обозначения в (13) удобно представить в виде разложения в ряд Фурье. В этом случае решение дифференциального уравнения (13) можно получить в общей форме через соответствующие коэффициенты разложения. Для данных в цифровой форме при этом наиболее подходит дискретное преобразование Фурье [16, с. 483-484]. Напишем, как выглядит решение уравнения (13), полученное с помощью этого преобразования.

Пусть функция W(p) задана своими дискретными отсчётами W(pi) с равномерным шагом по радиальному профилю; i = 0,1,2,..., N - i; p0 = 0; pN = 2п. Тогда она однозначно представима в виде конечного ряда Фурье:

[N/2]

W [pt) = £ a j ■ cos(j ■ p¡) + в j ■ sinj ■ pt), (16)

j=i

где pt = 2nN, [■ ] — целая часть. Коэффициенты этого ряда Фурье вычисляются по формулам:

2 N -1 2 N -1 1 N -1

aj = Vr ■£ W (pO ■ cos(j ■ pO’ Pj = Tt ^ W ■ sin(j ■ p^’ aN/2 = ту ■ Y, (-1)i ■W M’ (17)

N i=0 N i=0 N i=0

где i < j < N/2. Тогда искомый радиальный профиль представлен дискретными отсчётами вы-

соты H(pi) в N точках: при тех же значениях pi, что и W (pi):

f [N/2] 1 1

H [p¡) = R0 ■< Y. -2--2 ■ (n ■ an - a ■ Pn) ■ sin (n ■ p¡) - (a ■ an + n ■ Pn) ■ cos (n ■ p¡) >, (18)

[ n=1 a + n J

где i =0,1, 2 ,N - 1.

Таким образом, найдена алгебраическая связь между параметрами гармоники Фурье-спек-тра функций e и у и соответствующей гармоники (то есть, с тем же числом волн на оборот поверхности) Фурье-спектра профиля высот контролируемой поверхности.

Выражения (15)- (18) представляют собой функцию преобразования измерительной схемы. При получении радиального профиля используются данные, собранные только с него самого. В заключение сделаем замечание относительно случая, когда контролируемый объект имеет внутреннюю рабочую поверхность. Описанная конфигурация схемы измерений пригодна в этом случае, если поверхность имеет достаточный диаметр, чтобы внутри неё поместился фотоприёмник, с учётом необходимого расстояния от него до поверхности. На практике в ряде случаев это условие, к сожалению, не выполняется. При этом требуется изменить конфигурацию схемы таким образом, чтобы, луч, отражённый от контролируемой поверхности, не был бы перекрыт контролируемой деталью. Это можно сделать, например, отклоняя падающий луч в плоскости XOY. Или поступить иначе — наклонить ось вращения объекта в неизменной схеме, то есть сделать эту ось пересекающейся с осью OZ. Тогда задача сводится к рассмотренной. Нужно только после получения профиля сделать преобразование координат к удобному

виду. Но в обоих вариантах это потребует изменений алгоритмов и программ картирования поверхности, хотя их основа — система векторных уравнений — останется прежней. Отметим, что здесь могло бы помочь решение рассмотренной задачи в общем случае — при произвольной ориентации падающего и отражённого лучей.

На основании изложенного можно сделать следующие выводы.

Разработанные алгоритмы, основанные на регистрации координат сканирующего пучка света, зеркально отражённого от контролируемой поверхности, могут быть пригодны для картирования поверхностей вращения разнообразных форм.

Полученная аналитическая функция преобразования измерительной схемы, используемая в указанных алгоритмах, является корректной, что доказано конструктивным путём, следуя математическим правилам, в рамках выбранного физического приближения, в качестве которого использовано геометрооптическое приближение тонких пучков света.

Разработанные алгоритмы обладают следующими характерными свойствами:

- на основе одних и тех же данных они позволяют совместить измерения как среднего радиуса, так и высот профиля поверхности;

- позволяют автоматически настраиваться на контролируемые объекты широкого диапазона размеров;

- отвечают одному из основных направлений развития современных стандартов, связанному с усилением роли статистической обработки результатов отдельных измерений.

Автор благодарит В.Н. Белопухова за большую помощь в работе.

БИБЛИОГРАФИЧЕСКИЙ СПИСОК

1. Осипович, И. Р. Интерферометрический метод контроля формы асферических поверхностей качения прецизионных подшипников [Текст] / И. Р. Осипович, Д. Т. Пуряев // Вестник МГТУ. Сер. Приборостроение. — 1999. — №. 3. — С. 65 -75.

2. Bristow, T. C. Surface topography goes beyond just rms roughness [Text] / T.C. Bristow, R. A. Auriemma // Laser Focus World. — 1992. — Vol. 28, No. 3. — P. 148-150.

3. Sato, Atsushi Высокоточный профилометр типа Maxim 3D-5700 [Текст] / Atsushi Sato // Кэйсоку гидзюцу. (Яп.) = Instruments and Automation. — 1991. — Vol. 19, No. 2. — P. 54-58.

4. Morrison, E. The development of a prototype high-speed stylous profilometer and its application to rapid 3D surface measurement [Text] / E. Morrison // Nanotechnology. — 1996. —Vol. 7, No. 1. —P. 37-42.

5. Lu, S. A composite probe for surface topography measurement [Text] / S. Lu, Z. Li // International Journal of Advanced Manufacture Technology. — 20O1. — Vol. 18, No. 11. —P. 831-835.

6. Каршилов, Е. Метрололия в машиностроении. Измерительное оборудование на выставке «Метрмаш’ 2007» [Электронный ресурс]/ Е. Каршилов // Снабженец. — 2007. — 25 (573).—С. 44-47. —Режим доступа: http://www. snab.ru/stati/25_3.pdf

7. Technological system for control of bearing ring runway surface [Text] / V.N. Belopukhov, V.G. Volostnikov, V. N. Podvigin, O. A. Zayakin // Proc. SPIE. — 1996. — Vol. 27, No. 13. — P. 461-463.

8. Задача контроля профиля тел вращения при угловом сканировании и одномерной регистрации [Текст] / В. Н. Белопухов, О. А. Заякин, В. А. Катулин, В.Н. Подвигин // Труды ФИАН. — 1993. —Т. 217.—C. 170-176.

9. Способ измерения геометрической формы тел вращения с отражающей поверхностью [Текст]: пат. 2109250 Рос. Федерация: МПК6 G 01 B 11/24 / В. Н. Белопухов, С. И. Бесталанный, О. А. Заякин; заявитель и патентообладатель Самарский филиал фИаН. —№ 95100536/28 (001062); заявл. 12.01.95; опубл. 20.04.98, Бюл. № 11.

10. Заякин О. А. Информационно-измерительная система для контроля деталей подшипников, экспериментальная оценка точности восстановления микрорельефа рабочих поверхностей [Электронный ресурс] / О. А. Заякин // Электронный журнал «Исследовано в России». — 2004. — 187. — С. 1992-2001. — Режим доступа: http : //zhurnal. ape.relarn.ru/articles/2004/187.pdf.

11. Заякин О. А. Информационно-измерительная система контроля деталей подшипников на основе двумерной лазерной триангуляции [Текст]: дисс. ... к-та техн. наук: 05.11.16: защищена 24.10.2005: утв. 10.02.2006 / Олег Александрович Заякин. — Самара, 2005. — 178 с.

12. Устройство для контроля параметров криволинейных поверхностей [Текст]: пат. 2025659 Рос. Федерация: МПК6 G 01 B 11/24 / Ю. Я. Андрейченко, В. А. Волошинов, Д. В. Волошинов, В. В. Самсонов; заявитель и патентообладатель Ленинградский государственный технический университет. — № 4944595/28; заявл. 29.05.1994; опубл. в 1994, Бюл. № 24.

13. Baba, Mitsuru Determining 3-D shapes of hybrid surfaces by laser rangefinder [Text] / Mitsuru Baba, Kozo Ohtani // Proc. SPIE. — 2001. — Vol. 4398. — P. 86-97.

14. Борн, М. Основы оптики [Текст] / М. Борн, Э. Вольф; изд. 2-е. — М.: Наука, 1973. — 720 с.

15. Бугров, Я. С. Высшая математика. Дифференциальные уравнения. Кратные интегралы. Ряды. Функции комплексного переменного [Текст] / Я. С. Бугров, С. М. Никольский. — М.: Наука, 1981. —448 с.

16. Каханер, Д. Численные методы и математическое обеспечение [Текст] / Д. Каханер, К. Моулер, С. Нэш.— М.: Мир, 1998. — 575 с.

Физический институт им. П. Н. Лебедева РАН, Самарский филиал Поступила 07.11.2006

olegz@fian.smr.ru