Говорухин В., Цибулин Б. Компьютер в математическом исследовании (1185927), страница 73
Текст из файла (страница 73)
Для интегрирования при помощи функции ццаб8 оказалось достаточным 33 точек. С появлением в МАТЕАВ пакета ЯушЬо[!с Ма(Ь Тоо!Ьох на основе символьного ядра Мар!е стало возможным использовать аналитическое вычисление интегралов, см, главу 17 «Расширения МАТЮКАВ». Теперь перейдем к операциям численного дифференцирования.
Вычисление конечных разностей представлено командами, перечисленными в табл. 15.5. Для одномерного массива у вызов функции б1тт(у,1) или б)тт(у) порождает массив с числом элементов, меньшим на единицу. Для двумерного массива разности вычисляются по столбцам. Функция бт гг(у, и) находит конечные разности поряд- 392 Глава 18. Численный анализ в ИАТ[дВ ка п. По заданным значениям аргумента х и функции у производную порядка и можно оценить при помснци отношения 6! ГГ(у, и ) /() ! ГГ(х. п).
Таблица 15.5. Вычисление конечных разностей Назначение Команда Вычисление конечных разностей и-порядка для массива у Вычисление градиента функции двух переменных, заданной массивом Е Вычисление разностного оператора Лапласа для двумерного массива 0 о(гг(у,п) [РХ УУ)-дгаа)епС(Р,НХ.Нт) Ое)2(О.НХ.НУ) В качестве примера использования команды б(((проведем проверку массива точек на монотонность. Напомним, что команда а)) выведет 1 (истина), если условие выполнено для всех элементов, и О (ложь) в противном случае: » т=гао)(4,1)'. 4))(сц/Г(т)>0) т" 0.3784 0.8600 0.8537 0.6936 апз- 0 Ое)2(О,Н1...Нп) Назначения параметров Н1,, Ни аналогичны рассмотренным ранее для функции дгаг)! епС, а выходной массив имеет ту же размерность, что и массив О. Приведем пример вычисления конечных разностей для квадратичной функции ху (седловая функция).
Введем равномерную сетку по координате х и неравномер- ную — по у: » п-10; х--1: 1. /п: 1: у-з!п(р! /2*х); Подготовим массив значений функции: '" ' » [Х.т)-мезидг(о(х.у); 2 Х.*У; Вычислим лапласиан и убедимся, что полученный массив состоит из элементов, близких к значению константы погрешности ерв, это не удивительно при двукрат- ном численном дифференцировании данной функции: Градиент функции п переменных, заданной значениями в виде п-мерного массива Г, вычисляется при помощи команды [Г1О..рп)-дгагдепС(Г,Н1, .Нп) Необязательные параметры Н1, ..., Нп служат для определения шагов сетки по коор- динатам.
В случае равномерной сетки это сами шаги (скалярные величины), а для неравномерной сетки параметры Н1, ..., Нп задают одномерные массивы координат. При обращении дгаг) ! епС(Г) шаги сетки считаются равными едийице, а вызов функ- ции дгаг)(епС(Г.Н) с единственным скалярным параметром Н соответствует вычис- лению градиента с одинаковым шагом по всем коордннатам. Функция г)е) 2 вычисляет одну четвертую разностного оператора Лапласа от за- данной таблицей функции при помощи аппроксимации второго порядка: р)нтерполяцил и приближение функций 393 » 02-0е)2(2,х,у); вах(аах(02)), а)п(ви'п(02)) апв- 2.1301е-014 апв- -2.1301е-014 Построим линии уровня исходной функции и векторное поле градиента функции, см. рис.
15,4. Используем для этого графические команды ьрЬр)от, ссптонг и Сн)тгег (поле направлений): » врЬр)о( П21), соп(оцг(х,у.2), ецЬр1от(122), [П!А2)-чгаг)(епт(2); со)тег(х,у,о1.02) 0.0 0 1 Рис. 15.4. Линии равного уровня (слева) н градиент (справа) седловой функции Отметим, что команда соптонг построила линию уровня с нулевым значением, обой- дя точку х-у-0, Интерполяция и приближение функций Интерполяция позволяет оценить значение заданной таблично функции в тех точках, для которых данных нет. Полиномиальная интерполяция применяется для работы с произвольными данными, а интерполяция, основанная на быстром преобразовании Фурье (БПФ), — для периодических наборов данных. Для обработки экспериментальных данных и в некоторых других случаях нужны команды, позволяющие аппроксимировать функции методом наименьших квадратов. В табл.
15.6 перечислены команды интерполяции. При их описании использованы следующие обозначения: х, У, 7 — заданные массивы (точки и значения функции), г Хт и У1 — массивы точек, в которых должны быть вычислены значения интерполяционного полинома, М вЂ” описание метода интерполирования (линейная, кубическая или сплайн-интерполяция). Для проведения одномерной интерполяции по данным, организованным в таблицу из тОчек Х и значений функций 2, используется; команк(а 21-1 псе гр1(Х, 2.
Х). ' НЕТН00' ) Значения интерполируемой функции будут определены в точках Кт, и выходной массив 2т будет того же размера, что н массив Хтц Дополйительный параметр НЕТН)0 задав г метод интерполирования: по ближайшей точке (пеагеят), линейная интер- 394 Глава 35. Численный анализ в МАТ(Д8 поляция (Бпеаг), которая действует по умолчанию, кубическая (сцЫс) или интерполяция кубическим сплайном (зрйпе). Для всех методов величины элементов Х должны меняться монотонно.
Для узлов Х1, лежащих вне интервала, определенного Х, соответствующим значениям выходного массива будут присвоены нечисловые константы Иан. Если параметр Х отсутствует, то считается, что значения функции вычислены на равномерной целочисленной сетке Х=1:И, где Иев12е(7,1). Таблица 35.6. Коианды интерполяции Команда Описание Одномерная интерполяция Двумерная интерполяция Трехмерная интерполяция (3 массива координат) п-мерная интерполяция (и массивов координат) Интерполяция периодической функции Интерполяция на неравномерной сетке Интерполяция кубическим сплайном Определение методом наименьших квадратов полинома заданной степени И,аппроксимируюк(его таблицу чисел !птегр1(Х.7.Х1,'М') тптегр2(Х.Т,Х.Х!,т1,'М') !т\Сегрз(Х ., 7,Х! ..., М ) !птегрп(Х,....7.Х! .....'М') (пСЕгртт(7,И) 9г!паата(Х.Т,7,Х),т!) зр)зпе(Х.7.Х1) ро)ут)С(Х,7,И) Приведем простой пример.
Интерполируем функцию синуса на промежутке 10,21 и вычислим значения для точек О, 1, 2, 3. Для последней точки числового значения нет, поскольку точка х=З вне интервала, на котором задана таблица данных: » х-0:.5:2; х!-1птегр1(х.з1п(х).(0:31.'соо!с') з! 0 0.8415 0.9093 Иан Многомерный случай отличается от одномерной интерполяции большим числом массивов для задания координат. Рассмотрим команду и-мерной интерполяции: !птегрп(Х1,...,Хп.7,Хз1, ,Хтп,'МЕТМ00') В качестве параметра МЕТНОО используются те же имена, что и для одномерной интерполяции: йпеаг, сцо(с, пеагезт и зрйпе. Каждая и-мерная матрица из списка Х1, ..., Хп определяет значения одной координаты функции и переменных, представленной массивом значений 7.
Координаты Х1, ..., Хп должны изменяться моиотонно, причем в случае равномерного расположения узлов для ускорения интерполирования можно указывать в качестве параметра МЕТНОО значения ')1пеаг, *сц))1с, *пеагезт. Параметр со звездочкой означает, что вычисления ведутся по формулам с постоянным шагом, что и приводит к ускорению. Для двумерных матриц и массивов большей размерности имеется возможность быстро получить расширецный массив данных; вычислив дрполнительнырзнвчения для новых узлов, располагающихся между старыми узлами. При вызове 1птегр2(7,1), 1п1егрЗ(7,1) или 1птегрп(7,1) массив значений будет расширен по каждой размерности за счет вставки промежуточных интерполированных точек. При вызове 1птегр2(7,М), 1птегрЗ(7.И) или 1псегрп(7,И) подобное расширение рекуррентно будет проведено И раз.
Интерполяция и приближение функций 395 Продемонстрируем действие данной команды для двумерного массива ЗхЗ: »2 2 1 2 0 В результате выполнения команды мером 5х5: » 1п(егр2(2, 1. ' спЫ с ' ) дП5 1.0000 1.8750 2.0000 1.3750 0 при М-1 получаем расширенную матрицу раз- 5.0000 7.3750 8.0000 6.8750 4.0000 2.0000 3.0000 4,0000 3.2500 4.6250 6.0000 3.5000 5.0000 6.5000 2,7500 4.1250 5.5000 1.0000 2.0000 3.0000 При многомерной интерполяции массивы точек создаются при помощи команды ппдг! (( (аналог двумерной команды пезпды ((). Время, необходимое для интерполя- ции, зависит от длины массива.
Задавая в качестве дополнительного параметра обрабатываемое число точек, можно ускорить процесс обработки данных. Напри- мер, можно указать число, являющееся степенью двойки или просто выбирая чис- ло, меньшее числа точек. Если указывается число, превосходящее реально имею- щееся, то массивы будут пополнены нулями. Функция дг)((пдгд приМеняЕтСя для ДВУмерной интерполяции данных, опреДелен- ных в точках, координаты которых даются векторами Х и У, и получения массива значений 21 для упорядоченных массивов координат Х1, У1. Отметим, что в этом случае значения не обязательно упорядочены: 21-дг!Ода!а(Х.7.2.Х!.У() Приведем пример. Зададим поверхность сорока точками, вычисленными случай- ным образом: » х-2«гало(40,1)-1: у 2*гало(40,1)-1, г«чип(х «у); Построим равномерную сетку и интерполируем в ее узлах данные: » С -1:.1;1; [Х1Л()-аедндг!О(1); 2-дг!Опаса(х,у,Х,Х(,у!); Построим полученную интерполированием поверхность и нанесем исходные точки: » еедн(Х)Л1,2), По)О("оп"), Р)оСЗ(х.г.х.'ох') Обратим внимание, что поверхность на рис.
15.5 построена недда всех точек, опре- деленных массивами х1, 71. это связано с тем, что соответствующие элементы мас- Не приводя сам массив, укажем размерность результата для той же исходной матрицы прн вызове » д!хе(!п(егр2(2.3,'сиЫс')) дП5 17 17 396 Глава 15. Численный анализ в МАТ(АО сива имеют нечисловой характер (МОМ). Действительно, выведем несколько злементов второй строки: » 7(2,1:5) длз МаМ Ман -0.6503 -О.б025 Ман -1 -1 Рис. 15.5. Набор случайных тачек и результат действия коианды Онддата Команда зр11пе(Х, У, Х1 ) выведет значения кубического сплайна для узлов Хт по таблице (Х.У) аналогично тому, что делает команда 1псегр1(Х,У,Х1.'зр)тпе'). Кроме того, вариант атой команды Рреар11пе(Х, У1) позволяет получить сам сплайн (так называемая рргоглт форма) и в дальнейшем вычислять его значения при помощи функции ррчз).