Если функция f (х) является полиномом, то все его корни можно определить, используя встроенную функцию
polyroots(v),
где v — вектор, составленный из коэффициентов полинома.
Поскольку полином N-й степени имеет ровно N корней (некоторые из них могут быть кратными), вектор v должен состоять из N+I элемента. Результатом действия функции poiyroots является вектор, составленный из N корней рассматриваемого полинома. При этом нет надобности вводить какое-либо начальное приближение, как для функции root. Пример поиска корней полинома четвертой степени иллюстрируется листингом 8.5.
Листинг 8.5. Поиск корня полинома
Коэффициенты рассматриваемого в примере полинома
f (х) = (х-3 ) (х-1)3=х4-6х3 + 12х2-10х+3 (1)
записаны в виде вектора в первой строке листинга. Первым в векторе должен идти свободный член полинома, вторым — коэффициент при х1 и т. д. Соответственно, последним n+1 элементом вектора должен быть коэффициент при старшей степени Xn.
Иногда исходный полином имеется не в развернутом виде, а, например, как произведение нескольких полиномов. В этом случае определить все его коэффициенты можно, выделив его и выбрав в меню Symbolics (Символика) пункт Expand (Разложить). В результате символьный процессор Mathcad сам преобразует полином в нужную форму, пользователю надо будет только корректно ввести ее в аргументы функции polyroots.
Во второй строке листинга 8.5 показано действие функции poiyroots. Обратите внимание, что численный метод вместо двух из трех действительных единичных корней (иными словами, кратного корня 1) выдает два мнимых числа. Однако малая мнимая часть этих корней находится в пределах погрешности, определяемой константой TOL, и не должна вводить пользователей в заблуждение. Просто нужно помнить, что корни полинома могут быть комплексными, и ошибка вычислений может сказываться как на действительной, так и на комплексной части искомого корня.
Для функции poiyroots можно выбрать один из двух численных методов — метод полиномов Лаггера (он установлен по умолчанию) или метод парной матрицы.
Для смены метода:
Вызовите контекстное меню, щелкнув правой кнопкой мыши на слове polyroots. В верхней части контекстного меню выберите либо пункт LaGuerre (Лаг-гера), либо Companion Matrix (Парная матрица). Щелкните вне действия функции polyroots — если включен режим автоматических вычислений, будет произведен пересчет корней полинома в соответствии с вновь выбранным методом.Для того чтобы оставить за Mathcad выбор метода решения, установите флажок AutoSelect (Автоматический выбор), выбрав одноименный пункт в том же самом контекстном меню.
Рассмотрим решение системы N нелинейных уравнений с M неизвестными
f1(x1, ... ,хм) = 0,
. . . (1)
fn(x1, ... ,хм) = 0,
Здесь f1(x1, ... ,хм) , ..., fn(x1, ... ,хм) — некоторые скалярные функции от скалярных переменных хцх2/ ... /хм и, возможно, от еще каких-либо переменных. Уравнений может быть как больше, так и меньше числа переменных. Заметим, что систему (1) можно формально переписать в виде
f(x)=0, (2)
где х — вектор, составленный из переменных x1,х2, ... ,хм, a f (х) — соответствующая векторная функция.
Для решения систем имеется специальный вычислительный блок, состоящий из трех частей, идущих последовательно друг за другом:
Given — ключевое слово; система, записанная логическими операторами в виде равенств и, возможно, неравенств; Find(x1... ,хм) — встроенная функция для решения системы относительно переменных хх,..., хм.Вставлять логические операторы следует, пользуясь панелью инструментов Boolean (Булевы операторы). Если Вы предпочитаете ввод с клавиатуры, помните, что логический знак равенства вводится сочетанием клавиш <Ctrl>+<=>. Блок Given/Find использует для поиска решения итерационные методы, поэтому, как и для функции root, требуется задать начальные значения для всех х1, ... ,xм. Сделать это необходимо до ключевого слова Given. Значение функции Find есть вектор, составленный из решения пс каждой переменной. Таким образом, число элементов вектора равно число аргументов Find.
В листинге 8.6. приведен пример решения системы двух уравнений.
Листинг 8.6. Решение системы уравнений
В первых двух строках листинга вводятся функции, которые определяют систему уравнений. Затем переменным х и у, относительно которых она будет решаться, присваиваются начальные значения. После этого следует ключевое слово Given и два логических оператора, выражающих рассматриваемую систему уравнений. Завершает вычислительный блок функция Find, значение которой присваивается вектору v. Следующая строка показывает содержание вектора v, т. е. решение системы. Первый элемент вектора есть первый аргумент функции Find, второй элемент — ее второй аргумент. В последних двух строках осуществлена проверка правильности решения уравнений.
Часто бывает очень полезно проверить точность решения уравнений, вычислив значения образующих их функций в найденных вычислительным процессором корнях, как это сделано в конце листинга 8.6.
Отметим, что уравнения можно определить непосредственно внутри вычислительного блока. Таким образом, можно не определять заранее функции f (x,y) и д(х,у), как это сделано в первых двух строках листинга 8.6, а сразу написать:
Given
х4 + у2 =3
х+ 2 у = 0
Такая форма представляет уравнения в более привычной и наглядной форме, особенно подходящей для документирования работы.
Графическая интерпретация рассмотренной системы представлена на рис. 8.3. Каждое из уравнений показано на плоскости XY графиком. Первое — сплошной кривой, второе — пунктиром. Поскольку второе уравнение линейное, то оно определяет на плоскости XY прямую. Две точки пересечения кривых соответствуют одновременному выполнению обоих уравнений, т. е. искомым действительным корням системы. Как нетрудно убедиться, в листинге найдено только одно из двух решений — находящееся в правой нижней части графика Чтобы отыскать и второе решение, следует повторить вычисления, изменив начальные значения так, чтобы они лежали ближе к другой точке пересечения графиков, например x=-1, y=-1.
Рис. 8.3. Графическое решение системы двух уравнений
Пока мы рассмотрели пример системы из двух уравнений и таким же числом неизвестных, что встречается наиболее часто. Но число уравнений и неизвестных может и не совпадать. Более того, в вычислительный блок можно добавить дополнительные условия в виде неравенств. Например, введение ограничения на поиск только отрицательных значений х в рассмотренный выше листинг 8.6 приведет к нахождению другого решения, как это показано в листинге 8.7.
Листинг 8.7. Решение системы уравнений и неравенств
Обратите внимание, что, несмотря на те же начальные значения, что и в листинге 8.6, мы получили в листинге 8.7 другой корень. Это произошло именно благодаря введению дополнительного неравенства, которое определено в блоке Given в предпоследней строке листинга 8.7.
Если предпринять попытку решить несовместную систему, Mathcad выдаст сообщение об ошибке, гласящее, что ни одного решения не найдено, и предложение попробовать поменять начальные значения или значение погрешности.
Вычислительный блок использует константу CTOL в качестве погрешности выполнения уравнений, введенных после ключевого слова Given. Например, если CTOL=0.001, то уравнение х=10 будет считаться выполненным и при х=10.001, и при х=9.999. Другая константа TOL определяет условие прекращения итераций численным алгоритмом (см. разд. 8.4). Значение стоъ может быть задано пользователем так же как и TOL, например, CTOL:=0.01. По умолчанию принято, что CTOL=TOL=0.001, но Вы по желанию можете переопределить их.
Особенную осторожность следует соблюдать при решении систем с числом неизвестных большим, чем число уравнений. Например, можно удалить одно из двух уравнений из рассмотренного нами листинга 8.6, попытавшись решить единственное уравнение g(х,у)=о с двумя неизвестными х и у. В такой постановке задача имеет бесконечное множество корней: для любого х и, соответственно, у=-х/2 условие, определяющее единственное уравнение, выполнено. Однако, даже если корней бесконечно много, численный метод будет производить расчеты только до тех пор, пока логические выражения в вычислительном блоке не будут выполнены (в пределах погрешности). После этого итерации будут остановлены и выдано решение. В результате будет найдена всего одна пара значений (х,у), обнаруженная первой.
О том, как найти все решения рассматриваемой задачи, рассказывается в разд. 8.7.
Вычислительным блоком с функцией Find можно найти и корень уравнения с одним неизвестным. Действие Find в этом случае совершенно аналогично уже рассмотренным в данном разделе примерам. Задача поиска корня рассматривается как решение системы, состоящей из одного уравнения. Единственным отличием будет скалярный, а не векторный тип числа, возвращаемого функцией Find. Пример решения уравнения из предыдущего раздела приведен в листинге 8.8.
Листинг 8.8. Поиск корня уравнения с одним неизвестным с помощью функции Find
В чем же отличие приведенного решения от листинга 8.1 с функцией root? Оно состоит в том, что одна и та же задача решена различными численными методами. В данном случае выбор метода не влияет на окончательный результат, но бывают ситуации, когда применение того или иного метода имеет решающее значение.
Если Вы решаете "хорошие" уравнения, как все те, которые были приведены в предыдущих разделах, то можете никогда не задумываться, как именно Mathcad ищет их корни. Однако даже в этом случае полезно представлять, что происходит "за кадром", т. е. какие действия совершаются в промежутке между введением необходимых условий после ключевого слова Given и получением результата после применения функции Find. Это важно хотя бы с позиций выбора начальных значений переменных перед вычислительным блоком. Рассмотрим в данном разделе некоторые особенности численных методов и возможности установки их различных параметров, которые предоставляет Mathcad.
Функция Find реализует градиентные численные методы. Покажем их основную идею на примере уравнения с одним неизвестным f(x)=o для функции f (х)=х2+5х+2, график которой показан на рис. 8.4. Основная идея градиентных методов состоит в последовательных приближениях к истинному решению уравнения, которые вычисляются с помощью производной от f (х). Приведем наиболее простую форму алгоритма, называемого методом Ньютона:
За нулевую итерацию принимается введенное пользователем начальное значение хО=х. В точке хо методом конечных разностей вычисляется производная f'(х0). Пользуясь разложением Тейлора, можно заменить f (х) в окрестности хо касательной — прямой линией f(x)=f(х0)+f '(х0)(х-х0). Определяется точка x1, в которой прямая пересекает ось х (см. рис. 8.4). Если f (x1)<TOL, то итерации прерываются, и значение x1 выдается в качестве решения. В противном случае x1 принимается за новую итерацию, и цикл повторяется: строится касательная к f (х) в точке x1, определяется х2 — точка ее пересечения с осью х и т. д.Рис. 8.4. Иллюстрация метода Ньютона
Модификация алгоритма Ньютона для решения системы нескольких уравнений заключается в линеаризации соответствующих функций многих переменных, т. е. аппроксимации их линейной зависимостью с помощью частных производных. Например, для нулевой итерации в случае системы двух уравнений используются выражения.
Чтобы отыскать точку, соответствующую каждой новой итерации, требуется приравнять оба равенства нулю, т. е. решить на каждом шаге полученную систему линейных уравнений.
Mathcad предлагает три различных вида градиентных методов. Чтобы поменять численный метод:
Щелкните правой кнопкой мыши на названии функции Find. Наведите указатель мыши на пункт Nonlinear (Нелинейный) в контекстном меню. В появившемся подменю (рис. 8.5) выберите один из трех методов: Conjugate Gradient (Сопряженных градиентов), Quasi-Newton (Квази-Нью-тоновский) или Levenberg-Marquardt (Левенберга).Рис. 8.5. Смена численного метода
Чтобы вернуть автоматический выбор типа численного метода, в контекстном меню надо выбрать пункт AutoSelect (Автоматический выбор) Если установлена опция автоматического выбора (о чем говорит флажок, установленный в пункте AutoSelect), то текущий тип численного метода можно узнать, вызвав то же самое подменю и посмотрев, который из них отмечен точкой Два последних метода являются квази-Ньютоновскими, основная идея которых была рассмотрена выше Первый из них, метод сопряженных градиентов, является двухшаговым — для поиска очередной итерации он использует как текущую, так и предыдущую итерации. Алгоритм Левенберга подробно описан в справочной системе Mathcad, а детальную информацию о методах Ньютона и сопряженных градиентов можно найти в большинстве книг по численным методам
Помимо выбора самою метода, имеется возможность устанавливать их некоторые параметры Для этого нужно вызвать с помощью того же контекстного меню диалоговое окно Advanced Options (Дополнительные параметры), выбрав в контекстном меню пункты Nonlinear / Advanced options (Нелинейный / Дополнительные параметры) В этом диалоговом окне (рис. 8.6.) имеется пять групп переключателей, по два в каждой
В первой строке Derivative estimation (Аппроксимация производной) определяется метод вычисления производной Forward (Вперед) или Central (Центральная) Они соответствуют аппроксимации производной либо правой (двухточечная схема "вперед"), либо центральной (трехточечная симметричная схема) конечной разностью
Обратите внимание, что вычисление производной в градиентных численных методах решения уравнений производится более экономичным способом, нежели при численном дифференцировании (см гл 7)
Во второй строке Variable estimation (Аппроксимация переменных) можно определить тип аппроксимации рядом Тейлора Для рассмотренного нами в этом разделе случая аппроксимации касательной прямой линией выберите переключатель Tangent (Касательная), для более точной квадратичной аппроксимации (параболой) выберите Quadratic (Квадратичная) Следующая группа переключателей Linear variable check (Проверка линейности) позволяет в специфических задачах сэкономить время вычислений Если Вы уверены, что нелинейности всех функций, входящих в уравнение, мало сказываются на значениях всех их частных производных, то установите переключатель Yes (Да) В этом случае производные будут приняты равными константам и не будут вычисляться на каждом шаге
С осторожностью изменяйте параметры численных методов Пользуйтесь ими, когда решение не находится при выставленных по умолчанию параметрах или когда расчеты занимают очень продолжительное время.
Рис. 8.6. Диалоговое окно Advanced Options
Пара переключателей Multistart (Сканирование) задает опцию поиска глобального или локального минимума или максимума Если выставлен переключатель Yes (Да), Mathcad будет пытаться найти наиболее глубокий экстремум из области, близкой к начальному приближению Эта опция предназначена, в основном, для настройки (тех же самых, градиентных) алгоритмов поиска экстремума, а не для решения алгебраических уравнений
Наконец, последний переключатель Evolutionary (Эволюционный алгоритм), если установить его в положение Yes (Да), позволяет использовать модификацию численного метода для решения уравнений, определяемых не обязательно гладкими функциями Как мы убедились в этом разделе, все градиентные методы, реализованные в функции Find, требуют многократного вычисления производных Если Вы работаете с достаточно гладкими функциями, то градиентные методы обеспечивают быстрый и надежный поиск корня. Для поиска корня недостаточно гладких функций одной переменной, следует либо выставить данную опцию функции Find, либо использовать метод секущих (функцию root). Помните, что правильный выбор численного метода и его параметров может помочь при решении нестандартной задачи, которая при стандартных установках может и не поддаваться решению.
Иногда приходится заменять задачу определения корней системы уравнений задачей поиска экстремума функции многих переменных. Например, когда невозможно найти решение с помощью функции Find, можно попытаться потребовать вместо точного выполнения уравнений условий минимизировать их невязку. Для этого следует в вычислительном блоке вместо функции Find использовать функцию Minerr, имеющую тот же самый набор параметров. Она также должна находиться в пределах вычислительного блока:
x1:=C1 ... хм: =Cм — начальные значения для неизвестных. Given — ключевое слово. Система алгебраических уравнений и неравенств, записанная логическими операторами. Minerr (x1,... ,хм) — приближенное решение системы относительно переменных х1;... ,хм, минимизирующее невязку системы уравнений.В функции Minerr реализованы те же самые алгоритмы, что и в функции Find, иным является только условие завершения работы численного метода. Поэтому пользователь может тем же самым образом, с помощью контекстного меню (см. разд. 8.4), выбирать численный алгоритм приближенного решения для функции Minerr.
Пример использования функции Minerr показан в листинге 8.9. Как видно, достаточно заменить в вычислительном блоке имя функции на Minerr, чтобы вместо точного (с точностью до TOL) получить приближенное решение уравнения, заданного после ключевого слова Given.
Листинг 8.9. Приближений решение уравнения, имеющего корень (x=0, y=0)
Листинг 8.9 демонстрирует приближенное решение уравнения kx2+y2:=0, которое при любом значении коэффициента k имеет единственный точный корень (х=0,у=0). Тем не менее, при попытке решить его функцией Find для больших k, порядка принятых в листинге, происходит генерация ошибки "No solution was found" (Решение не найдено). Это связано с иным поведением функции f (x,y)=kx2+y2 вблизи ее корня, по сравнению с функциями, приводимыми в качестве примеров выше в этой главе (см. рис. 8.1, 8.2). В отличие от них, f (х,у) не пересекает плоскость f (х,у)=о, а лишь касается ее (рис. 8.7) в точке (х=0,у=0). Поэтому и найти корень изложенными в предыдущем разделе градиентными методами сложнее, поскольку вблизи корня производные f (х,у) близки к нулю, и итерации могут уводить предполагаемое решение далеко от корня.
Ситуация еще более ухудшается, если наряду с корнем типа касания (см. рис. 8.7) имеются (возможно, весьма удаленные) корни типа пересечения. Тогда попытка решить уравнение или систему уравнений с помощью функции Find может приводить к нахождению корня второго типа, даже если начальное приближение было взято очень близко к первому. Поэтому если Вы предполагаете, что система уравнений имеет корень типа касания, намного предпочтительнее использовать функцию Minerr, тем более, что всегда есть возможность проверить правильность решения уравнений простой подстановкой в них полученного решения (см. листинг 8.6).
Рис. 8.7. График функции kx2+y2
В листинге 8.9 мы рассмотрели пример нахождения существующего решения уравнения. Приведем в заключение пример нахождения функцией Minerr приближенного решения несовместной системы уравнений и неравенств (листинг 8.10). Решение, выдаваемое функцией Minerr, минимизирует невязку данной системы.
Согласно своему математическому смыслу, функция Minerr может применяться для построения регрессии серии данных по закону, заданному пользователем (см. разд. 15.2).
Листинг 8.10. Приближенное решение несовместной системы уравнений и неравенств
Как видно из листинга, в качестве результата выдаются значения переменных, наилучшим образом удовлетворяющие уравнению и неравенствам внутри вычислительного блока. Внимательный читатель может обнаружить, что решение, выдаваемое функцией Minerr в рассматриваемом примере, не является единственным, поскольку множество пар значений (х,у) в равной степени минимизирует невязку данной системы уравнений и неравенств. Поэтому для различных начальных значений будут получаться разные решения, подобно тому как разные решения выдаются функцией Find в случае бесконечного множества корней (см. обсуждение листинга 8.6 в разд. 8.3). Еще более опасен случай, когда имеются всего несколько локальных минимумов функции невязки. Тогда неудачно выбранное начальное приближение приведет к выдаче именно этого локального минимума, несмотря на то, что другой (глобальный) минимум невязки может удовлетворять системе гораздо лучше.
Задачи поиска экстремума функции означают нахождение ее максимума (наибольшего значения) или минимума (наименьшего значения) в некоторой области определения ее аргументов. Ограничения значений аргументов, задающих эту область, как и прочие дополнительные условия, должны быть определены в виде системы неравенств и (или) уравнений. В таком случае говорят о задаче на условный экстремум.
Для решения задач поиска максимума и минимума в Mathcad имеются встроенные функции Minerr, Minimize и Maximize. Все они используют те же градиентные численные методы, что и функция Find для решения уравнений. Поэтому Вы можете выбирать численный алгоритм минимизации из уже рассмотренных нами численных методов (см. разд. 8.4).
Поиск экстремума функции включает в себя задачи нахождения локального и глобального экстремума. Последние называют еще задачами оптимизации. Рассмотрим конкретный пример функции f(x), показанной графиком на рис. 8.8 на интервале (-2,5). Она имеет глобальный максимум на левой границе интервала, глобальный минимум, локальный максимум, локальный минимум и локальный максимум на правой границе интервала (в порядке слева направо).
В Mathcad с помощью встроенных функций решается только задача поиска локального экстремума. Чтобы найти глобальный максимум (или минимум), требуется либо сначала вычислить все их локальные значения и потом выбрать из них наибольший (наименьший), либо предварительно просканиро-вать с некоторым шагом рассматриваемую область, чтобы выделить из нее подобласть наибольших (наименьших) значений функции и осуществить поиск глобального экстремума, уже находясь в его окрестности. Последний путь таит в себе некоторую опасность уйти в зону другого локального экстремума, но часто может быть предпочтительнее из соображений экономии времени.
Рис. 8.8. График функции f(х)=х4+5х3-10х
Для поиска локальных экстремумов имеются две встроенные функции, которые могут применяться как в пределах вычислительного блока, так и автономно.
Minimize (f, x1, ... ,хм) — вектор значений аргументов, при которых функция f достигает минимума; Maximize (f, х1, ... ,хм) — вектор значений аргументов, при которых функция f достигает максимума; f (x1, ... , хм,...) — функция; x1, ... , xм — аргументы, по которым производится минимизация (максимизация).Всем аргументам функции f предварительно следует присвоить некоторые значения, причем для тех переменных, по которым производится минимизация, они будут восприниматься как начальные приближения. Примеры вычисления экстремума функции одной переменной (рис. 8.8) без дополнительных условий показаны в листингах 8.11- 8.12. Поскольку никаких дополнительных условий в них не вводится, поиск экстремумов выполняется для любых значений.
Листинг 8.11. Минимум функции одной переменной
Листинг 8.12. Максимум функции одной переменной
Как видно из листингов, существенное влияние на результат оказывает выбор начального приближения, в зависимости от чего в качестве ответа выдаются различные локальные экстремумы. В последнем случае численный метод вообще не справляется с задачей, поскольку начальное приближение х=-10 выбрано далеко от области локального максимума, и поиск решения уходит в сторону увеличения f (х).
В задачах на условный экстремум функции минимизации и максимизации должны быть включены в вычислительный блок, т. е. им должно предшествовать ключевое слово Given. В промежутке между Given и функцией поиска экстремума с помощью булевых операторов записываются логические выражения (неравенства, уравнения), задающие ограничения на значения аргументов минимизируемой функции. В листинге 8.13 показаны примеры поиска условного экстремума на различных интервалах, определенных неравенствами. Сравните результаты работы этого листинга с двумя предыдущими.
Листинг 8.13. Три примера поиска условного экстремума функции
Не забывайте о важности выбора правильного начального приближения и в случае задач на условный экстремум. Например, если вместо условия - 3<х<0 в последнем примере листинга задать -5<х<0, то при том же самом начальном х=-10 будет найден максимум Maximize(f,x) =-0.944, что неверно, поскольку максимальное значение достигается функцией f (х) на левой границе интервала при х=-5. Выбор начального приближения х=-4 решает задачу правильно, выдавая в качестве результата Maximize (f ,x) =-5.
Вычисление экстремума функции многих переменных не несет принципиальных особенностей по сравнению с функциями одной переменной. Поэтому ограничимся примером (листинг 8.14) нахождения максимума и минимума функции, показанной в виде графиков трехмерной поверхности и линий уровня на рис. 8.9. Привлечем внимание читателя только к тому, как с помощью неравенств, введенных логическими операторами, задается область на плоскости (х,y) .
Листинг 8.14. Экстремум функции двух переменных
Рис. 8.9. График функции f (х, у) и отрезок прямой х+у=10
Дополнительные условия могут быть заданы и равенствами. Например, определение после ключевого слова Given уравнения х+у=10 приводит к такому решению задачи на условный экстремум.
Как нетрудно сообразить, еще одно дополнительное условие привело к тому, что численный метод ищет минимальное значение функции f(x,y) вдоль отрезка прямой, показанного на рис. 8.9.
Поиск минимума можно организовать и с помощью функции Minerr. Для этого в листинге 8.14 надо поменять имя функции Minimize на Minerr, а после ключевого слова Given добавить выражение, приравнивающее функции f (х,у) значение, заведомо меньшее минимального, например f (х,у) =0.
Задачи поиска условного экстремума функции многих переменных часто встречаются в экономических расчетах для минимизации издержек, финансовых рисков, максимизации прибыли и т. п. Целый класс экономических задач оптимизации описывается системами линейных уравнений и неравенств. Они называются задачами линейного программирования. Приведем характерный пример т. н. транспортной задачи, которая решает одну из проблем оптимальной организации доставки товара потребителям с точки зрения экономии транспортных средств (листинг 8.15).
Листинг 8.15. Решение задачи линейного программирования
Модель типичной транспортной задачи следующая. Пусть имеется N предприятий-производителей, выпустивших продукцию в количестве b0,... ,bn-1, тонн. Эту продукцию требуется доставить м потребителям в количестве a0,... ,aм-1 тонн каждому. Численное определение векторов а и ь находится в первой строке листинга. Сумма всех заказов потребителей ai равна сумме произведенной продукции, т. е. всех bi (проверка равенства во второй строке). Если известна стоимость перевозки тонны продукции от 1-го производителя к j-му потребителю cij то решение задачи задает оптимальное распределение соответствующего товаропотока xij с точки зрения минимизации суммы транспортных расходов. Матрица с и минимизируемая функция f (х) матричного аргумента х находятся в середине листинга 8.15.
Условия, выражающие неотрицательность товаропотока, и равенства, задающие сумму произведенной каждым предприятием продукции и сумму заказов каждого потребителя, находятся после ключевого слова Given. Решение, присвоенное матричной переменной sol, выведено в последней строке листинга вместе с соответствующей суммой затрат. Обратите внимание, что для получения результата пришлось увеличить погрешность CTOL, задающую максимальную допустимую невязку дополнительных условий. В строке, предшествующей ключевому слову Given, определяются нулевые начальные значения ддя х простым созданием нулевого элемента матрицы xn-1, м-1
Если взять другие начальные значения для х, решение, скорее всего, будет другим! Возможно, Вы сумеете отыскать другой локальный минимум, который еще больше минимизирует транспортные затраты. Это еще раз доказывает, что задачи на глобальный минимум, к классу которых относится линейное программирование, требуют аккуратного отношения в смысле выбора начальных значений. Часто ничего другого не остается, кроме сканирования всей области начальных значений, чтобы из множества локальных минимумов выбрать наиболее глубокий.
Некоторые уравнения можно решить точно с помощью символьного процессора Mathcad. Делается это очень похоже на численное решение уравнений с применением вычислительного блока. Присваивать неизвестным начальные значения нет необходимости. Листинги 8.16 и 8.17 демонстрируют символьное решение уравнения с одним неизвестным и системы двух уравнений с двумя неизвестными соответственно.
Листинг 8.16. Символьное решение алгебраического уравнения с одним неизвестным
Листинг 8.17. Символьное решение системы алгебраических уравнений
Как видно, вместо знака равенства после функции Find в листингах следует знак символьных вычислений, который можно ввести с панели Symbolic (Символика) или, нажав клавиши <Ctrl>+<.>. He забывайте, что сами уравнения должны иметь вид логических выражений, т. е. знаки равенства нужно вводить с помощью панели Booleans (Булевы операторы). Обратите внимание, что в листинге 8.17 вычислены как два первых действительных корня, которые мы уже находили численным методом (см. разд. 8.3), так и два других мнимых корня. Эти два последних корня чисто мнимые, т. к. множитель, входящий в них.
С помощью символьного процессора решить уравнение с одним неизвестным можно и по-другому:
Введите уравнение, пользуясь панелью Booleans (Булевы операторы) или нажав клавиши <Ctrl>+<> для получения логического знака равенства, например х2+2(х-4)=0. Щелчком мыши выберите переменную, относительно которой Вы собираетесь решить уравнение. Выберите в меню Symbolics (Символика) пункт Variable / Solve (Переменная / Решить).После строки с уравнением появится строка с решением или сообщение о невозможности символьного решения этого уравнения.
В данном примере после осуществления описанных действий появляется вектор, состоящий из двух корней уравнения
Символьные вычисления могут производиться и над уравнениями, в которые, помимо неизвестных, входят различные параметры. В листинге 8.18 приведен пример решения уравнения четвертой степени с параметром а. Как видите, результат получен в аналитической форме.
Листинг 8.18. Символьное решение уравнения, зависящего от параметра
В следующем разделе мы рассмотрим более подробно, как с помощью Mathcad можно численными методами решать подобные задачи.
Решение "хороших" нелинейных уравнений и систем типа тех, которые были рассмотрены в предыдущих разделах этой главы, представляет собой несложную, с вычислительной точки зрения, задачу. В реальных инженерных и научных расчетах очень распространена более сложная проблема: решение не одного уравнения (или системы), а целой серии уравнений, зависящих от некоторого параметра (или нескольких параметров). Для таких задач существуют очень эффективные методы, которые называются методами продол-жения. Эти методы непосредственно не встроены в Mathcad, но могут быть легко запрограммированы с помощью уже рассмотренных нами средств. Будем далее говорить об одном уравнении, имея в виду, что всегда возможно обобщение результатов на случай системы уравнений.
Пусть имеется уравнение f (а,х)=0, зависящее не только от неизвестного х, но и от параметра а. Требуется определить зависимость его корня х от параметра а, т. е. х(а). Простой пример такой задачи был приведен в предыдущем разделе в листинге 8.18. Тогда нам повезло, и решение в общем виде было найдено с помощью символьных вычислений. Рассмотрим еще один, чуть более сложный пример, аналитическое решение которого также заранее известно, но с которым символьный процессор, тем не менее, не справляется. Обобщим задачу, добавив зависимость от параметра а следующим образом:
sin(ax)=0. (1)
Аналитическое решение этого уравнения находится из соотношения ах =Npi, где N=0,1,2,..., т е. имеется бесконечное количество (для каждого N) семейств решений xN=Npi/a. Несколько из семейств для N=1,2,3 показаны на рис. 8.10 тремя сплошными кривыми. Кроме того, следует иметь в виду, что для N=0, т. е. х=0, решением является прямая, совпадающая с осью х. Заметим (листинг 8.19), что символьный процессор Mathcad выдает в качестве решения только эту серию х=0.
Листинг 8.19. Символьное решение уравнения sin (а-х)=0
Забудем на время, что аналитическое решение известно, и подойдем к уравнению (1) как к любой другой новой задаче. Решим его методом секущих, применяя для этого встроенную функцию root. Самый простой, но далеко не лучший способ иллюстрируется листингом 8.20. Корень уравнения (1) требуется определить численно для каждого значения параметра а. Для этого первые две строки листинга создают ранжированную переменную 1, с помощью которой определяется вектор значений параметра а„ для которых будут производиться расчеты. Его элементы пробегают значения от о.ою до 0.025 с шагом o.ooos (эти числа взяты ради примера, Вы можете поэкспериментировать с другими значениями). Последняя строка листинга присваивает элементам еще одного вектора у вычисленные с помощью функции root значения корней уравнения (1) для каждого ах. Но, для того чтобы функция root заработала, необходимо предварительно задать начальное приближение к решению, что сделано в третьей строке листинга. Ключевой момент метода, примененного в листинге 8.20, заключается в том, что одно и то же начальное значение х=зоо использовано для всех ах.
Рис. 8.10. Решение уравнения sin(ax) =0 (см. листинг 8.21)
Листинг 8.20. Один из способов численного решения
Результат расчетов у± показан на том же рис. 8.10 точками, причем самая левая точка отвечает начальному значению у0=300. Обратите внимание, что, по мере увеличения а, кривая корней уравнения сначала идет по семейству решений у=я/а, а потом (в районе а«0.17) "перепрыгивает" на другое семейство у=2-я/а. С вычислительной точки зрения такая ситуация чаще всего крайне неблагоприятна, поскольку хотелось бы отыскать непрерывное семейство решений. Скачки зависимости у (а) могут вводить пользователя в заблуждение, вовсе скрывая от него существование семейства решений.
Почему же происходят эти скачки с одного семейства решений на другое? Конечно, причина кроется в выборе начального значения для вычисления каждого из корней. Линия начальных значений у=зоо обозначена на рис. 8.10 в виде пунктирной горизонтальной прямой. Для a0=o.oi, и вообще для нескольких первых а1 начальное значение у=зоо находится ближе всего к нижнему семейству решений у=я/а. Поэтому неудивительно, что численный метод находит именно эти корни. В правой части графика на рис. 8.10 к линии начальных значений ближе второе семейство решений у=2-я/а, к ним-то и приводит численный метод.
Приведенные соображения диктуют очень простой рецепт избавления от скачков и нахождения одного из семейств непрерывных решений. Для этого требуется при поиске каждого (i+D-го корня взять начальное значение, по возможности близкое к отыскиваемому семейству. Неплохим вариантом будет выбор приближения в виде предыдущего 1-го корня, который был найден для прошлого значения параметра аг. Возможный вариант воплощения этого метода, называемого продолжением по параметру, приведен в листинге 8.21. В нем функция root применена внутри функции пользователя f uo,a), определенной в самом начале листинга с помощью средств программирования. Назначение функции f uo,a) заключается в том, что она выдает значение корня для заданного значения параметра а и начального приближения к решению хо. В остальном, смысл листинга 8.21 повторяет предыдущий, за исключением того, что задается явно (в его предпоследней строке) только начальное значение УО=ЗОО только для поиска уг. Для всех последующих точек, как следует из последней строки листинга, взято начальное значение, равное предыдущему корню У!.
Листинг 8.21. Решение уравнения sin(ax)=0 методом продолжения по параметру
Результат вычислений, приведенный на рис. 8.11, разительно отличается от предыдущего. Как видно, столь малое изменение идеологии применения численного метода привело к определению непрерывного семейства корней. Отметим, что получить результат рис. 8.10 (без продолжения по параметру) в терминах введенной нами функции f (х0,а) можно, изменив ее первый аргумент в последней строке листинга 8.21 на константу: f oodaj.
Рис. 8.11. Решение уравнения sin(a-x)=0 методом продолжения для у0=300 (листинг 8.21)
Чтобы найти другое семейство решений, нужно взять соответствующее первое начальное значение у0, например у0=600. Результат действия листинга 8.21 для этого случая показан на рис. 8.12. Если взять у0 ближе к третьему семейству решений (например у0=300), то оно и будет найдено вычислительным процессором Mathcad, и т. д.
С помощью метода продолжения можно решать и соответствующие задачи оптимизации, зависящие от параметра. Идеология в этом случае остается точно такой же, но вместо функций решения нелинейных уравнений root или Find вам следует применить одну из функций поиска экстремума Min-err, Maximize или Minimize.
Мы привели основную идею и один из возможных способов реализации метода продолжения по параметру. Безусловно, Вы можете предложить иные как математические, так и программистские решения этой проблемы. В частности, для выбора очередного начального приближения к корню можно использовать результат экстраполяции уже найденной зависимости х(а), придумаъ более сложные алгоритмы для ветвящихся сем решений и т. д.
Рис. 8.12. Решение уравнения sin (ах) =0 методом продолжения по параметру для у0=600