Г Л А В А   8

 

ПРОГРАММИРОВАНИЕ
ЧИСЛЕННЫХ МЕТОДОВ в
Mathcad

 

 

§ 1. Элементы программирования в Mathcad

 

В настоящей главе проиллюстрированы возможности программирования, имеющиеся в пакете Mathcad, на примере реализации некоторых, из рассмотренных в настоящем пособии, алгоритмов численных методов.

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

Изначально универсальные математические пакеты, и в частности Mathcad, создавались именно для того, чтобы освободить человека, занимающегося вычислительной работой, от необходимости программирования на традиционных алгоритмических языках, таких как, Фортран, Паскаль или Си. С этой с целью в математические пакеты встраиваются наиболее популярные алгоритмы решения различных математических задач, мощные графические средства, а также разрабатывается входной язык, который, в той или иной степени, приближен к традиционному языку математической символики.

Однако практический опыт использования математических пакетов показал, что, несмотря на наличие встроенных вычислительных алгоритмов, достаточно часто возникает необходимость как в их модификации  (адаптации к специфике решаемой задачи), так и реализации собственных алгоритмов для решения специфических задач. Очевидно, что для этого математический пакет должен, во-первых, иметь средства программирования, позволяющие при необходимости реализовать любой вычислительный алгоритм, а, во-вторых, обладать открытым кодом встроенных функций, что обеспечивало бы возможность их модификации. Математический пакет Mathcad не обеспечивает доступа к алгоритмам встроенных функций, т.е. обладает закрытым кодом, однако, начиная с версии 6.0, он обрел возможность встраивать в свой документ программные модули, содержащие стандартные алгоритмические структуры: ветвление и  циклы. Более того, современные версии этого пакета допускают интеграцию своих документов с программами на традиционных алгоритмических языках, в частности Си.

 

Программные модули или программы-функции

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

Ниже приведен пример программы функции, определяющий алгоритм вычисления кусочно-заданной функции:

 


Обзор программных операторов

Количество программных операторов входного языка Mathcad невелико. Но, тем не менее, их набор достаточен для эффективной реализации любого алгоритма. Ниже приводится их перечнь:

 

AddLine               создает новую программную строку

                        оператор присваивания

if                          условный оператор

otherwise               оператор, задающий альтернативную ветвь

условного опреатора

for                         оператор цикла с параметром

while                     оператор цикла с предусловием

break                     опертор дострочного прекращения цикла или

программы

continue                оператор перехода к следующей итерации

return                    оператор, определяющий возвращаемое значение

on error                 оператор, определяющий возвращаемое значение в

случае возникновения ошибки

 

Подробный синтаксис этих операторов и примеры их использования можно найти в любом справочном пособии по пакету Mathcad, в том числе, и в интегрированной справочной системе самого пакета.

Мы же проилюстируем возможности программирвоания в пакете Mathcad на примере реализации разнообразных алгоритмов тех численных методов, которые были рассмотрены в предыдущих главах.

 

 

§ 2. Программы поиска и уточнения корней
нелинейных уравнений

 

В этом параграфе мы рассмотрим примеры программ-функций поиска и уточнения корней нелинейных уравнений с одним неизвестным. Заметим, что для решения задачи уточнений корней трансцендентных и алгебраических уравнений в пакете Mathcad существуют встроенные функции root и polyroots, реализующие методы секущих, Ридерса, Лаггера и метод сопровождающей матрицы. Мы же приведем программы для методов половинного деления и Ньютона. Начнем с метода половинного деления.

 

 

Функция bs_root (bisection root) содержит четыре формальных параметра: f – функция, нули которой ищутся; a и b – границы отрезка локализации корня; e – абсолютная погрешность. Основной алгоритм реализован с помощью циклического оператора while, в теле которого происходит вычисление средней точки отрезка , проверка знаков функции на отрезке  и выбор нового отрезка.  Для подсчета числа итераций используется переменная it. В качестве результата функция возвращает приближенное значение корня и количество итераций. Поскольку в программах-функциях пакета Mathcad в качестве результата можно указывать только одну переменную, то для получения нескольких значений, в качестве итоговой переменной может служить вектор или матрица. В данном примере формируется вектор-строка ans, первый (нулевой) элемент которой содержит приближенное значение, а второй ­­– количество итераций.

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

 

 

Функция n_root имеет те же, формальные параметры, что и функция bs_root. Поясним основные шаги алгоритма. После определения переменных it и it_max, с помощью операторов if и otherwise осуществляется выбор начального приближения к корню . В качестве начального приближения выбирается одна из границ отрезка локализации, в которой совпадают знаки функции и ее второй производной. Далее с помощью оператора while организован «бесконечный» цикл. Единица после слова while соответствует логическому значению «истина». В теле цикла происходит вычисление очередного приближения по формуле Ньютона. Цикл прекращается, если выполняется условие достижения заданной точности или же в случае превышения максимально допустимого числа итераций. Для прекращения цикла используется оператор break. Так, с помощью операторов while и break, организован аналог цикла с постусловием. Если количество итераций превысило максимально допустимое, то с помощью функции error будет выведено соответствующее сообщение. Далее происходит формирование вектора ans, содержащего приближенное значение и количество итераций.

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

Приведем пример использования рассмотренных функций для уточнения корня уравнения  с точностью . В качестве отрезка локализации выберем отрезок .

 

 

В результате вызова функций bs_root и n_root мы получили приближенное значение корня 3.693441 с шестью верными знаками после запятой и количество итераций 22 и 6 для метода половинного деления и Ньютона, соответственно. Здесь же показан пример вызова встроенной в Mathcad функции root.

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

Программа roots реализует алгоритм табличного поиска корней и содержит следующие параметры: f – функция; a и b – границы интервала поиска; h – шаг поиска; e – точность. В теле программы используются следующие переменные: n_root – номер найденного корня (нумеруются с нуля); x1 и x2 – границы текущего отрезка локализации; y1 и y2 – значения функции  на границах отрезка локализации; bs – приближенное значение корня, полученное методом половинного деления.

 

 

Поясним основные шаги алгоритма. На начальном этапе проверяется корректность исходных данных. В случае ошибки работа программы прекращается с выдачей соответствующего сообщения. Поиск корней начинается с точки x1 = a. В этой точке вычисляется значение функции y1. В теле цикла while вычисляется правая граница отрезка локализации x2 и значение функции в этой точке y2. Далее проверяются знаки функции на концах отрезка [x1, x2]. Если знаки противоположные (произведение y1 и y2 меньше нуля), то вызывается функция уточнения корня bs_root и формируется вектор результатов ans. Номер корня увеличивается на единицу. Для перехода к следующему отрезку поиска осуществляется переприсваивание x¬ x2, y¬ y1. Цикл прекращается если x2 превысит значение правой границы интервала поиска b.

Для иллюстрации работы функции roots найдем все корни уравнения  на отрезке от -10 до 10.

 

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

 

 

§ 3. Программы интерполирования зависимостей

 

В четвертой главе данного пособия были рассмотрены различные способы полиномиальной интерполяции зависимостей, в частности: глобальная интерполяция, кусочная интерполяция и интерполяция кубическими сплайнами. Встроенные в Mathcad функции позволяют осуществлять некоторые виды интерполяции, например кусочно-линейную и кубическую сплайн-интерполяцию. Мы же в качестве примеров рассмотрим программы кусочно-квадратичной интерполяции и интерполяции квадратичными сплайнами, которые не реализованы в пакете Mathcad.

Кусочно-квадратичная интерполяция. Алгоритм кусочно-квадратичной интерполяции достаточно прост и состоит из двух шагов. На первом шаге необходимо определить три узла, ближайших к точке интерполирования, а на втором шаге вычислить значение интерполяционного полинома второй степени (например, по формуле Ньютона). Способ реализация первого шага будет зависеть от того, как расположены узлы интерполирования – регулярно или нет. Для большей общности мы будем предполагать произвольное расположение узлов. В этом случае определение ближайших узлов можно выполнить циклом, в котором значение точки интерполирования x последовательной сравнивается с правыми границами отрезков интерполирования . И, в случае, если выполняется условие , то для интерполирования выбираются (k-1)-й, k-й и (k+1)-й узлы. В противном случае происходит увеличение k на единицу.

 

 

Программа-функция pinterp содержит три параметра: vx, vy – вектора, задающие исходные данные, и x – точка, в которой вычисляется значение интерполяционного полинома. Поясним алгоритм её работы. Сначала, с помощью функции last, определяется количество узлов интерполирования. Затем, переменная k, обозначающая номер локального отрезка интерполирования (отрезка ), принимает значение 1. Далее стартует цикл while, в котором  переменная k наращивается на единицу. Цикл завершается, если значение vxk окажется больше или равно значению x, или, если k станет равным (n-1). При этом для всех точек x, которые попадают в последний и предпоследний отрезки, номер k будет иметь одно и тоже значение – (n-1). Далее вычисляются коэффициенты полинома Ньютона А0, А1, А2, и, после этого, – значение полинома y.

Пример работы функции pinterp показан ниже.

 

 

 

 

 

 

 

 

 

 

 

 

 

 


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

Итак, на каждом частичном отрезке интерполирования  будем строить функцию вида:

.                                 (8.1)

Коэффициенты сплайнов  и  будем искать из следующих условий:

а) условия Лагранжа

,                                 (8.2)

б) непрерывность первой производной в узловых точках

,     .                    (8.3)

Условия (8.2) и (8.3) дают  уравнений, в то время как количество неизвестных коэффициентов – 3n. Недостающее уравнение может быть получено из дополнительных условий, накладываемых на поведение сплайна. Например, можно потребовать, чтобы значение первой производной сплайна  в точке  было бы нулевым, т.е.

.                                           (8.4)

Подстановка выражения (8.1) в (8.2), (8.3) и (8.4) приводит к следующим уравнениям

                       (8.5)

где введено обозначение .

Выразим из второго уравнения коэффициенты , предварительно подставив в него значения коэффициентов  из первого уравнения:

,                                 (8.6)

Затем,  подставив это выражение в третье уравнение системы (8.5), получим простое рекуррентное соотношение для коэффициентов

 

,   .                         (8.7)

 

Теперь алгоритм определения коэффициентов сплайнов становится достаточно очевидным. Сначала по формуле (8.7) определяем значения всех коэффициентов , с учетом того, что . Затем, по формуле (8.6), вычисляем коэффициенты . Коэффициенты  определяются из первого уравнения системы (8.5). При этом процедуру вычисления коэффициентов сплайнов нужно проделать лишь один раз.

После того, как коэффициенты рассчитаны, для вычисления самого сплайна достаточно определить номер интервала, в который попадает точка интерполирования, и воспользоваться формулой (8.1). Для определения номера интервала будем применять алгоритм, аналогичный тому, что был использован в предыдущем примере для кусочно-квадратичной интерполяции.

Рассмотрим пример программы-функции, вычисляющей коэффициенты параболических сплайнов.

 

 

В качестве входных параметров функции coeff_ps выступают векторы vx и vy, содержащие исходные данные. Алгоритм функции прост и не нуждается во дополнительном пояснении за исключением последнего шага, в котором с помощью встроенной функции augment формируется матрица коэффициентов сплайнов С, составленная из столбцов а, b и c.

Программа-функция для вычисления сплайнов может быть, например, такой:

Функция ps содержит следующие параметры: cs – матрица коэффициентов сплайнов, vx – вектор узловых точек,  x – точка, в которой вычисляется значение сплайна. Ниже приведен пример интерполяции параболическими сплайнами с использованием рассмотренных функций.

§ 4. Программы прямых и итерационных методов решения систем линейных алгебраических уравнений

 

В этом параграфе мы рассмотрим три программы, реализующие, соответственно алгоритмы метода Гаусса, метода прогонки и итерационного метода Зейделя для нахождения решений систем линейных алгебраических уравнений. Эти алгоритмы подробно рассмотрены в третьей главе пособия.

Метод Гаусса. Напомним, что данный метод содержит два этапа: прямой ход, в процессе которого матрица системы приводится к треугольному виду, и обратный ход – решение системы уравнений с треугольной матрицей. Пример программы-функции для метода Гаусса приведен ниже.

Функция gauss содержит два параметра: a – матрица системы уравнений и с – вектор свободных членов. Алгоритм функции содержит два основных цикла по переменной i. Первый цикл реализует прямой ход, а второй – обратный ход метода. В теле первого цикла предусмотрена процедура перестановки уравнений в случае, если ведущий элемент оказывается равным нулю. В качестве результата функция возвращает вектор x – вектор решения системы. Для правильной работы программы значение встроенной переменной ORIGIN должно быть равно единице.

Для иллюстрации работы функции решим следующую систему:

 

 

 

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

Ниже приведен текст программы-функции progonka, реализующий алгоритм метода прогонки. Функция содержит четыре параметра: вектор a  содержит коэффициенты при неизвестных, стоящие под главной диагональю, вектор b – на главной диагонали, вектор c – над главной диагональю. Вектор d содержит свободные члены (правая часть системы уравнений). В качестве результата функция возвращает вектор x – вектор решения системы. Для правильной работы программы значение встроенной переменной ORIGIN должно быть равно единице.

 

 

 

Проиллюстрируем работу функции progonka на примере решения следующей системы уравнений:

Итерационный метод Зейделя. Алгоритм метода Зейделя подробно описан в третьей главе. Ниже приведен текст программы zeidel, которая, содержит три параметра: a – матрица системы уравнений,  с – вектор свободных членов и e – точность вычислений.

Обратите внимание на условие окончания итераций . Здесь запись  означает вычисление евклидовой нормы вектора, т.е. расстояния между точками в n-мерном пространстве. Вычисления проводятся по формуле

.

 

Применим функцию zeidel для поиска приближенного решения системы, матрица которой является матрицей с диагональным преобладанием, что обеспечивает сходимость метода Зейделя.