4.3.
Функция трех и более аргументов
В предыдущей главе были даны рекомендации по «оцифровке» графиков – по созданию их «after ego»: таблиц, хранящих координаты точек на графиках. Кроме того, как уже отмечалось ранее, справочники сами содержат помимо «безтабличных» графиков и, если так можно выразиться, «безграфиковые» таблицы. В данной главе будет рассказано о том, как таблицы, хранящихся в справочниках или «выуженные» из графиков, превращаются в формулы по которым:
а) можно считать (оценивать) промежуточные (интерполяция) или «запредельные» (экстраполяция) данные;
б) допустимо построить сам график, если его не было, или воссоздать имеющийся в более четком исполнении (сравните рис. 3.4 и 3.10); кроме того, график можно воспроизвести в другой системе координат (см. рис. 3.10 и 3.11);
в) воссоздать функциональные зависимости в других программных средах и других программных продуктах, перенося формулы в них из соответствующих сайтов.
Кстати, о программных продуктах. Очень часто их разработчики не публикуют заложенные в них функциональные зависимости, считая это лишним или намеренно скрывая их, боясь появления на рынке конкурирующих продуктов.
Так вот, эти зависимости можно при необходимости воссоздать приемами, описанными в данной главе, предварительно выяснив, сколько у функции аргументов и не являются ли они в свою очередь сами функциями других аргументов и протабулировав ее, меняя исходные данные в программе и считывая ответ.
На рис. 4.1 показан Mathcad-документ, проводящий линию по заданным точкам и рассчитывающий значение функции в промежуточной точке. В документе через радиокнопки можно выбирать шесть из многих предназначенных для этой цели встроенных функций Mathcad. Три типичных и самых популярных метода на рис. 4.1 отображены графиками: линейная интерполяция (linterp), интерполяция сплайнами (*spline)[1] и аппроксимация полином n-й степени (regress).
В промежуточной точке (t = 15°С) оценки значения функции[2], рассчитанные разными методами, естественно, несколько отличаются. Это должно нас не очень смущать, т.к. значения функции в самих опорных точках (0, 10, 18, 25, 34 и 50°С) не являются «истиной в последней инстанции» и даны с погрешностью, связанной с методикой проведения замеров в экспериментах.
Коснемся функций, перечисленных на рис. 4.1. На рис. 4.2 показана скрытая область, хранящая эти операторы. Пять функций из шести работают в паре с универсальной встроенной функцией интерполяции interp. Функция regress возвращает вектор, первые три элемента которого – это служебная информация, а остальные – коэффициенты полинома, который формируется через оператор суммы.
Особо следует отметить задачу экстраполяции. Перечисленные на рис. 4.1 функции за исключением loess, могут выдавать ответ при аргументе, выходящем за рамки исходных значений. На рис. 4.2 показано, что оператор t := if(…), запрещает это делать, меняя значение аргумента с числа на строку “Error”. Ограничения здесь бывают двух сортов – физические и статистические, если так можно выразиться: ниже нуля градуса Цельсия мы опускаться не можем – вода превратиться в лед[3], график правее 50°С показан пунктиром. Тут (50°С < t < 60°С) cчет по формулам блокируется, но характер зависимости (экстраполяцию) можно видеть. Кстати, в среде Mathcad для этой задачи (предсказание поведения функции «в будущем») есть специальная встроенная функция predict, работу которой можно испробовать на сайте http://…
В Mathcad встроено еще ряд функций для решения задач аппроксимации, содержащих в своем названии корень fit (от англ. fitting — сглаживание): expfit, lgsfit, logfit, pwrfit, sinfit, genfit, linfit, lnfit и medfit. Их описание можно найти на сайте http://twt.mpei.ac.ru/ochkov/Mathcad_12/Гурский_Турбина_Вычисления_в_Маткад_Глава_16.pdf. А работу некоторых из этих функций можно увидеть на сайте http://twt.mpei.ac.ru/mas/worksheets/regress.mcd.
Задача поиска коэффициентов выражения, наилучшим образом описывающего
точки на плоскости – это типичная оптимизационная задача. На рис. 4.3
отображен WebSheet,
решающий данную задачу для произвольного выражения с одним, двумя или тремя
коэффициентами. На рис. 4.4 показано ядро этой программы. В данной задаче очень
важно задать правильные (удачные) значение первых приближений для коэффициентов
а,
Главный недостаток интерполяции состоит в том, что промах по единственной точке (неверные исходные данные или ошибка при их ручном вводе) сразу сильно искажает результат. Также сильно искажает результат преднамеренный или невольный выход за рамки интервала интерполяции. При аппроксимации же искажение будет не таким сильным, что дает основание считать аппроксимацию в какой-то мере и экстраполирующим решением, решением, сглаживающим последствия возможных промахов. С другой стороны, промах при вводе исходных данных легко видеть и затем исправить, если к этим данным приложена интерполяция, а не аппроксимация. Тут важно выводить на графики и кривые и сами исходные точки – см. рис. 4.1 и 4.3.
Второй недостаток интерполяции – из нее трудно вычленить вид искомой функции для ее использования, например, в другой программной среде, чего не скажешь про аппроксимацию, когда эта функция обычно задается заранее, а в расчете только уточняются ее коэффициенты
Вернемся к нашему графику зависимости скорости звука в морской воде от ее температуры и солености – см. рис. 3.10. Нам остается рассказать о последнем шаге «оцифровки» семейства кривых (поверхности) – о превращении таблицы в формулу, по которой и строится (воссоздается) график.
Встроенные в Mathcad функции lspline, pspline и cspline, о которых речь шла выше, годятся для сплайн-интерполяции табличных зависимостей функций и двух аргументов. В этом случае первым аргументом функций *spline и interp должен быть не вектор, а матрица из двух столбцов, первый из которых хранит дискретные значения первого аргумента (X), а второй – (У). Второй же аргумент функций *spline и interp остается вектором и будет хранить дискретные значения функции аргументов Х и У в соответствующих стоках. Но работа с этими функциями при двух аргументах не совсем удобна из-за того, что эти функции требуют «квадратных» исходных табличных данных, где число точек по первому аргументу должно быть равным числу точек по второму аргументу. Реальные же данные на графиках, как правило, ложатся некую в прямоугольную (далеко не квадратную) таблицу. Из-за этого приходится либо искусственно «оквадрачивать» исходные табличные данные: убирать лишние строки или столбцы, либо наоборот – добавлять недостающие с помощью той же интерполяции или аппроксимации, например. Альтренативное решение данной проблемы на рис. 4.5, где показана программа, сводящая задачу об интерполяции по двум аргументам к двум задачам по интерполяции по одному аргументу. Для этого из исходной матрицы-таблицы М с помощью функции submatrix вычленяются два вектора Х (боковик таблицы) и У («шапка» таблицы) и матрица Z (значение функции). Далее первой интерполяцией в цикл for формулируется вспомогательный вектор Z’ «вставляющий» в матрицу М недостающий промежуточный столбец. Второй интерполяцией в данном столбце определяется нужная искомая точка.
Второе преимущество двойной интерполяции (рис. 4.5) в том, что в ней можно комбинировать различные методы, обходя при этом разные «подводные камни». Так на рис. 4.6 показан случай так называемой осцилляции при решении конкретной интерполяционной задачи[4], когда промежуточная кривая (пунктир) «качается» (от лат. oscillo – качаюсь) – дает нежелательный «всплеск» (см. рис. 4.6а). «Выправить» задачу удалось за счет замены на втором шаге интерполяции сплайном (cspline) на интерполяцию линией (linterp) – см. рис. 4.6б. Решение, конечно, можно считать, несколько грубоватым – ломаная линия имеетступенчатую первую производную. Но не стоит забывать о том, что сама количественная оценка значений аргументов будущей функции (шероховатость!) также весьма грубая. Такие замены (сочетания разных методов интерполяции и аппроксимации по разным аргументам в разных диапазонах) бывают весьма продуктивны, если этим занимается сам специалист-прикладник в конкретной области науки и техники. Здесь нужно помнить (или забывать) о том, что сами исходные кривые были получены с довольно значительной погрешностью (5-10-20 и более процентов). В этом контексте спор о преимуществах тех или иных способов проведения линий через точки или вблизи них теряет свою остроту: по таким линиям можно судить о неких качественных сторонах тех или иных физических явлений, но считать по ним весьма проблематично. Другое дело специальные номограммы, предназначенные именно для расчетов, или собственно таблицы в справочниках, требующие интерполяции или аппроксимации – см. рис. 4.8 и 4.9. На этих рисунках отображено решение типовой задачи из справочника – дана таблица различны параметров вещества (турбинного масла в данном случае) от температуры (см. матрицу Data на рис. 4.9 и комментарии к ней: t – температура, Ср – удельная изобарная теплоемкость, λ – теплопроводность, ν – кинематическая вязкость, μ – динамическая вязкость, а – температуропроводность и Pr – число Прандтля[5]). Необходимо: а) рассчитать значение параметров в промежуточных точках; в) показать на графике характер изменения выбранного параметра от температуры и с) выводить на дисплей выражение, аппроксимирующее табличные данные.
Для решения данной задачи задействована функция linfit, возвращающая коэффициенты аппроксимирующего выражения, заданного в виде функции-вектора F(x) с компонентами – элементарными функциями и их комбинациями. В нашем примере (рис. 4.8 и 4.9) семь параметров турбинного масла разбиваются на две группы с линейной (а0 + а1·х) и нелинейной зависимостью (а0 + а1·х + а2/х + а3/х2) в заданном диапазоне температур 5 – 100° С. Оригинал таблицы можно увидеть здесь.
Другая типичная организация табличных данных показана на рис. 4.10 (оригинал таблицы можно увидеть здесь). Тут параметров не семь, а только два – один в виде скаляра, а другой – в виде матрицы с двумя строками. Но эти параметры относятся не к одному, а к нескольким материалам (вещества). В этом случае табличные данные целессобразно хранить не матрицей (рис. 4.9), а программой (рис. 4.10), возвращающей составной массив – матрицу-строку с двумя элементами, первый из которых скаляр, а второй – матрица новая (вложенная) с двумя строками и переменным числом столбцов.
Еще один типичный случай организации табличных данных можно увидеть на рис 4.11, а его графическое отображение на рис. 4.12 (исходный же график можно увидеть здесь). Это неполностью заполненная таблица или разреженная матрица. Такая матрица в среде Mathcad недопустимы. Отсутствующие элементы в матрице Data на рис. 4.10 – это переменные типа NaN (Not a Number) с белым цветом шрифта (невидимая переменная) – на рис. 4.10 она имела привычный черный цвет в строке «Дихлорметан». Недостаток двумерного регрессионного анализа, использованного на рис. 4.10, нами уже отмечался – это приложение одной степени полинома к двум аргументам, что иногда имеет нежелательные последствия (см. рис. 4.6). На рис. 4.13 этот недостаток устранен за счет уже использовавшегося нами приема замены одного двухмерного инструмента на два одномерных.
Наш подход к снижению порядка задачи, показанный на рис. 4.5, 4.6 и 4.13 можно применить и для решения уже трехмерной задачи — генерации с помощью сплайн-интерполяции функции трех аргументов (рис. 4.14).
Для хранения исходных табличных данных четырехмерной задачи нужна уже не плоская, а объемная матрица со
строками, столбцами и слоями. В языках программирования такая
конструкция есть (трехмерный массив: A(i, j, k)
— BASIC, A[i, j, k] — Pascal и C), но в
Mathcad — нет. Один из выходов — использование составного массива,
[1] lspline строит сплайн с линейным продолжением, pspline – с параболическим, а cspline – с кубическим, но внутри заданных точек значения, возвращаемые этими тремя функциями почти не отличаются.
[2] У нас это электропроводимость чистой воды, дискретные значения которой при разных температурах введены в расчет через таблицу Mathcad-командой Insert/Data/Table. Эти данные можно также считывать с диска.
[3] По умолчанию считается, что таблица, приведенная на рис. 4.1, дана для атмосферного давления.
[4] Строится зависимость влияния шероховатость поверхности внутри трубы и числа Re на гидравлическое сопротивление.
[5] Некоторые данные в столбцах таблицы являются функциями данных других столбцов, и об этом часто забывают.