ИКТ

5.1. Интегрирование

Интегрирование в Mathcad реализовано в виде вычислительного оператора. Допускается вычислять интегралы от скалярных функций в пределах интегрирования, которые также должны быть скалярами. Несмотря на то что пределы интегрирования обязаны быть действительными, подынтегральная функция может иметь и комплексные значения, поэтому и значение интеграла может быть комплексным. Если пределы интегрирования имеют размерность (см. разд. "Размерные переменные" гл 4), то она должна быть одной и той же для обоих пределов.

5.1.1. Операторы интегрирования Интегрирование, дифференцирование, как и множество других математических действий, устроено в Mathcad по принципу "как пишется, так и вводится". Чтобы вычислить определенный интеграл, следует напечатать его обычную математическую форму в документе. Делается это с помощью панели Calculus (Вычисления) нажатием кнопки со значком интеграла или вводом с клавиатуры сочетания клавиш <Shift>+<7> (или символа "&"). Появится символ интеграла с несколькими местозаполнителями (рис. 7.1), в которые нужно ввести нижний и верхний интервалы интегрирования, подынтегральную функцию и переменную интегрирования. Можно вычислять интегралы с одним или обоими бесконечными пределами. Для этого на месте соответствующего предела введите символ бесконечности, воспользовавшись, например, той же самой панелью Calculus (Вычисления). Чтобы ввести -«> (минус бесконечность), добавьте знак минус к символу бесконечности, как к обычному числу. Рис. 7.1. Оператор интегрирования Чтобы получить результат интегрирования, следует ввести знак равенства или символьного равенства. В первом случае интегрирование будет проведено численным методом, во втором — в случае успеха, будет найдено точное значение интеграла с помощью символьного процессора Mathcad. Эти два способа иллюстрирует листинг 7.1. Конечно, символьное интегрирование возможно только для небольшого круга несложных подынтегральных функций. Листинг 7.1. Численное и символьное вычисление определенного интеграла Подынтегральная функция может зависеть от любого количества переменных. Именно для того чтобы указать, по какой переменной Mathcad следует вычислять интеграл, и нужно вводить ее имя в соответствующий местозаполнитель. Помните, что для численного интегрирования по одной из переменных предварительно следует задать значение остальных переменных, от которых зависит подынтегральная функция и для которых вы намерены вычислить интеграл (листинг 7.2). Листинг 7.2. Интегрирование функции двух переменных по разным переменным Оператор интегрирования может использоваться точно так же, как и другие операторы: для определения функций, в циклах и при вычислении ранжированных переменных. Пример присваивания пользовательской функции д(х) значения определенного интеграла и вычисления нескольких ее значений приведен в листинге 7.3. Листинг 7.3. Использование оператора интегрирования в функции пользователя

5.1.2. Об алгоритмах интегрирования

Результат численного интегрирования — это не точное, а приближенное значение интеграла, определенное с погрешностью, которая зависит от встроенной константы TOL. Чем она меньше, тем с лучшей точностью будет найден интеграл, но и тем больше времени будет затрачено на расчеты. По умолчанию TOL=0.001. Для того чтобы ускорить вычисления, можно установить меньшее значение TOL.

Если скорость расчетов имеет для Вас принципиальное значение, например при многократном вычислении интеграла внутри цикла, проявите осторожность, выбирая значение точности Обязательно поэкспериментируйте на тестовом примере с характерной для Ваших расчетов подынтегральной функцией. Посмотрите, как уменьшение константы TOL сказывается на погрешности интегрирования, вычислив интеграл для разных ее значений и выбрав оптимальное, исходя из соотношения точность / скорость вычислений.

Отдавайте себе отчет в том, что при вводе в редакторе Mathcad оператора численного интегрирования, Вы, фактически, создаете самую настоящую программу. Например, программой является первая строка листинга 7.1, просто большая часть ее скрыта от Вашего взора разработчиками компании MathSoft. В большинстве случаев об этом не приходится специально задумываться, можно полностью положиться на Mathcad. Но иногда может потребоваться умение управлять параметрами этой программы, как мы уже рассмотрели на примере выбора константы TOL. Кроме нее, пользователь имеет возможность выбирать сам алгоритм численного интегрирования. Для этого:

• Щелкните правой кнопкой мыши в любом месте на левой части вычисляемого интеграла.
• В появившемся контекстном меню выберите один из четырех численных алгоритмов (рис. 7.2).

Рис. 7.2. Выбор алгоритма численного интегрирования

Обратите внимание, что, перед тем как один из алгоритмов выбран впервые, как показано на рис. 7.2, флажок проверки в контекстном меню установлен возле пункта AutoSelect (Автоматический выбор). Это означает, что алгоритм определяется Mathcad, исходя из анализа пределов интегрирования и особенностей подынтегральной функции. Как только один из алгоритмов выбран, этот флажок сбрасывается, а избранный алгоритм отмечается точкой.

Разработчиками Mathcad 11 запрограммированы четыре численных метода интегрирования:

• Romberg (Ромберга) — дли большинства функций, не содержащих особенностей;
• Adaptive (Адаптивный) — для функций, быстро меняющихся на интервале интегрирования;
• Infinite Limit (Бесконечный предел) — для интегралов с бесконечными пределами ();
• Singular Endpoint (Сингулярная граница) — для интегралов с сингулярностью на конце. Модифицированный алгоритм Ромберга для функций, не определенных на одном или обоих концах интервала интегрирования.

Старайтесь все-таки оставить выбор численного метода за Mathcad, установив флажок AutoSelect (Автоматический выбор) в контекстном меню. Попробовать другой метод можно, например, чтобы сравнить результаты расчетов в специфических случаях, когда у Вас закрадываются сомнения в их правильности.

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

• Сначала строится несколько интерполирующих полиномов, которые заменяют на интервале интегрирования подынтегральную функцию f(x). В качестве первой итерации полиномы вычисляются по 1,2 и 4 интервалам. Например, первый полином, построенный по 1 интервалу, — это просто прямая линия, проведенная через две граничные точки интервала интегрирования, второй — квадратичная парабола и т. д.
• Интеграл от каждого полинома с известными коэффициентами легко вычисляется аналитически. Таким образом, определяется последовательность интегралов от интерполирующих полиномов: ilf i2, i4,... Например, по правилу трапеций i1=(b-a) (f (a)+f (b) )/2 и т. д.
• Из-за интерполяции по разному числу точек вычисленные интегралы ii, i2,... несколько отличаются друг от друга. Причем, чем больше точек используется для интерполяции, тем интеграл от интерполяционного полинома ближе к искомому интегралу, стремясь к нему в пределе бесконечного числа точек. Поэтому определенным образом осуществляется экстраполяция последовательности ilf i2, I4,... до нулевой ширины элементарного интервала. Результат этой экстраполяции J принимается за приближение к вычисляемому интегралу.
• Осуществляется переход к новой итерации с помощью еще более частого разбиения интервала интегрирования, добавления нового члена последовательности интерполирующих полиномов и вычисления нового (N-ГО) приближения Ромберга JN.
• Чем больше количество точек интерполяции, тем ближе очередное приближение Ромберга к вычисляемому интегралу и, соответственно, тем меньше оно отличается от приближения предыдущей итерации. Как только разница между двумя последними итерациями | JN-j11"11 становится меньше погрешности TOL или меньше TOL-|JN|, итерации прерываются, и JN появляется на экране в качестве результата интегрирования.

Об алгоритме полиномиальной сплайн-интерполяции см .гл. 15.

Если интеграл расходится (равен бесконечности), то вычислительный процессор Mathcad может выдать сообщение об ошибке, выделив при этом оператор интегрирования, как обычно, красным цветом. Чаще всего ошибка будет иметь тип "Found a number with a magnitude greater than 10^307" (Найдено число, превышающее значение 10³⁰⁷) или "Can't converge to a solution" (Не сходится к решению), как, например, при попытке вычислить интеграл. Тем не менее, символьный процессор справляется с этим интегралом, совершенно правильно находя его бесконечное значение (листинг 7.4).

Листинг 7.4. Символьное вычисление расходящегося интеграла

Символьный процессор предоставляет замечательные возможности аналитического вычисления интегралов, в том числе зависящих от параметров и неопределенных интегралов, как показано в листингах 7.5 и 7.6. Об этом и о вычислении интегралов с помощью меню Symbolics (Символика), упоминалось в гл. 5.

 

Листинг 7.5. Символьное вычисление интеграла с переменным пределом

Листинг 7.6. Символьное вычисление неопределенного интеграла

При попытке численного решения задачи из листинга 7.4 методом, отличным от алгоритма вычисления интегралов с бесконечными пределами (Infinite Limit), будет получено неверное решение (листинг 7.7) — вместо бесконечности выдано большое, но конечное число, немного не дотягивающее до численной бесконечности, являющейся для вычислительного процессора просто большим числом 10³⁰⁷ (см. разд. "Встроенные константы" гл. 4). Отметим, что Mathcad в режиме автоматического выбора алгоритма (AutoSelect) предлагает именно алгоритм Infinite Limit.

Листинг 7.7. Плохо выбранный численный алгоритм неверно находит расходящийся интеграл

5.1.4. Кратные интегралы

Для того чтобы вычислить кратный интеграл:

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

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

Листинг 7.8. Символьное и численное вычисление кратного интеграла

Аккуратнее вводите в редакторе Mathcad кратные интегралы, если они имеют различные пределы интегрирования по разным переменным. Не перепутайте пределы, относящиеся к разным переменным. Если Вы имеете дело с такого рода задачами, обязательно разберитесь с листингом 7.9, в котором символьный процессор вычисляет двукратный интеграл. В первой строке пределы интегрирования [а, b] относятся к переменной у, а во второй строке — к переменной X.

Листинг 7.9. Символьное вычисление кратных интегралов

5.2. Дифференцирование С помощью Mathcad можно вычислять производные скалярных функций любого количества аргументов, от о-го до 5-го порядка включительно. И функции, и аргументы могут быть как действительными, так и комплексными. Невозможно дифференцирование функций только вблизи точек их сингулярности. Вычислительный процессор Mathcad обеспечивает превосходную точность численного дифференцирования. Но больше всего пользователь оценит возможности символьного процессора, который позволяет с легкостью осуществить рутинную работу вычисления производных громоздких функций, поскольку, в отличие от всех других операций, символьное дифференцирование выполняется успешно для подавляющего большинства аналитически заданных функций. В Mathcad 11 для ускорения и повышения точности численного дифференцирования функций, заданных аналитически, автоматически задействуется символьный процессор (см. разд. "Оптимизация вычислений"гл. 3).

5.2.1. Первая производная

Для того чтобы продифференцировать функцию f (х) в некоторой точке:

• Определите точку х, в которой будет вычислена производная, например х:=1.
• Введите оператор дифференцирования нажатием кнопки Derivative (Производная) на панели Calculus (Вычисления) или введите с клавиатуры вопросительный знак <?>.
• В появившихся местозаполнителях (рис. 7.3) введите функцию, зависящую от аргумента х, т. е. f(х), и имя самого аргумента х.
• Введите оператор <=> численного или < -> > символьного вывода для получения ответа.

Рис. 7.3. Оператор дифференцирования

Пример дифференцирования функции f(x)=cos(x)ln(x) приведен в листинге 7.10.

Листинг 7.10. Численное и символьное дифференцирование

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

Рис. 7.4. Ошибка в применении оператора дифференцирования

Конечно, можно, как и при использовании других операторов, предварительно определить функцию в отдельном выражении, а затем посчитать ее производную (см. листинг 7.11); или применить оператор дифференцирования для определения собственных функций пользователя (см. листинг 7.12).

 

Листинг 7.11. Символьное и численное дифференцирование функции пользователя

Листинг 7.12. Определение функции через оператора дифференцирования

В обоих листингах первой строкой определяется функция f (x)=1/x. Во второй строке листинга 7.11 с помощью символьного процессора находится аналитическое выражение ее производной, а в оставшейся части, подобно листингу 7.10, сначала численно, а затем аналитически определяются значения этой производной в точке х=0.1. В листинге 7.12 через производную от f (х) определяется еще одна пользовательская функция д(х) и затем находится ее конкретное значение в той же точке х=0.1.

Как Вы заметили, оператор дифференцирования, в основном, соответствует его общепринятому математическому обозначению. Однако в некоторых случаях при его вводе следует проявить осторожность. Рассмотрим один показательный пример, приведенный в листинге 7.13. Его первые две строки вычисляют производную sin(x) в точке х=0.5. Последняя строка демонстрирует неправильное применение оператора дифференцирования. Вместо вычисления производной sin(x) в той же точке, как этого можно было ожидать, получено нулевое значение. Это случилось потому, что аргумент функции sin(x) введен не в виде переменной х, а в виде числа. Поэтому Mathcad воспринимает последнюю строку как вычисление сначала значения синуса в точке х=0.5, а затем дифференцирование этого значения (т. е. константы) также в точке х=0.5, в соответствии с требованием первой строки листинга. Поэтому ответ, на самом деле, неудивителен — в какой точке ни дифференцируй константу, результатом будет ноль.

Листинг 7.1З. Приер правильного и неправильного применения дифференцирования

Для численного дифференцирования Mathcad применяет довольно сложный алгоритм, вычисляющий производную с колоссальной точностью до 7-8-го знака после запятой. Этот алгоритм (метод Риддера) описан во встроенной справочной системе Mathcad, доступной через меню Help (Справка). Погрешность дифференцирования не зависит от констант TOL или CTOL, в противоположность большинству остальных численных методов, а определяется непосредственно алгоритмом.

Исключение составляют функции, которые дифференцируются в окрестности сингулярной точки; например для рассмотренной нами функции f (x)=1/x это будут точки вблизи х=о. При попытке найти ее производную при х=о будет выдано сообщение об одной из ошибок деления на ноль "Can't divide by zero" (Деление на ноль невозможно) или "Found a singularity while evaluating this expression. You may be dividing by zero" (Найдена сингулярность при вычислении этого выражения. Возможно, Вы делите на ноль). Если попробовать численно определить производную очень близко к нулю, например, при х=10^-100, то может появиться сообщение об ошибке "Can't converge to a solution" (Невозможно найти решение). Встретившись с одной из упомянутых ошибок, присмотритесь повнимательнее к дифференцируемой функции и убедитесь, что Вы не имеете дело с точкой сингулярности.

5.2.2. Производные высших порядков

Mathcad позволяет численно определять производные высших порядков, от 0-го до 5-го включительно. Чтобы вычислить производную функции f (х) N-го порядка в точке х, нужно проделать те же самые действия, что и при взятии первой производной (см. разд. 7.2.1), за тем исключением, что вместо оператора производной необходимо применить оператор N-й производной (Nth Derivative). Этот оператор вводится с той же панели Calculus (Вычисления) либо с клавиатуры нажатием клавиш <Ctrl>+<?> и содержит еще два местозаполнителя, в которые следует поместить число N. В полном соответствии с математическим смыслом оператора, определение порядка производной в одном из местозаполнителей приводит к автоматическому появлению того же числа в другом из них.

"Производная" при N=0 по определению равна самой функции, при N=1 получается обычная первая производная. Листинг 7.14 демонстрирует численное и символьное вычисление второй производной. Обратите внимание, что, как и при вычислении обычной производной, необходимо перед оператором дифференцирования присвоить аргументу функции значение, для которого будет вычисляться производная.

Листинг 7.14. Численное и символьное вычисление второй производной

Убедиться в том, что символьный процессор Mathcad в последней строке листинга 7.14 дает тот же результат, что и вычислительный процессор в предыдущей строке, можно, упростив его. Для этого следует выделить полученное последнее выражение и выбрать в меню Symbolics (Символика) пункт Symplify (Упростить). После этого ниже появится еще одна строка с численным результатом выделенного выражения.

Чтобы вычислить производную порядка выше 5-го, следует последовательно применить несколько раз оператор N-й производной, подобно тому как вводились операторы кратного интегрирования (см. разд. 7.1.4). Однако для символьных вычислений этого не потребуется — символьный процессор умеет считать производные порядка выше 5-го. Сказанное иллюстрирует листинг 7.15, в котором сначала численно, а затем символьно вычисляется седьмая производная синуса в точке х=0.1.

Листинг 7.15. Численное и символьное вычисление седьмой производной

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

Из сказанного ясно, что падение точности при численном расчете высших производных может быть очень существенно. В частности, если попытаться определить девятую производную синуса, подобно идее листинга 7.15, то в качестве результата будет выдан нуль, в то время как истинное значение девятой производной равно cos (0.1).

7.2.3. Частные производные

С помощью обоих процессоров Mathcad можно вычислять производные функций любого количества аргументов. В этом случае, как известно, производные по разным аргументам называются частными. Чтобы вычислить частную производную, необходимо, как обычно, ввести оператор производной с панели Calculus (Вычисления) и в соответствующем местозаполнителе напечатать имя переменной, по которой должно быть осуществлено дифференцирование. Пример приведен в листинге 7.16, в первой строке которого определена функция двух переменных, а в двух следующих строках символьным образом вычислены ее частные производные по обеим переменным — х и у — соответственно. Чтобы определить частную производную численным методом, необходимо предварительно задать значения всех аргументов, что и сделано в следующих двух строках листинга. Последнее выражение в листинге снова (как и в третьей строке) определяет символьно частную производную по у. Но поскольку переменным х и у уже присвоено конкретное значение, то в результате получается число, а не аналитическое выражение.

 

Листинг 7.16 Символьное и численное вычисление частных производных

Частные производные высших порядков рассчитываются точно так же, как и обычные производные высших порядков (см. разд. 7.2.2). Листинг 7.17 иллюстрирует расчет вторых производных функции из предыдущего примера по переменным х, у и смешанной производной.

Листинг 7.17. Вычисление второй частной производной

Возможно, Вы обратили внимание, что в обоих листингах 7.16 и 7.17 оператор дифференцирования записан в форме частной производной. Подобно тому как существует возможность выбирать вид, например оператора присваивания, можно записывать операторы дифференцирования в виде обычной или частной производной. Запись оператора не влияет на вычисления, а служит лишь более привычной формой представления расчетов. Для того чтобы изменить вид оператора дифференцирования на представление частной производной, следует:

• Вызвать контекстное меню из области оператора дифференцирования нажатием правой кнопки мыши.
• Выбрать в контекстном меню верхний пункт View Derivative As (Показывать производную как).
• В появившемся подменю (рис. 7.5) выбрать пункт Partial Derivative (Частная производная).

Рис. 7.5. Изменение вида оператора дифференцирования

Чтобы вернуть вид производной, принятый по умолчанию, выберите в подменю пункт Default (По умолчанию) либо, для представления ее в обычном виде — Derivative (Производная).

Завершим разговор о частных производных двумя примерами, которые нередко встречаются в вычислительной практике. Программная реализация первого из них, посвященная вычислению градиента функции двух переменных, приведена в листинге 7.18. В качестве примера взята функция f (x,y), определяемая в первой строке листинга, график которой показан в виде линий уровня на рис. 7.6. Как известно, градиент функции f(х,у) является векторной функцией тех же аргументов, что и f (х,у), определенной через ее частные производные, согласно второй строке листинга 7.18. В оставшейся части этого листинга задаются ранжированные переменные и матрицы, необходимые для подготовки графиков функции и ее градиента.

Листинг 7.18. Вычисление градиента

Векторное поле рассчитанного градиента функции f(x,y) показано на рис. 7.7. Как можно убедиться, сравнив рис. 7.6 и 7.7, математический смысл градиента состоит в задании в каждой точке (х,у) направления на плоскости, в котором функция f (х,у) растет наиболее быстро.

До сих пор в данной главе мы рассматривали скалярные функции, к которым, собственно, и можно применять операторы дифференцирования. Часто приходится иметь дело с вычислением производных векторных функций. Например, в различных областях математики (см. разд. "Жесткие системы ОДУ" гл. 11) мы сталкиваемся с проблемой вычисления якобиана (или матрицы Якоби) — матрицы, составленной из частных производных векторной функции по всем ее аргументам. Пример вычисления якобиана векторной функции f (х) векторного аргумента х приведен в листинге 7.19. В нем для определения частных производных якобиана каждый i-й скалярный компонент f (x)_i дифференцируется символьным процессором Mathcad.

Рис. 7.6. График линий уровня функции f (х,у) (листинг 7.18)

Рис. 7.7. Векторное поле градиента функции f (х,у) (листинг 7.18)

Тот же самый якобиан можно вычислить и несколько по-другому, если определить функцию не одного векторного, а трех скалярных аргументов f(x,y,z) (листинг 7.20).

Листинг 7.19. Вычисление якобиана векторной фенкции векторного аргумента

Листинг 7.20. Вычисление якобиана векторной функции трех скалярных аргумента

Не забывайте, что для численного определения якобиана необходимо сначала определить точку, в которой он будет рассчитываться, т. е. вектор х в терминах листинга 7.19, или все три переменных х,у, z в обозначениях листинга 7.20.