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

Гибридом декартова графика (рис. 1.17 и 1.18) и графика поверхности (рис. 1.21) является так называемый трехмерный точечный график (

Элементы матрицы (или вектора) можно уподобить столбикам, расставить их на плоской поверхности двух аргументов и получить трехмерную столбчатую диаграмму (

О графических возможностях пакета Mathcad можно написать целую книгу (что входит в перспективные планы автора). Мы же отсылаем читателя к приложению 1, где перечислены некоторые команды настройки графиков. Видов графиков, как уже было отмечено, семь, а инструментов работы с ними – три. Первый инструмент – это форматирование графика. Объемные графики (вернее, псевдообъемные – видимость объема имитируется на все том же плоском дисплее) – это скорее способ украшения
документов, а не средство научного анализа. Настоящая работа ведется с плоскими (двухмерными) графиками. Для них в среде Mathcad припасены два других инструмента: лупа (Zoom – рис. 1.26) и трассировка (Trace – рис. 1.27).
Если человеку необходимо более подробно рассмотреть какой-либо объект, он берет в руки лупу ¾ она и изображена на одной из кнопок панели Graph (см. рис. 1.26). Нажав на эту кнопку, пользователь вызывает диалоговое окно X-Y Zoom с четырьмя полями и пятью кнопками. Протяжкой на рассматриваемом графике можно выделить прямоугольную область (см. рис. 1.26), отмеченную пунктиром. Координаты области отображаются в диалоговом окне X-Y Zoom. Нажатием кнопки Zoom эту выделенную область можно увеличить до размеров исходного графика ¾ и таким образом под лупой рассмотреть левый корень системы уравнений, решаемой на рис. 1.8. Можно все вернуть в исходное положение (кнопка Unzoom), выбрать новую область и раскрыть ее окончательно (Full View).
На график можно посмотреть как бы через прицел снайперской винтовки, если нажать на последнюю, девятую кнопку графической панели (рис. 1.27), появится диалоговое окно X-Y Trace. Если после этого курсором мыши прощупывать кривую на графике, то появятся два волоска (crosshair), координаты пересечения которых отображаются в окне Trace. Эти координаты можно скопировать (CopyX и CopyY) в буфер обмена, а потом перенести в Mathcad-документ и использовать в качестве первого приближения при поиске корня системы уравнений.
Естественно, семь графиков далеко не всегда могут удовлетворить потребности и фантазию пользователя. Разработчики Mathcad не стали утяжелять пакет «экзотической» графикой, а создали новый пакет Axum, предназначенный для более тонкой и изощренной визуализации данных из среды Mathcad. Кроме того, данные из Mathcad-документа с помощью инструмента MathConnex (см. приложение 7) можно перенести в среду Excel или MatLab и там воспользоваться графическими возможностями этих пакетов.
Встроенные операторы вводятся через нажатие кнопок с их изображением – см. рис. 1.3.
Функции (встроенные и пользовательские) с одним или двумя аргументами можно вводить в Mathcad-документ через нажатие кнопок «fx», «xf», «xfy» и «xfy» (нижний ряд на панели Evaluation на рис. 1.3). При этом появляются заготовки постфиксного
и префиксного операторов с одним операндом («fx» и «xf» – первый из них и был задействован на рис. 1.16), инфиксного «xfy» и древовидного «xfy» операторов с двумя операндами. Древовидный вызов оператора хорошо проиллюстрирован на рис. 3.14 в этюде 3.
Встроенные операторы отличаются от встроенных функций (не только в среде Mathcad, но и в программировании), во-первых, тем, что функции все равны между собой, оператор же умножения «главнее» оператора сложения: 2+2*2 равно шести, а не восьми. Кроме того, в среде Mathcad встроенную функцию с одним или двумя аргументами можно непосредственно вызвать в виде оператора (a sin, например). Со встроенными операторами не все так просто. На рис 1.29 даны иллюстрации работы функций и операторов в среде Mathcad.
Второе отличие функции от оператора в среде Mathcad в том, что оператор имеет фиксированное количество операндов: один (n! – факториал, например), два (сложение, вычитание, степень, дифференциал, неопределенный интеграл), три (сумма и произведение элементов вектора, производная высокого порядка и др.), четыре (сумма и произведение ряда, определенный интеграл (см. выше) и т.д.).
Некоторые же функции (Find, MinErr, Minimize, Maximize, например) способны иметь дело с переменным числом аргументов. Третье отличие функции от оператора в среде Mathcad в том, что встроенную функцию можно переопределить. Если, например, пользователя не устраивает то, что аргумент синуса должен быть в радианах, он может заставить синус «глотать» угловые градусы (рис. 1.29).
Внешне для пользователя функция отличается от оператора тем, что у функции есть имя (это обычно слово или сокращение слова), а у оператора – символ. Правда, некоторые операторы вообще не имеют ни имени, ни символа: xy (x в степени y) и Vi (i-й элемент вектора V). Функций-анонимов в среде Mathcad нет.
Пара операторов Mathcad может иметь один и тот же символ, но прописанный разным стилем: сравните светлое равно (вывод числового значения) и полужирное равно (булево равенство). В среде Mathcad можно оперировать одноименными переменными и функциями пользователя с различным шрифтом, с помощью которого отмечаются совершенно различные переменные и функции. Так, в расчете основные переменные можно прописать шрифтом размером 14, а вспомогательные – 10. Традиционные языки программирования такого «безобразия» не допускают.
Одна из причин популярности Mathcad заключается в том, что пользователь вправе вставлять в документы либо функцию, либо оператор в зависимости от того, к чему он привык, изучая математику в школе или в институте. Благодаря этому Mathcad-документ максимально похож на лист с математическими выкладками, написанными от руки или созданными в среде какого-либо текстового процессора (Scientific Word, ChiWriter и др.).

Между числовыми и текстовыми переменными нет глубокого водораздела. Это проиллюстрировано в пункте 2 на рис. 1.30, где моделируется бросание монетки: орел – это 1, а решка – это решка. Переменная a принимает то числовое, то текстовое значение. Об этом мы еще поговорим в этюде 6 (рис. 6.46 и 6.47). В связи с этим в среду Mathcad 8 введены специальные булевы функции, возвращающие 0 или 1 – в зависимости от типа аргумента (см. пункт 4 на рис. 1.30).
[1] Покойная матушка автора умела вычислять на счетах квадратный корень.
[2]
На рисунках книги шрифт переменных и констант Arial Cyr (см. рис. 1.29), а шрифт комментариев – Times New Roman Cyr.
[3]
Там может быть и оператор, например Vi или M<n>, локализующий элемент массива (V) или матрица (M), куда заносится соответствующее значение.
[4]
Проблему русских имен переменных в среде Mathcad мы рассмотрим в главе 6 при раскладке пасьянса.
[5]
Между Given и Find могут быть записаны и неравенства. Но об этом позже.
[6]
Элементом вектора (матрицы) может быть новый вектор (матрица). В этом случае говорят о вложенном массиве (nested arrays).
[7] В математике принято говорить «элемент матрицы» и «компонента вектора», а не «элемент вектора». Но мы будем применять второй термин, так как в среде Mathcad нет принципиальной разницы между вектором и матрицей: вектор ¾ это матрица с одним столбцом.
[8]
Если при создании матрицы нажимать не OK, а Insert, то окно работы с матрицами будет оставаться на экране дисплея.
[9] В разработке функций, предназначенных для решения алгебраических уравнений и систем (Find, MinErr и др.) принимала участие и фирма Frontline System, Inc. (см. конец приложения 1 с указанием авторских прав). Эта же фирма поставила Решатель (Solver) для электронных таблиц Excel.
[10]
Это не совсем так, вернее совсем не так – см. начало этюда 3.
[11]
Символьная математика Mathcad умеет не только выяснять, является ли выражение полиномом или нет, но и вычисляет коэффициенты полинома (см. этюд 7).
Требуется найти угол вырезки a (альфа), при котором объем ведра будет максимальным.
Эту оптимизационную задачу можно решить аналитически: см. пункт 3.3 на рис. 2.1. Но, как понимает читатель, далеко не всякую математическую задачу можно решить аналитически. Иначе бы не было таких научных дисциплин, как «Прикладная математика», «Программирование», «Численные методы» и др. Поэтому мы рассмотрим численное
решение задачи.
И при численном, и при аналитическом решении задачи мы должны вывести зависимость объема ведра V от угла вырезки a. Далее при аналитическом решении можно взять первую производную от этой функции, приравнять ее к нулю и найти корень полученного уравнения. Не обойтись тут и без второй производной, если нужно убедиться, что найденное решение – максимум, а не минимум или точка перегиба, где, как помнит читатель из курса матанализа, первая производная также равна нулю. Жестянщик, которому поручат сделать пожарное ведро, скорее всего, незнаком с дифференциальным счислением, азы которого мы только что изложили. Но в среде Mathcad поставленная задача вполне окажется по плечу «компьютеризированному» жестянщику.
Рисунок заготовки ведра и самого ведра в пункте 1 на рис. 2.1
сделан с помощью графического редактора Paintbrush и перенесен в Mathcad-документ через буфер обмена Clipboard.
В пункт 2 на рис. 2.1 скопированы данные о геометрии конуса из стандартного справочника Mathcad, который удобен тем, что входит в состав пакета и всегда находится под рукой. Справочник открывается командой Open Book… в меню Help. Перенос данных из справочника в Mathcad-документ также автоматизирован, что исключает их искажение – списывая формулу из книги немудрено и ошибиться.
Пункт 3 на рис. 2.1 – это «воспоминание о будущем». Там записаны операторы символьных (аналитических) преобразований (см. этюд 7). В пункте 3.1 оператором solve решается алгебраическое уравнение, позволяющее сформулировать функции пользователя с именами r, h и V (см. пункт 3.2). Зависимости выводятся из несложной геометрии круга и конуса: длина дуги выкройки (2×p×R-2×p×R×a/360) становится длиной окружности в основании конуса (2×p×r), а высота конуса h, радиус его основания r и радиус заготовки R ¾ это стороны прямоугольного треугольника, длины которых связаны теоремой Пифагора (см. пункт 3.2 на рис. 2.1). Оператор substitude позволяет заменять в выражениях подвыражение на другое и вывести еще одну формулу V(a) без вызова вспомогательных функций r(a) и h(a) – см. конец пункта 3.2 на рис. 2.1.
В пункте 3.3 берется производная от V(a), которая сразу упрощается (оператором символьного преобразования simplify). Несложный анализ вида производной показывает, что оптимальный угол вырезки – это один из корней квадратного уравнения a2-720a+43200. Но мы пока на это не обращаем внимания и начинаем численное решение задачи.
Пусть радиус заготовки R равен одному метру (см. начало рис. 2.2). Это можно было бы и не оговаривать[1], так как значение оптимального угла вырезки не зависит от радиуса заготовки (см. конец рис. 2.1), но для численного
решения задачи о пожарном ведре это необходимо. Правда, можно было написать проще – R:=1, не привязываясь к метрам (литрам, галлонам – см. ниже). Но единицы физических величин позволяют нам дополнительно вести контроль правильности формул через соответствие размерностей[2].
Далее записана (вернее, скопирована из рис. 2.1) цепочка функций пользователя, формирующая нашу анализируемую функцию V(a), в которую вложены другие функции – r(a) и h(a). В последнюю в свою очередь также вложена функция r(a). Механизм вложения функций и операторов (встроенных и пользовательских) ¾ это мощный инструмент не только Mathcad, но и других программных сред, позволяющий быстро и изящно решать довольно сложные задач. Вложенные функции просты по виду, а механизм их формирования открыт, чего не скажешь о функции V(a) в пункте 3.2 на рис. 2.1.
Прежде чем искать максимум функции, необходимо убедиться, что он есть. Лучший же способ увидеть максимум – просмотреть график функции. В среде Mathcad есть семь видов графиков (см. этюд 1), первый из которых (X-Y-график в декартовых координатах) отображен на рис. 2.2. Здесь график построен по «двухшаговой» технологии: задается вид функции и сразу отдается команда на вставку графика в Mathcad-документ. По умолчанию аргумент меняется от минус 10 до плюс 10 с пятьюдесятью точками на графике. После построения наброска графика его нужно будет отформатировать – изменить разброс аргумента и др.
На графике в районе 60-70 градусов отчетливо виден максимум функции. Как его уточнить?
Для решения такой задачи в Mathcad 8 встроена новая функция Maximize, возвращающая координаты максимума анализируемой функции вблизи точки начального приближения. Если из заготовки вырезать сектор с углом в 66 с чем-то градусов, то такое ведро будет иметь максимальную вместимость ¾ 403 литра (106 с половиной американских или 88 с половиной британских галлонов – продолжение темы единиц измерения физических величин, начатой в этюде 1).
На рис. 2.3 показано численное решение «двухведерной задачи» в среде Mathcad. В основном оно повторяет решение, показанное на рис. 2.2, но имеет такие особенности:
Сделано допущение, что радиус заготовки равен единице – переменная R в расчете отсутствует.
Функция объема ведра не опирается на вложенные функции. Это несколько усложняет понимание сути задачи, но ускоряет расчет.
Объем второго ведра рассчитывается также через функцию V, но аргумент a при этом сдвинут на 360 градусов: второму ведру достаются обрезки от первого ведра.
На первом графике выведена не одна кривая, а три: объем первого ведра V(a), объем второго ведра V(360 - a) и сумма объемов обоих ведер SV(a)[4]. Кроме того, шкала оси функции начата не с нуля, а со значения 0.38 для того, чтобы пользователь отчетливо увидел два максимума у функции SV и выиграл пари: заготовку нужно разрезать не по диаметру, а несколько иначе, чтобы получить два разных ведра, но с максимальным суммарным объемом.
Построен график производной функции SV, на котором видны три точки пересечения кривой с осью абсцисс, свидетельствующих о двух симметричных максимумах и об одном локальном минимуме в середине. При форматировании первый график был обрамлен (умолчание), а на втором прорисованы оси X и Y. Второй график строится намного дольше первого из-за того, что значение производной в каждой точке графика приходится высчитывать, используя алгоритм численного дифференцирования. А это сама по себе довольно сложная задача. Можно, конечно, из рис. 2.1 скопировать в рис. 2.2 выражение для производной и работать уже с ней, но мы договорились решить задачу только численными методами.
Максимумы для разнообразия найдены не через функцию Maximize, как на рис. 2.2, а через поиск корней производной функции SV. Для этого в расчет включена встроенная функция root, возвращающая корень уравнения и тоже требующая первого приближения к решению. Изменили первое приближение с 120 на 240 угловых градусов ¾ и ответ иной (второй). Кроме того, пришлось изменить c 10-3 на 10-6
У квадратной жестянки по углам вырезаются четыре квадрата. Полученная таким образом крестообразная заготовка сгибается по пунктиру в прямоугольную призму без верхней крышки, а четыре шва свариваются (паяются). Требуется рассчитать размер сторон вырезаемых квадратов (a ¾ отношение размеров квадратов), при котором объем нашего «квадратного ведра» (коробки) будет максимальным.
На рис. 2.4 показано численное решение задачи. Оно отличается от решения по пожарному ведру (см. рис. 2.2) только видом анализируемой функции и методом оптимизации: использована функция MinErr[6] в паре с ключевым словом Given, между которыми булево равенство и ограничение, заставляющее систему искать локальный, а не глобальный максимум.
Несложный анализ функции V(a) показывает, что она имеет максимум при a=1/6. Это позволяет оценить точность использованного нами метода численной оптимизации.
Продолжение задачи о коробке также похоже на продолжение задачи о пожарном ведре: обрезки идут на изготовление новых четырех коробок, новые обрезки (их уже будет 16) тоже пойдут в дело, и так до бесконечности. Но до нее (до бесконечности) мы доберемся только в этюде 7, сейчас же рассмотрим только первые семь шагов раскроя квадратной заготовки, приводящих к формированию 5461 (1+4+16+43+44+45+46) коробок-«матрешек» (рис. 2.5).


На рис. 2.8 помещен протокол решения «трехведерной» задачи. Поиск максимума начат опять же с формирования функции пользователя и с ее графического анализа. Поверхность функции двух переменных строилась так, как показано в пункте 3 на рис. 2.8. Создается сетка с 1681 узлом: 41 линия (0 до 40) по оси переменной a и столько же по оси переменной b. Эта сетка кладется на плоскость a-b, а затем ее узлы (от 0 до 360 градусов с шагом 9 градусов) поднимаются по оси функции V на «подобающую» каждому узлу высоту. Далее эта ажурная конструкция с помощью диалогового окна форматирования поворачивается[10]
вокруг осей V, a и b так, чтобы человек смог увидеть то, что ему нужно, – максимумы, минимумы и др. Всю эту работу машина берет на себя. От человека требуется только наметить узлы сетки (i := 0.. 40, j := 0.. 40), дать «угловое» значение ее узлам (ai := 9×i, bj
:= 9×j – углы меняются от 0 до 360 с шагом 9), заполнить матрицу M значениями «трехведерной» функции ¾ Mi,j:= V(ai, bj), перевести курсор на свободное место и отдать команду построения трехмерного графика, отображающего элементы матрицы М.
Основной недостаток трехмерной графики Mathcad заключается в том, что область изменения аргументов должна быть прямоугольной. Но в нашей «трехведерной» задаче она треугольная, так как аргументы связаны ограничением a+b£360. В пункте 2 рис. 2.8 функция V строится так, чтобы ее значения, выходящие за рамки треугольника, приравнивались к нулю (метод штрафных санкций[11]). Из-за этого задняя грань поверхности на рис. 2.8 получилась зубчатой. Тем не менее видны максимумы на сторонах треугольника области существования аргументов a и b и провисание в центре. В трех вертикальных сечениях просматривается «двухведерный» верхний график из рис. 2.3. Пакет Mathcad не смог решить двухмерную «трехведерную» оптимизационную задачу по методике, представленной на рис. 2.3 (взятие частных производных по переменным a и b и поиск корней полученной системы алгебраических уравнений). Пакет Mathcad пытался искать максимум «у края обрыва» и «сваливался» в него. Не справилась с этой задачей и специально введенная в Mathcad 8 функция Maximize при всех трех начальных приближениях (см. пункты *.1[12]): мы «танцевали» к максимуму от трех «печек» ¾ из центра треугольника (120 и 120 градусов ¾ см. пункт 4), из одного угла треугольника (0 и 0 градусов ¾ см. пункт 5) и от одной из сторон треугольника (150 и 210 градусов ¾ см. пункт 6). С поиском максимума справилась «старая добрая» функция MinErr[13] ¾ см. пункты *.2. Старая в буквальном смысле ¾ она ведет свою родословную еще с DOS-версий Mathcad и поэтому хорошо отлажена.
Данная задача относится к широкому классу задач под названием задачи линейного программирования: необходимо установить план (программу!) выпуска изделий (у нас это стулья), ориентируясь на целевую функцию (у нас их две ¾ общее количество и общая стоимость стульев) и принимая во внимание ограничения
(ресурсы по доскам, ткани и человеко-часам). Из рис. 2.9 видно, что можно выпускать не более 150 стульев (максимум первой целевой функции). А вот с максимизацией их общей цены не получилось: Mathcad не умеет решать задачу целочисленного линейного программирования. План выпуска стульев, максимизирующий их количество, вышел целочисленным случайно.
Спрашивается, как нужно организовать перевозки (найти значения переменных с1м1, с1м2, с2м1 и с2м2), чтобы затраты были минимальны. На рис. 2.10 дан ответ. Парадокс задачи в том, что по самому дешевому маршруту (со второго склада в первый магазин – 800 у.е.) ничего не возится (с2м1 = 0). Этот парадокс мы также обыграем в следующем этюде.
Задачи на рис. 2.9 и.2.10 простенькие, но очень, если так можно выразиться, жизненно важные. На каждом шагу приходится что-то оптимизировать (расходы, например), принимая во внимание всякого рода ограничения (доходы!). Возвращаясь к сноске 17, можно привести такой пример. После часа пик (зимнее утро, к примеру) расход электроэнергии падает и необходимо снижать нагрузку электрогенераторов. Как это делать? Можно отключить отдельные турбогенераторы, а можно оставить их в работе, изменив нагрузку. Диспетчер энергосистемы дает соответствующие команды, ориентируясь на некие целевые функции: средний расход топлива по системе, выброс с дымовыми газами вредных веществ в атмосферу, износ оборудования, степень готовности электростанций и дальше менять нагрузку и т.д. Переменные такой оптимизации могут быть и вещественными (мощность отдельного энергоблока, которая меняется, естественно, в разумных пределах, определяемых техническими условиями – ограничения в задаче) и целочисленными (число работающих блоков). Эта задача очень сложная, но и очень эффективная – здесь речь идет о высвобождаемых составах с топливом.
Вот еще примеры. Когда нужно убирать пшеницу? Пораньше – зерно еще не вызрело. Попозже – часть зерна уже осыпалась. Сколько и каких акций стоит купить на ограниченную сумму денег, чтобы будущий дивиденд был максимален? В каких средствах массовой информации стоит размещать рекламу на выделенные по смете деньги, чтобы эффект от нее был максимален?
Разговор об оптимизации мы продолжим в этюде 3 в несколько ином ключе.
[1] Тем более что, переменная R уже занята под хранение градусов Ренкина. В этом можно убедиться, набрав R= и получив R=0.556 K (градус Ренкина в градусах Кельвина). Присвоением R:=1m мы «испортили» данную системную переменную, что может выйти нам боком, если в расчет придется вводить температуру. Отсюда вытекает хорошее правило работы в среде Mathcad: «Никогда не пользуйтесь оператором «:=». Для присваивания значения переменной лучше работать с оператором «=», автоматически превращающимся в оператор «:=», если соответствующая переменная не занята пользователем или системой.
Теперь мы подшутим
над функцией root, и если не получим от этого удовольствия, то хотя бы уясним себе, какие «подводные камни» могут нас ожидать.
Если в качестве первого приближения (опорной точки) принять не минус 50, а плюс 5 (пункт 5), то функция root выкинет «белый флаг»: сообщение «Can’t converge to a solution», отказываясь решать поставленную задачу, хотя плюс 5 намного ближе к корню, чем минус 50. Вот вам и первое приближение! Но это еще полбеды. Настоящая беда случается тогда, когда функция root (как, впрочем, и некоторые другие функции и операторы Mathcad) не отказывается решать поставленную задачу, выдавая при этом неверный результат (феномен медвежьей услуги).
Если функцию y(x) умножить на константу, например на 10-5, то ее корни останутся на старых местах. Но это утверждение в среде Mathcad не является истиной – см. пункт 6 рис. 3.1.
Не то беда, Авдей Флюгарин,
Что родом ты не русский барин,
Что на Парнасе ты цыган...
Беда, что скучен твой роман.
Не то беда, что ты, Mathcad, неправильно решил простейшую задачу: и более мощные специализированные пакеты, ориентированные только на решение алгебраических уравнений и систем, делают тут промах. Беда в том, что Mathcad в этом честно не признался, как это было в пункте 5 на рис. 3.1.
Философский смысл любой шутки заключается в том, что у шутника и у того, над кем подшучивают, разные понятия о природе вещей. У пользователя и у среды Mathcad разные понятия о том, что такое корень уравнения: человек считает, что корень – это то значение аргумента, при котором выражение равно нулю; функция же root «считает», что корень – это то значение аргумента, при котором значение выражения по модулю не превышает значения системной переменной TOL, которая по умолчанию равна 10-3. Отсюда и путаница в пункте 6 на рис. 3.1. Чтобы функция root там сработала правильно, необходимо переменной TOL присвоить новое значение (10-7, например), заменив им предопределенное. Можно поступить и по-другому – умножить в пункте 6 правую часть функции на 100000, убрав тем самым коэффициент 0.00001. Метод балластных (нормирующих) коэффициентов особенно эффективен при решении алгебраических систем (что уже было отмечено в этюде 1). Он позволяет уравнять все уравнения (прошу простить за тавтологию) по отношению к точности, с которой система решается через блок Given-Find. На нижнем графике рисунка 3.1 ось X «утолщена» до значения TOL (пунктир). Корень там, где кривая касается этой «толстой» оси.
Проанализировав эту программу, можно понять, почему корнем другого уравнения Y = X2 - 4 при нулевом (симметричном – проблема Буриданова осла) начальном приближении оказывается плюс, а не минус 2:
x := 0 x := root(X2 - 4, X) x = 2
Хотя в BASIC-программе полностью реализован метод секущих, описанный в Руководстве пользователя, но...
Во-первых, в вышеприведенную BASIC-программу не заложено никаких сообщений об ошибках, а функция root его выдает (см. пункт 5 рис. 3.1). Во-вторых, BASIC-программа, в отличие от функции root, прекрасно нашла корень уравнения y(x), заданного в пункте 1, от начального приближения, равного плюс 5. И, в-третьих, в BASIC-программу заложен не метод секущих, а комбинированный метод Ньютона-секущих. Классический метод секущих требует задания не одного, а двух начальных значений аргумента: через одну точку проводится касательная (метод Ньютона), а через две – секущая. Вот так мы разобрались с функцией root! Хотели над ней подшутить, а она сама над нами посмеялась.
В этюде 6 читатель может найти функции пользователя для поиска корней уравнений методом Ньютона, методом половинного деления и методом секущих, ориентированные по точности не на значение анализируемой функции (Do ... Loop Until Abs(FnY(X1)) <= TOL – см. программу на рис. 3.2), а на значение аргумента (Do ... Loop Until Abs(X1 – X0) <= TOL, например). На такую подмену приходится идти, решая реальную задачу. Так, например, в коллекции автора есть учебная программа определения значения рН раствора, где рН ¾ это корень уравнения электронейтральности раствора (баланс катионов и анионов). Но нас мало интересует дисбаланс ионов (отклонение анализируемой функции от нуля) ¾ главное, чтобы новое значение рН отличалось от предыдущего не более чем на величину 10-3
(обычная точность рН-метра).
Двухмерная экспоненциальная функция имеет минимум (нулевое значение) при x=1 и y=10. Это хорошо видно в пункте 1 на рис. 3.3, где показаны поверхность и линии уровня вблизи минимума[2]. Эти графики несложно построить, если, конечно, знаешь, где находится минимум, охватываемый переменными x1, x2, y1 и y2. Поверхность развернута так (см. ее координаты на фрагменте окна редактирования), чтобы линии уровня являлись проекцией поверхности на горизонтальную плоскость.
В пункте 2 при поиске минимума в качестве начальных приближений давались точки, расположенные по углам прямоугольника ¾ области существования аргументов на графиках. Так испытывалась сходимость метода поиска минимума. Зафиксирована всего лишь одна осечка: при начальном приближении x=1.2 и y=14 не был найден корень системы, состоящей из приравненных к нулю частных производных двухмерной экспоненциальной функции[3].
Функция Розенброка примечательна тем, что ее минимум невозможно увидеть ни на линиях уровня, ни на поверхности (пункт 1 на рис. 3.4). Скажем осторожнее – автору не удалось этого сделать. На графиках просматривается типичный овраг, где анализируемая функция в одном направлении изменяется круто, а в перпендикулярном – слабо. Этим пытаются как бы дезориентировать программу поиска минимума – на то и создаются тестовые функции. Графики строились по новой, третьей технологии ¾ задавался центр (x и y) квадрата со стороной 2D области существования аргументов на графиках. Вторая технология ¾ это когда задаются координаты углов прямоугольника существования аргументов на графиках (см. пункт 1 на рис. 3.3). Первая технология (кстати, самая неудобная для понимания) была использована, например, в рис 2.8 – там область размаха поверхности завуалирована в значениях переменных i и j.
В пункте 2 на рис. 3.4 поиск минимума велся от трех начальных точек (10-10, 100-10 и 100-100) тремя способами. Итоги «соревнования»: на первом месте по-прежнему функция MinErr, которая выдавала абсолютно точный результат (пару единиц). На втором месте функция Minimize c относительно правильными ответами. Функция же Find сошла с трети дистанции – только при первом приближении был выдан относительно точный результат. Здесь, по-видимому, с функцией Find овраг сыграл злую шутку.
Функция Розенброка имеет минимум (нуль) при x=1 и y=1. Через эту точку можно провести секущие плоскости и показать на декартовых графиках, что функция там минимальна (см. пункт 3 на рис. 3.4), ее частные производные равны нулю, а частные производные второго порядка положительны.
Функция Пауэла (рис. 3.5) имеет четыре аргумента. Следовательно, в трехмерном мире – даже виртуальном – никакими поверхностями и линиями уровня минимум (нуль) функции Пауэла не локализовать (не визуализировать). Проверить соответствие найденной точки (четыре нуля) минимуму можно либо через декартовы графики сечений по технологии рисунка 3.3), либо через линии уровня, что и было сделано в пунктах 2.1 и 2.2 на рис. 3.5. Для этого две координаты фиксируются на нулевых значениях, а две другие изменяются вокруг нуля (нули – это координаты минимума). Таких топограмм четырехмерной функции можно построить шесть штук следующих пар аргументов: a0-a1, a0-a2 (см. рис. 3.5), a0-a3, a1-a2, a1-a3 и a2-a3 (по четырем последним парам читателю предлагается построить топограммы самому). Решение по функции Пауэла более-менее удачно велось с помощью функций MinErr и Minimize. Поиск корня системы четырех алгебраических уравнений – частных производных функции Пауэла – здесь применить нельзя, так как переменные анализируемой функции у нас не скалярные величины (с именами a, b, c и d, например), а одиночная переменная-вектор с именем a. По индексным переменным производная в среде Mathcad не высчитывается. Читатель при желании может переопределить функцию Пауэла (задействовать в ней четыре переменные-скаляра) и решить задачу с помощью частных производных. Так, кстати, решалась эта задача в предыдущих изданиях книги.
Выводы по испытаниям трех функций – в конце этюда.
Есть более простой способ статистического расчета числа p, чем тот, который использовал Бюффон. Можно нарисовать квадрат, вписать в него круг («квадратура круга») и бросать туда камешки. Так как в площади круга запрятано число p («пи эр в квадрате» – вспомним старый анекдот о том, почему у поезда колеса стучат), то через подсчет попаданий в круг можно оценить число p. Бюффон этот метод не использовал наверное из-за того, что трудно добиться равномерного попадания камешков в квадрат.
На рис. 3.6 зафиксирована оценка площади круга, вписанного в квадрат, и сделана попытка расчета значения числа p. Как это делалось?
В пункте 1 на рис 3.6 формируются два вектора X и Y, элементы которых (а их 7000) ¾ случайные вещественные числа в интервале от минус до плюс единицы. X и Y – это по своей сути координаты точек случайного падения камешков в квадрат размером 2 на 2. Далее в пункте 2 эти камешки сортируются на «чистых и нечистых»: координаты точек, попавших в круг, дублируются в векторах Xo и Yo (o ¾ попал), а не попавших ¾ в векторах Xx и Yx (x ¾ промах). После этого несложно визуализировать попадание точек в круг с помощью параметрического декартового графика. Достаточно при форматировании графика указать, что линий нет, а есть одни точки. При этом точки, попавшие в круг (вектора Xo и Yo), более толстые. Теперь для оценки числа p можно подсчитать число попаданий камешков в круг (пункт 4) ¾ число ненулевых элементов вектора Xo (или Yo). В пункте 5 число p рассчитывается и сравнивается с его точным значением.
Функция rnd в среде Mathcad имеет аргумент, отличающийся от своих аналогов на языках программирования. В среде BASIC или Pascal аргумент функции rnd связан с приставкой псевдо- в ее названии, а не с диапазоном генерируемых случайных чисел. Да, числа генерируются случайные, но ряд этих чисел псевдослучаен, так как его в любой момент можно повторить. За такой повтор на языках программирования отвечает аргумент функции rnd, а в среде Mathcad – число (по умолчанию это 1), хранящееся в окне Seed value for random numbers (инициализация генератора случайных чисел) ярлыка Built-in Variables окна Math Options, которое вызывается на дисплей командой Options в меню Math:

На рис. 3.7 не только рассчитан объема конуса методом Монте-Карло, но и дано изображение самого конуса. В изобразительном искусстве есть направление, называемое пуантилизм
(от французского pointiller – писать точками). Его последователи формируют картину отдельными мазками правильной точечной или прямоугольной формы. Мы уже нарисовали в такой манере «Черный круг в черном квадрате». «Натюрморт с конусом», выполненный в этом же стиле, можно увидеть в пункте 3 рис. 3.7. Здесь «поработал» Scatter Plot – трехмерный точечный график. Так незаметно мы перешли от двухмерной к трехмерной графике – от живописи к… скульптуре.
Рисовался (вернее, ваялся) конус так. В прямоугольный параллелепипед размером 2 на 2 на 1 бросались случайным образом «камешки» – формировались три вектора X, Y и Z, элементы которых – случайные числа в диапазоне минус единица – плюс единица (векторы X и Y) и ноль – единица (вектор Z см. пункт 1 на рис. 3.7). Далее точки, не вписывающиеся в конус, отсеивались: соответствующим элементам вектора Z присваивалось значение бесконечности (пункт 2), и они... «улетали на небо», формируя линию над конусом (пункт 2).
Пункты 3 и 4 – это подсчет объема конуса методом Монте-Карло и сравнение его с истинным объемом конуса[6]. Сомнения остались, но очень незначительные – в пределах ошибки метода Монте-Карло.
Один великий художник сказал, что ваять очень просто – берется глыба мрамора и от нее отсекается все лишнее. «Всякий талант неизъясним. Каким образом ваятель в куске каррарского мрамора видит сокрытого Юпитера и выводит его на свет, резцом и молотом раздробляя его оболочку?»[7]
А вот еще одно парадоксальное высказывание: «Играть на фортепиано очень просто. Достаточно в нужный момент нажимать нужные клавиши!» Все это можно рассматривать как некую браваду или как попытки Художника подурачить публику, пристающую с «глупыми» вопросами о тайне творчества. А можно принять за руководство к действию. С живописью мы разобрались (см. рис. 3.6) – возьмемся теперь за скульптуру.
На рис. 3.7 у нас есть «глыба мрамора» размером 2 на 2 на 1, от которой также отсекается лишнее (осколки падают не на пол, а улетают на небо). Пока у нас «изваялся» простой конус, но ничто не мешает нам усложнить формулы в пункте 2 на рис. 3.7 и слепить Венеру Милосскую или Мыслителя Родена… На рис. 3.8 дается методика ваяния человеческой фигуры.
Автор пока поступил проще. Он через команды форматирования графика окрасил точки конуса на рис. 3.7 в зеленый цвет, подбросил еще несколько точек красного цвета и большего диаметра. Поучилась неплохая новогодняя елка.
Сейчас компьютер широко используется как рабочий инструмент художника (интеллектуальная кисть или что-то в этом роде). Распечатки цветных принтеров оправляются в рамки и выставляются в, так сказать, реальных и виртуальных компьтерно-художественных салонах (см., например, журнал «КомпьюАрт»).
Но автору хотелось бы обратить внимание уважаемых читателей на другое – на проблему эстетического вида не просто компьютерных рисунков, а листингов программ и, в частности, на проблему соответствия (или противопоставления) формы листинга содержанию программы.
Программисты, которым не чуждо образное мышление, давно уже подметили, что процедуры и функции имеют свое собственное «лицо», по которому она безошибочно узнается на экране дисплея или на бумаге принтера. Одна процедура как ухоженная крестьянская лошадка круглая и гладкая – работает себе спокойно, перекачивая, например, данные из одного формата в другой. И внешне она неприметна – взгляд на ней не останавливается. Другая процедура все время норовит выкинуть какой-нибудь фортель, настолько она неотлаженна (необъезженна). И своими очертаниями она походит на скакуна, в седле которого сидит герой многочисленных живописных полотен и скульптур. Третья процедура так и просится, чтобы ее оправили в раму и повесили на стену, настолько она хороша и закончена, а главное, ее форма полностью отвечает ее содержанию. Она передает не только мысли, но даже и настроение
художника, пардон, программиста, ее создавшего.
Автор далеко не искусствовед и не смеет особо распространяться на эту тему.
В комментариях к публикуемым компьютерным рисункам, как правило, подчеркивается, что их авторы – компьютерные художники. В прилагательных к существительным очень часто таится некая ущербность или по, крайней мере, двусмысленность: не просто математика, а «Прикладная математика» – см. главку «Mathcad и Maple» в этюде 7. Термин компьютерный художник
В решении задачи о пожарных ведрах на рис. 3.9 в пунктах 1 (одноведерная задача) и пункт 2 (двухведерная задача) формируются векторы a, V и V2 по 361 элементу в каждом[9]. Ключевой оператор решения использует встроенную функцию max, возвращающую максимальное значение своего аргумента – вектора (матрицы). Функции, возвращающей номер максимального элемента вектора (матрицы), в среде Mathcad нет – с рассуждениями по этому поводу читатель может ознакомиться в этюде 4. Поэтому на рис. 3.9 (как и на рис. 3.6 и 3.7) задействован оператор суммы, перебирающий все варианты ведер и запоминающий угол вырезки для изготовления одного ведра максимального объема – 66 градусов (пункт 1) или двух ведер с максимальным суммарным объемом – 117 и 243 (360-117) градусов (пункт 2). Наш метод перебора опасен тем, что если у вектора (матрицы) два и более максимальных значений (как у вектора V2), то в лучшем случае появится сообщение об ошибке, а в худшем – неправильный ответ (см. пункт 2 на рис. 3.9 с aопт[10]=360). При решении трехведерной задачи (пункт 3) на область максимума была как бы наброшена сетка, в узлах которой просчитаны значения «трехведерной» функции. Далее были определены координаты сетки, где данная функция максимальна. Такой контрольный расчет перебором еще раз показал, что третье ведро лишнее – мы получили уточненное решение двухведерной задачи из пункта 2.
Наше решение выглядит несколько извращенным – в функцию max, составляющую ядро расчета на рис. 3.9, и в другие, ей подобные, также заложен перебор: что-то другое здесь вряд ли придумаешь.
Хотя как сказать. Представим такую житейскую ситуацию. Садовод собрал на своем дачном участке урожай яблок и решил похвастаться самым крупном
плодом перед соседями. Будет ли он перебирать все яблоки, чтобы выбрать предмет гордости? Конечно, нет. Самое большое яблоко никому не нужно. Нас интересует самое большое яблоко с определенной степенью вероятности. Кроме того, в термин «большое» мы вкладываем не вес яблока и не его геометрические размеры, а его зрительный образ.
Оказывается, при переборе всех вариантов выпуска стульев (а их не так уж много – на рис. 3.10 мы просчитали 150 на 150 = 22 500 вариантов) можно найти два оптимальных плана, причем самый оптимальный и по цене, и по количеству (20+112=132 стула стоимостью 1504 у.е.) – это не округление дробного ответа, полученного на рис. 2.9. Возвращаясь к теме враждебности задачи, можно так подобрать ее условия, что ответ окажется совсем вдалеке от первоначального дробного…
Это был эпиграф, приступаем к рассказу.
У Михаила Жванецкого часто спрашивают, откуда он берет темы для своих миниатюр. «Выглядываю в окно и прислушиваюсь к разговорам на улице», – таков ответ великого сатирика. «А как Вы все это запоминаете?», – следует новый вопрос. «Да я забыть не могу!»
Житейские сюжеты стоит коллекционировать и для написания компьютерных этюдов, что является хобби автора этой книги.
По профессии автор – преподаватель вуза (Московского энергетического института), где он читает курс лекций по информатике и смежным дисциплинам, а также руководит группой технологов и программистов, разрабатывающих обучающие программы и компьютерные тренажеры для ТЭС и АЭС[12]. Электростанциям и энергообъединениям нужны наши программы, но их приобретению мешает пресловутый кризис неплатежей. Вот какой компьютерный этюд имел место в марте 1997 года.
Акционерное общество Тамбовэнерго, не имея свободных денег, тем не менее изъявило желание приобрести наши программы. Котовскому лакокрасочному заводу (ЛКЗ – Тамбовская обл.) для производства нужна электроэнергия. Московскому энергетическому институту для ремонта аудиторий требуется краска. Нашей научной группе необходимо новое компьютерное «железо», инструментальные средства и, естественно, зарплата. Для решения подобных проблем человечество еще на заре цивилизации придумало деньги[13]. Переход же нашей страны от непонятно чего к рынку возродил натуральный обмен – бартер[14]. В вышеописанной товарной цепочке не хватало одного звена, чтобы она замкнулась. К счастью, в МЭИ поступила партия компьютеров, парочку которых мы договорились обменять на краску. В этой комбинации заключается только часть описываемого компьютерного
Протокол «контрольного взвешивания» краски в среде Mathcad приведен на рис. 3.11. Комментарии поясняют, что происходит в формулах. Во-первых, функция Maximize, как и ожидалось, дала дробный ответ (см. пункт 2) – маленьких банок можно не брать, если можно брать дробное количество больших. Пришлось, вспомнив эпиграф и название этюда, перейти к перебору вариантов. В Mathcad-документе формируются две матрицы с именами Об (пункт 3.2) и Ст (пункт 3.3), элементы которых (их 1088 – у матриц 17 столбцов и 64 строки) хранят значения объема (Об) и стоимости (Ст) краски в зависимости от комбинаций расфасовки. Далее (пункт 3.4) некоторым элементам матрицы Об присваиваются нулевые значения, если данные комбинации расфасовки не проходят по стоимости. Остальное – ловкость рук и никакой математики: в пункте 3.5 определяется номер строки (переменная N_15) и номер столбца (N_55) матрицы Об, на пересечении которых находится элемент с максимальным значением. Ответ (6 маленьких барабанов и 15 больших) неприятно удивил Олю. Она невольно обманывала меня на 5 литров краски и на 139 тыс. руб.
Метод поиска координат точки максимума, реализованный на рис. 3.11 (двойная сумма), имеет существенное ограничение: в анализируемой матрице (у нас это Об) должен быть только один максимальный элемент. Если их несколько, то ответ будет неверен: в переменные N_15 и N_55 будут записаны суммы координат точек с максимальным элементом. Мы это наблюдали в пункте 2 на рис 3.9.
Так Mathcad сэкономил мне почти полторы сотни тысяч рублей. Деньги не такие уж большие, но если присовокупить к ним новый компьютерный этюд в книгу, новую тему лекции и новую лабораторную работу по информатике, а также гонорар за эту книгу, то игра стоила свеч.
Вернувшись из Тамбова домой в Москву, я в спокойной обстановке у своего родного компьютера еще раз проанализировал задачу. И вот что получилось.
Во-первых, заставить Решатель Excel правильно «разъяснять» задачу о краске можно было, изменив начальные установки Решателя. А для этого нужно было не полениться и нажать на кнопку Параметры... в диалоговом окне Поиск решения. В новом диалоговом окне Параметры поиска решения достаточно было допустимое отклонение уменьшить с 5 до 1%. После этого правильное решение (15 больших и 6 маленьких барабанов) было бы найдено. Честно говоря, в Excel плох не Решатель, а его начальные установки. Очень мало пользователей Excel, прибегающих к услугам Решателя, нажимают кнопку Параметры... Тот же, кто разбирается в сути установок оптимизации, как правило, с Excel не работает. Отсюда и недоразумения.
|
Вариант расфасовки (число маленьких барабанов/число больших барабанов) |
2/16 |
6/15 |
13/13 |
37/6 |
|
Объем краски (л) |
910 |
915 |
910 |
885 |
|
Остаток невыбранных денег (руб.) |
186 000 |
47 000 |
12 000 |
11 000 |
На рис. 3.12 помещен протокол решения задачи о поиске места для ларька на дачном участке. Критерий поиска – минимум суммы расстояний от ларька до всех остальных домиков. Их координаты X и Y – случайные числа в интервале от 0 до 20 (наша задача учебная). Ларек может быть либо встроен в один из домиков (пункт 1), либо стоять отдельно (пункт 2). На рис. 3.12 координаты встроенного ларька определяются перебором. Затем эти координаты (Xiопт и Yiопт) становятся первым приближением для уточнения местоположения отдельно стоящего ларька. Заканчивается расчет графическим описанием и сути и решения задачи. Здесь главное – правильно отформатировать точки на графике. Поэтому выведено окно форматирования графика.
Вернемся к тестовым задачам на рис. 3.3-3.5.
На рис. 3.3 и 3.4 формируются матрицы М, элементы которых – значения анализируемых функций в узлах сетки, накрывающей точку минимума. Элементы матрицы М средствами Mathcad превращаются в линии уровня и в поверхность. Но эти матрицы могут сослужить нам и другую службу – координаты их минимальных элементов могут стать точками первого приближения к минимуму. Нащупать минимум (максимум) функции более чем двух аргументов (например, функции Пауэла – см. рис. 3.5) можно средствами программирования (см. этюд 6).
Транспортная задача, решенная нами на рис. 2.10, кочует из одного учебника в другой. Везде отмечается такой ее парадокс – по самому дешевому маршруту при минимизации затрат на перевозки ничего не возится. Если бы мы решали эту задачу вручную без компьютера, то сначала полностью загрузили бы дешевый маршрут, а потом все остальные. Этим парадоксом может воспользоваться хозяин транспортного предприятия, максимизировав свою прибыль от перевозок (см. рис. 3.13):
В этом случае второй по дороговизне маршрут (1200 у.е.) остается свободным – внешне все выглядит прилично.
Но парадокс задачи не в этом. Вернее, не только в этом. Дело в том, что она недостойна не только функций Minimize и Maximize, но даже и грубого перебора, так как сводится к решению простейшего уравнения:
с1м1+с2м1=40
Если с1м1=0, то с2м1=40, и затраты на перевозки максимальны. Если с2м1=0, то с1м1=40, и затраты максимальны. Рис. 2.10 и 3.13 – это чистой воды извращения. Или неумение либо нежелание подумать как следует над задачей.
Во-первых, булевы функции в среде Mathcad можно определить. Математика (см., например, «Справочник по математике для научных работников и инженеров» Корн Г. и Корн Т.) оперирует одной
одноместной (с одним аргументом) и семью
двухместными (с двумя аргументами) булевыми функциями. Все они определены в пунктах 0-7 на рис. 3.14. Если булеву переменную уподобить выключателю с двумя позициями («вкл» и «выкл»), то конъюнкция
– это последовательное соединение выключателей (пункт), а дизъюнкция – параллельное. Электрический аналог эквиваленции (пункт 3) может очень пригодиться для освещения длинного коридора, свет в котором зажигается и тушится независимо в двух его концах двумя выключателями.
В пункте 8 на рис. 3.14 сформирована трехместная булева функция с именем Решение, возвращающая вердикт жюри из трех человек: решение проходит, когда «за» голосуют двое или трое. Воздерживаться или уклоняться от голосования нельзя.
Функция Решение (программная реализация процедуры голосования) в пункте 8 на рис. 3.14 также имеет электрический аналог (аппаратная реализация – машинка для голосования) – комбинацию выключателей, соединенных последовательно и параллельно.
Двухместные булевы функции (пункты 1-7 на рис 3.14) имеют четыре (22) комбинации значений аргументов (таблица истинности), трехместные – уже 8 (23), одноместная, естественно, только две (21) – 0 или 1. Самих же двухместных булевых функций может быть 16 (42 – мы описали только семь), трехместных уже 64 (43 – мы описали только одну). Одноместных булевых функций четыре (41 – мы описали только одну). Вот другие три одноместные булевы «функции» y2(x):=1, y3(x):=0 и y4(x):=x. Но никакой практической ценности в программировании они не имеют: первые две (y2 и y3) – это константы, а y4 – это просто сам аргумент. Ненаписанные нами остальные девять (16-7) двухместные булевы функции (там тоже есть константы) имен не имеют и, как правило, ни в математике, ни в программировании не применяются.
В математике булева функция выдает два значения (0 – 1, да – нет, истина – ложь и т.д.), в программировании же – минимум три: да, нет и... аварийный останов, связанный с ошибкой (один или несколько аргументов не определены). Такую ошибку можно обработать (в Mathcad для этого служит оператор on error) и пустить расчет по третьему пути. Одноместных булевых функций может быть больше четырех. Как понравится такая функция: y5(x):=if(rnd(1)>0.7, 1, 0), возвращающая единицу с вероятностью 30%, и нуль – с вероятностью 30%.
Что происходит в функции Победитель?
В начале дуэли все участники живы: все три элемента вектора Статус принимают значение “жив”[50]. Далее проводится жеребьевка: определяется направление очередности выстрелов (если переменная Очередь равна единице, то очередность идет в таком направлении ...0®1®2®0®1®2.., если минус единице – ...0®2®1®0®2®1) и определяется первый стреляющий (переменная Стрелок). Кроме того, обнуляется переменная Убийство, по которой прерывается цикл выстрелов в дуэли.
Математическая модель дуэли опирается на цикл с выходом из середины (while ... break…): дуэль продолжается, пока не будут сделаны два результативных выстрела. В теле цикла while определяется Цель – самый меткий противник, которого убивают (СтатусЦель
¬ “убит”), если, во-первых, не промахиваются (МеткостьСтрелок
> rnd(1)) и (And), во-вторых, не (Not) стреляют намеренно в воздух. Второе имеет место при хитрой тактике стреляющего (ТактикаСтрелок = 2) и (And), если метких противников более одного.
Определение следующего стреляющего ведется в цикле с постпроверкой (while ... break): цикл прерывается, когда, перебирая очередь, отмеченную выше (...0®1®2®0®1®2.. или ...0®2®1®0®2®1...), «натыкаются» на живого участника.
Возвращает функция Победитель номер участника дуэли (0, 1 или 2), оставшегося в одиночестве (значение переменной Стрелок по выходу из цикла).
Функция Победитель возвращает непредсказуемое целочисленное значение 0, 1 или 2, так как в ней в трех местах вызывается встроенная в Mathcad функция rnd, которая возвращает псевдослучайное число в интервале от нуля до значения аргумента функции rnd. Этот аргумент у нас равен либо единице (случайный выбор очередности выстрелов и имитация выстрела с вероятностью попадания, пропорциональной меткости стреляющего), либо трем (случайный выбор первого стреляющего – здесь дополнительно работает встроенная функция floor, возвращающая у положительного вещественного числа его «пол» (в смысле не «потолок» – по-английски a floor): floor(0.54), floor(1.82), floor(2.48) = 0...
В пункте 1 на рис. 4.1 переменным X и Y присваиваются транспонированные векторы-строки, а не просто векторы-столбцы. Это делается для компактности записи. Кроме того, элементы векторов имеют лишние нули – 4.0 вместо 4 и т.д. За счет этого выравниваются по вертикали пары значений Xi и Yi. Без такой маленькой хитрости рано или поздно пары собьются, что будет мешать их просмотру и редактированию. Альтернативное решение этой проблемы – хранение пар данных в матрице с двумя строками и с числом столбцов, равным числу пар[1].
Считается, что программиста от простого смертного можно отличить по простому тесту. Если программиста поставить в голову шеренги и приказать: «По порядку рассчитайсь!», то программист сначала уточнит, по какой системе нужно рассчитываться (двоичная, восьмеричная, шестнадцатеричная, десятеричная[2]...), а потом выкрикнет: «Нулевой!» В среде Mathcad по умолчанию номер первого элемента вектора (первого ряда и первого столбца матрицы) нулевой. Именно поэтому при семи экспериментальных точках, координаты которых заносятся в векторы X и Y, константа N равна шести (феномен программиста в строю). Номер первого элемента массивов и векторов хранится в системной переменной ORIGIN (an origin – начало, источник), значение которой (по умолчанию оно нулевое) в Mathcad-документе можно изменять (ORIGIN:=1, например). Допустимо менять и второе умолчание «шеренги» – систему счислений.
Ввести в среде Mathcad переменную-вектор можно двумя различными способами (см. этюд 1): отдачей команды Matrices из меню Math (Insert – Mathcad 7 и 8) либо нажатием на панели математических инструментов кнопки с изображением матрицы (щелкнув по ней курсором мыши) – см. рис. 1.7. Ввод за переменной ее индекса также допустим двумя способами: нажатием на панели математических инструментов на кнопку-иероглиф «Переменная с индексом» или набором за именем переменной символа открывающихся квадратных скобок (рудимент языков Pascal и C, где квадратные скобки означают индексную переменную).
Векторы X и Y совсем не обязательно вводить в Mathcad-документ вручную с клавиатуры. Если экспериментальный стенд оборудован средствами АСНИ (автоматизированной системой научных исследований) и данные с приборов заносятся на магнитный диск, то Mathcad-выражение X:=READPRN(имя файла) поможет считать их и оформить в виде Mathcad-вектора (матрицы) с именем X[3]. Кроме того, не следует забывать, что Mathcad – это полноценное Windows-приложение со встроенными средствами обмена в статике и динамике (Clipboard, DDE, OLE). Объемную задачу можно решить лишь тогда, когда голос Mathcad звучит в стройном хоре других приложений (графические, текстовые и табличные процессоры, базы данных, языки программирования и т.д.).
Задача аппроксимации, и не только линейной – это типичная оптимизационная задача (см. этюды 2 и 3 и линии уровня в пункте 3.2 на рис. 4.2), сводящаяся к поиску минимума целевой функции СКО (среднеквадратичное отклонение) двух переменных a и b. В свою очередь, линейное сглаживание сводится к решению системы линейных алгебраических уравнений (см. этюд 1 и пункт 3.3 на рис. 4.2), состоящей из приравненных к нулю частных производных функции СКО по коэффициентам a и b.
Решение, показанное в пункте 3.2, предпочтительней для сферы образования – оно как бы «кричит» о сути метода наименьших квадратов: в функции СКО фигурирует квадрат, а сама функция минимизируется.
Найденные тем или иным способом значения коэффициентов a и b сглаживающей функции y(x) = a + b× x позволяют построить на графике прямую с роящимися вокруг нее точками (у нас квадратиками – рис. 4.3). Подобным графиком на практике, как правило, завершают обработку экспериментальных данных: график, во-первых, даст наглядное представление о качестве сглаживания, а во-вторых, поможет в случае чего отловить допущенные ошибки ввода исходных данных (пропуск десятичной точки, например). Этой цели может служить и предварительная сортировка векторов (см. пункт 2 на рис. 4.1): ошибочные значения (промах эксперимента, неправильный ввод данных) часто всплывают на концах упорядоченного вектора. В-третьих, график сам по себе ценен. С помощью графика, то есть с другого конца, можно довольно быстро решить задачу линейного сглаживания. У автора в лаборатории есть сотрудница, у которой глаз-алмаз: она при помощи тонкой прозрачной линейки так проводит прямую вблизи экспериментальных точек, что по ней можно определить коэффициенты a и b с точностью не меньше трех процентов (толщина карандашной линии).
Несколько слов о графических возможностях Mathcad и других подобных пакетов. Если студент начнет строить график функции по технологии, заложенной в математическом пакете, то автор выгонит такого студента с занятий, да притом вослед будет улюлюкать и топать ногами. Все (скажем осторожнее – почти все) математические пакеты при построении графиков никак (почти никак) не используют элементы искусственного интеллекта, а просто сканируют значение аргумента и проставляют точки с заданным пользователем шагом или с шагом, определяемым разрешением дисплея (принтера). Студентов же учат совсем другому – анализу функции, поиску характерных точек (корней, минимумов, максимумов, точек перегиба и т.д.), опираясь на которые и строится график: парабола, гипербола или какая-нибудь там лемниската Бернулли (см. рис. 1.18 в этюде 1). Но Богу – Богово, кесарю – кесарево, машине – машиново. Беда многих преподавателей в том, что они относятся к математическим пакетам не как к инструментальным средствам, требующим определенной сноровки и навыка (и головы, конечно), а как к своим своеобразным коллегам, которых на пушечный выстрел нельзя подпускать к студентам, изучающим математику, дабы они (студенты) не набрались от них (от пакетов) разных глупостей. Дежурная фраза одной знакомой автора: «Mathcad – круглый дурак, Maple – полный кретин, Mathematica – законченная идиотка».
На рис. 4.4 формируются матрица A коэффициентов при неизвестных и вектор B свободных членов системы четырех линейных алгебраических уравнений, к которой сводится задача об аппроксимирующем полиноме третьей степени (кубическом). Сама система может решаться либо матричными операторами (a:=A-1×B), либо обращением к встроенной функции lsolve (см. пункт 3.3.2 рис. 4.2). Полноценный программный пакет всегда должен предоставлять пользователю возможность сделать какую-либо операцию двумя, а еще лучше тремя различными способами. Пользователь будет чувствовать себя комфортно в программной среде, если его право выбора не ограничено. На выбор же могут влиять не только сомнения в качестве тех или иных готовых инструментальных средств (см. выше) или нечеткое знание области их применения, но и вкусовые предпочтения пользователя. Решение системы линейных алгебраических уравнений (базовая задача любого Решателя) через матричные операторы неудобно тем, что пользователь все время напарывается на сообщения об ошибках, набирая R:= B / A или R := B × A-1. Крамольность первого выражения очевидна – в математике нет понятия матричного деления. Но почему нельзя писать R := B × A-1, знают далеко не все. Так или иначе, о KISS-принципе не следует забывать. Помнить следует и другое предостережение, зафиксированное в пословице: «Простота хуже воровства». Выкладки рис. 4.4, от которых рябит в глазах из-за обилия «сигм», можно упростить вводом переменных с индексом Ak,m и Bk (рис. 4.9). Но! Автор уже раз обжегся на суммировании в среде Mathcad, когда рассматривал систему искусственного интеллекта SmartMath (рис. 7.15 в этюде 7): за суммой, особенно тройной с пересекающимися индексами, может таиться ошибка.
Аргументы функции linfit – числовые векторы X и Y и вектор-функция f, хранящая набор пользовательских сомножителей (на рис. 4.5 – это 1, 1/x, 1/x2 и 1/x3), коэффициенты которых функция linfit возвращает в вектор a, опираясь на метод наименьших квадратов (см. рис. 4.5).
1. В среде Mathcad есть функция genfit (GENeral FITting – общее сглаживание), допускающая множество параметрических коэффициентов у сомножителей произвольного вида. Она, в отличие от функции linfit, решает нелинейную[7]
задачу аппроксимации. Пример использования функции genfit дан в пункте 7.1 на рис. 4.6. В функцию-вектор f необходимо занести не только само сглаживающее выражение (у нас это экспоненциальная функция), но и его частные производные[8] по искомым параметрам a0 и a1,.объединенным в вектор. Из-за того, что нелинейная задача может иметь более одного решения, функция genfit требует начального приближения (третий аргумент-вектор), вблизи которого и ищется одно из решений.
2. Без функции genfit задачу нелинейной аппроксимации, как правило, решают линеаризацией
исходных данных, подключая к работе уже известные нам функции intercept и slope. В нашем случае (пункт 7.2 на рис. 4.6) при исходном уравнении Y=a × x b встроенные функции intercept и slope прикладывают к новому уравнению Ln(Y)=Ln(a)+b×Ln(X). Но это приводит к некоторому искажению математической модели: сумма X и сумма логарифмов X – это далеко не одно и то же.
3. В ряде случаев при решении задачи нелинейной аппроксимации можно обойтись и без функции genfit (пункт 7.1 рис. 4.6), и без искажающей линеаризации (пункт 7.2), применив открытый алгоритм поиска значений коэффициентов a0 и a1, минимизирующих целевую функцию – сумму квадратов отклонений точек от кривой. Расчет, показанный в пункте 7.3 на рис. 4.6, кроме того, не требует знания частных производных сглаживающего уравнения.
В конце рис. 4.6 приведено графическое сравнение трех методов нелинейного сглаживания. Из графика видно, что линеаризация исходных данных (пункт 7.2) привела к существенным искажениям результата. Поиски же параметрических коэффициентов через функцию genfit (пункт 7.1) и через минимизацию (пункт 7.3) близки по результатам.
Практический совет по выбору наиболее подходящей формулы для сглаживания. В Mathcad-документе можно хранить набор формул, подключая одну из них по мере надобности к статистической обработке экспериментальных точек. В диалоговом окне Properties (свойства – оно вызывается командой Properties в меню Format) есть переключатель Disable Evaluation (Запретить вычисления), позволяющий превращать математические выражения[9]
в комментарий:

Данный переключатель превращает выражение в комментарий (что отмечено черным квадратиком правее выражения) и переопределяет тем самым функцию пользователя. О втором переключателе (Enable Optimization – Разрешить оптимизацию) будет рассказано в этюде 7 (раздел 7.3).
В Mathcad-документе можно записать множество функций пользователя, подключая одну из них к аппроксимации: линейной (рис. 4.2), экспоненциальной (рис. 4.6), логарифмической и т.д., и подбирая тем самым наилучшую модель обработки экспериментальных данных. Выбор формулы можно вести и без переключателя Disable Evaluation, опираясь на свойство среды Mathcad включать в расчеты только последнюю запись функции. Нужную функцию можно просто перетаскивать (технология drag-and-drop) в конец списка.
Отличия в линейной (lspline), параболической (pspline) и кубической (cspline) интерполяции сплайном заметно проявляются только на концах отрезка X1-XN. Эти области на рис. 4.8 рассмотрены «под лупой». Здесь нужно говорить уже не об интерполяции, а об экстраполяции (см. ниже рис. 4.14)
Несложно через точки провести полином N-й степени (рис. 4.9):
Но нанизать опытные точки на интерполяционный «шампур», напрочь игнорируя неизбежные ошибки эксперимента, может только совсем безграмотный исследователь. У интерполяции другие сферы применения. Расскажем об одной из них. При решении в среде Mathcad какой-либо задачи нередко образуется составная функция[10], обращение к которой вызывает длинную цепочку сложных вычислений, связанных с поиском корней уравнения, с дифференцированием, интегрированием и т.д. Работать с такой функцией становится невмоготу даже на мощном компьютере. Один из выходов – омолаживание «бабушки»: табулирование «тормозной» функции с последующей заменой ее на эрзац-функцию, опирающуюся на интерполяцию – линейную, нелинейную или сплайном. И что удивительно, «омоложенная» функция, хоть и теряет напрочь свою физику, но в особых условиях может возвращать более точное значение, чем ее прародительница, в которой накапливаются ошибки численных методов. Узлы же интерполяции можно просчитать на пределе точности.
Далее автор приводит несколько пронумерованных «Советов тем, кто работает с Mathcad», которые публикуются в журнале КомпьютерПресс – на бумаге и на прилагаемом к журналу лазерном диске.
Педагогический опыт автора[12]
говорит о том, что студенты, выполняющие термодинамические расчеты, очень часто ошибаются в размерностях: складывают, например, джоули с британской единицей теплоты, а ответ записывают в калориях, забывая о соответствующем пересчете. Мы этой темы уже подробно коснулись в этюде 1. Среда Mathcad позволяет правильно (с соответствующими пересчетами) оперировать размерными величинами и «ругается» только в крайних случаях, когда, например, складывается длина с массой. Наша пользовательская функция hss(T, P) возвращает размерную величину (удельная энтальпия перегретого водяного пара), а ее аргументы (давление и температура пара) – также размерные величины.
На рис. 4.13 приведен пример формирования с помощью одномерной сплайн-интерполяции еще одной важной термодинамической функции одного аргумента – зависимости удельного объема кипящей воды от ее температуры.
Сплайн-интерполяция особенно критична на границах отрезка (области) существования аргументов. Вне этих границ речь идет уже не об интерполяции, а об экстраполяции – о предсказании значений функции за границами отрезка.
На рис. 4.14 показана работа функции predict, с помощью которой предсказывается
поведение какой-либо зависимости.
У функции predict (prediction – предсказание) три аргумента: вектор имеющихся данных (у нас 100 значений функции y(x) на интервале от 0 до 99 с шагом 1), число последних данных, учитываемых в предсказании (у нас их 50) и число предсказываемых данных (100). На графике на рис. 4.14 толстая кривая – функция y(x) на интервале 0-99, тонкая кривая – функция y(x) на интервале 100-199, а квадратики – предсказанные значения (не все, а только кратные 5), которые где-то в районе 160-170 обрываются.
Задание читателю, подготавливающее его к чтению следующего этюда, где будет описана математическая модель одной финансовой операции, – предсказать курс доллара на будущую неделю, опираясь на данные предыдущего месяца или года.
На рис. 4.15 проиллюстрирована работа сглаживающей функции supsmooth. Решается такая задача – сгладить функцию стоимости, например, комплектующих компьютера (см. задачу на рис. 6.33-6-35) в зависимости от объема закупки (скидка оптовикам).
В пункте 1 на рис. 4.16 оценки экспертов введены в матрицу M, содержащую девять строк (число экспертов n, j=1..9) и десять столбцов (число качеств m, i=1..10).
Пятибалльная средневзвешенная оценка деловых качеств менеджера определяется по формуле:

где m – количество оцениваемых качеств, i := 1, 2.. m,
n – число экспертов, j := 1, 2.. n,
Мj,i – оценка j-м экспертом i-го качества в баллах,
ai – весовой коэффициент для i-го качества.
Весовые коэффициенты определяют относительную значимость качеств: если i-е качество представляется незначимым, то весовой коэффициент ai равен нулю. Придание коэффициенту
ai значения, равного единице, делает незначимыми все остальные качества. Весовые коэффициенты должны удовлетворять следующему условию:

Значения весовых коэффициентов могут быть установлены различными способами:
лицом, принимающим решения (ЛПР);
экспертом с соблюдением условия (4.2);
по результатам экспертных оценок качества, косвенно отражающих мнения экспертов о значимости оцениваемых характеристик.
В последнем случае весовые коэффициенты рассчитываются по формуле:

Вычисленные значения весовых коэффициентов в табличной и графической формах представлены в пункте 3 на рис. 4.17.
Средневзвешенная оценка качества менеджера в баллах (пункт 4 на рис. 4.17) ¾ Kм = 4.162.
Алгоритм определения весовых коэффициентов ai влияет на оценку качества, которая может быть выполнена по выбору ЛПР на основе гипотезы равной значимости весовых коэффициентов:

где ai = 1 / n.
В этом случае оценка качества значительно возрастает (в сравнении со средневзвешенной оценкой) – K1м
= 4.593.
Значения весовых коэффициентов могут быть установлены с учетом специфики конкретного заказа из условия приоритета качеств (требуется «высокий» профессионал, личность приятная во всех отношениях, но без особой склонности к руководству). В этом случае приоритеты могут быть установлены, например, на следующих значениях: для профессиональных качеств ¾ kп1=1.15, личностных ¾ kп2 =1.05, деловых ¾ kп3 =0.8. Значения весовых коэффициентов для групп качеств r определяются формулой:
(ar)i = kпr / n, где n = 9, r= 1, 2, 3,
и равны – (a1)i =
1.15 / 9 = 0.128; (a2)i = 1.05 / 9 = 0.117; (a3)i = 0.8 / 9 = 0.089.
Вектор-строка весовых коэффициентов имеет вид:
| i | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | ||||||||||
| (ar)I | (a1)i | (a1)i | (a1)i | (a2)i | (a2)i | (a2)i | (a3)i | (a3)i | (a3)i |
Средневзвешенная оценка качества менеджера в баллах (пункт 8 на рис. 4.17) ¾
Kм
= 4.496.
Итак, в рассмотренном примере в зависимости от способа определения весовых коэффициентов получены три оценки качества:
Km = 4.162 – при экспертной оценке a;
K1m
= 4.593 – случай равнозначимых весовых коэффициентов;
K2m = 4.496 – при приоритете (kп1=1.15) профессиональных качеств.
Выбор кандидата – прерогатива лица, принимающего решение, или конкурсной комиссии.
Рассмотренная методика применима для оценки качества продукции, результатов испытаний, экзаменов и других задач экспертной оценки.
Возможно, руководителю группы экспертов или заказчику экспертизы потребуется анализ результатов работы экспертов для оценки их компетентности, пристрастий или добросовестности.

Все это словесное описание модели эпидемии легко вмещается в Mathcad-документ (рис. 5.1) с двумя формулами (пункт 2), объединенными в вектор, что эквивалентно BASIC-конструкции:
For t = 1 To 13
Больные(t + 1) = Пр * Больные(t) * Здоровые(t)
Здоровые(t + 1) = Здоровые(t) - Больные(t)
Next
Если бы два выражения в пункте 2 на рис. 5.1 не были заключены в скобки, то их выполнение сразу бы прерывалось сообщением об ошибке. Система Mathcad пыталась бы сначала полностью заполнить вектор Больные, а уже потом – вектор Здоровые. Скобки изменяют порядок счета: он ведется не по строкам, а по столбцам: сначала заполняются вторые элементы векторов Больные и Здоровые (первые элементы заполняются в пункте 1 – начальные условия), а потом третьи и т.д. Этими скобками мы меняем естественный порядок выполнения операторов[1]
¾ они выполняются не слева направо и не сверху вниз, а крест накрест.
Результаты расчета графически отображены[2]
в пункте 3. На тринадцатый день (спад эпидемии) в городе было 105 больных. Критическая точка – девятый день (3972 больных), ради поиска которой и затевают весь этот расчетный сыр-бор: моделируя эпидемию, мы можем распланировать работу санитарных служб города – подвезти в аптеки лекарства, отозвать врачей из отпуска, выписать из больниц выздоравливающих и т.д. В оригинальном решении задачи дополнительно прослеживается динамика изменения числа умерших во время эпидемии. Но мы исключим эту печальную кривую. Тем более что она и математически не вписывается в задачу: динамика числа умерших просчитывается отдельно от динамики больных и здоровых.
В пункте 5 рис. 5.2 столбцы матрицы Z (время, число больных и число здоровых) разнесены по отдельным векторам и отображены графически, что позволяет проследить динамику развития эпидемии.
На рис. 5.2 получены несколько иные результаты, чем на рис. 5.1, хотя характер кривых сохранился: максимум больных (3119) наблюдается не на 9-й, а на 7-й день, на 13-й день мы имеем не 105, а 209 больных. Это объясняется и различными значениями точности расчетов (на рис. 5.1 делалось 13 шагов интегрирования, а на рис 5.2 ¾ 500) и различными примененными методиками (Эйлер против Рунге и Кутта[5]).
В функцию rkfixed заложен широко распространенный метод решения дифференциальных уравнений – метод Рунге ¾ Кутта[6]. Несмотря на то что это не самый быстрый метод, функция rkfixed почти всегда справляется с поставленной задачей. Однако есть случаи, когда лучше использовать более сложные методы. Эти случаи попадают под три широкие категории: система может быть жесткой[7]
(Stiffb, Stiffr), функции системы могут быть гладкими
(Bulstoer) или плавными (Rkadapt). Нередко приходится пробовать на одном дифференциальном уравнении (одной системе) несколько методов, чтобы определить, какой метод лучше (быстрее, точнее). Примерно так мы сравнивали в этюде 3 разные способы поиска оптимумов функции.
Когда известно, что решение гладкое, используется функция Bulstoer, куда заложен метод Булирша ¾ Штёра, а не Рунге ¾ Кутты, используемый функцией rkfixed. В этом случае решение будет точнее. Список аргументов и матрица, получаемая при работе с функцией Bulstoer, такие же, как и при работе с rkfixed.
Можно решить задачу более точно (более быстро), если уменьшать шаг (у нас это Dt) там, где производная меняется быстро, и увеличивать шаг там, где она ведет себя более спокойно. Для этого предусмотрена функция Rkadapt (adaption –адаптация). Но, несмотря на то что при решении дифференциального уравнения функция Rkadapt использует непостоянный шаг, она тем не менее представит ответ для точек, находящихся на одинаковом расстоянии, заданном пользователем. Аргументы и матрица, возвращаемая функцией Rkadapt, такие же, как при rkfixed.
Ответ (51 больной) получен за семь приемов: задается начальное число больных, которое корректируется в зависимости от того, какое число больных оказывается в конце эпидемии. На рис. 5.3 можно, конечно, не дублировать Mathcad-оператры, а просто вручную подправлять первое приближение.
С помощью функций Bustoer, bustoer, Rkadapt, rkadapt, rkfixed, Stiffb, stiffb, Stiffr и stiffr, решающих задачу Коши, краевую задачу можно также решить последовательными приближениями (см. рис. 5.4 с функцией rkfixed):