CC BY
56
Cloud of Science. 2017. T. 4. № 1 http:/ / cloudofscience.ru
Сравнение встроенных методов решения систем линейных алгебраических уравнений в MathCAD
А. И. Козлов
Чувашский государственный университет им. И. Н. Ульянова 428015, Чебоксары, Московский проспект, 15
e-mail: [email protected]
Аннотация. Рассмотрены прямые методы решения систем линейных алгебраических уравнений (СЛАУ) в системах компьютерной математики MathCAD Prime и MathCAD 14. Приведены примеры использования встроенных функций c оценкой их погрешности на примере глобальной интерполяции каноническим полиномом.
Ключевые слова: MathCAD, системы линейных алгебраических уравнений, интерполяция, матрица Вандермонда, тестирование.
1. Введение
Система MathCAD обладает довольно представительным набором средств для решения СЛАУ. При этом вопросы их конкретной применимости изучены явно недостаточно. В данном случае эти средства были использованы на примере полиномиальной интерполяции. Как известно, интерполяция как вид аппроксимации применяется для получения промежуточных значений функций, заданных чаще всего в табличной форме (при условии, что табличные значения не содержат погрешностей). Канонический полином довольно широко используется в целом ряде областей вычислительной математики (численное интегрирование и дифференцирование, решение нелинейных уравнений, строительная механика и т. д.). Из известных проблем полиномов высоких степеней можно упомянуть высокую вероятность ос-циллирования [1], а также осложнения при решении систем линейных алгебраических уравнений (СЛАУ) относительно искомых коэффициентов (как правило, плохая обусловленность матрицы и в ряде случаев близость ее к вырожденной) [2]. Поэтому для собственно интерполяции лучше использовать полиномы в другой форме (Лагранжа, Чебышева и т. д.).
2. Постановка задачи
Как правило, в задачах интерполяции истинная зависимость f (х) задана в табличной форме:
f ) = Уг , * = 0 П, x0, X, •••, X — узлы интерполяции.
Требуется построить функцию P(x), совпадающую в узлах xi с функцией
f (x) P(x ) = У, i = 0, n (условия Лагранжа). Для P(x) в виде канонического полинома
Р(х) = «0 +ахх-\-----ь апх" = У^акхк
(1)
k=0
получаем СЛАУ вида:
а0+а1х0+- + айх0 = у0, а0 +ajXj + --- + cinx" = уг,
а0+а1хй+- + айхй = у„
или в матричной форме
Mv • a = y,
X х2 ■■■ х"
I Л0 Л0 о
M
v1 xn xn
(2)
nj
где Му — матрица Вандермонда.
Для существования и единственности решения (2) необходимо, чтобы узлы интерполяции были различными xi Ф х. при г Ф
МаШСЛБ предоставляет несколько способов решения (2):
- традиционный матричный;
- обращением матрицы с помощью функции geninv;
- методом Гаусса (функция гге1);
- функцией Ьо^е;
- функцией ро1усоеГ1".
Рассмотрение реализованных в них вычислительных алгоритмов является отдельной задачей, поэтому ограничимся в дальнейшем чисто практическим тестированием перечисленных способов. Стоит сразу отметить, что функция ро1усоеГГ является узкоспециализированной и предназначена только для решения (2). Остальные способы предназначены для решения произвольных СЛАУ.
Тестирование проводилось для двух диапазонов изменения аргумента: -1 — х — 1 и 0 — х —100. Для упрощения был принято равномерное расположение
узлов. В качестве интерполируемой зависимости была выбрана затухающая синусоида F(x) = x • sin(20x). Степень полинома N последовательно принимала значения 10, 50 и 100. Разумеется, трудно представить себе практическое применение полинома даже 50-й степени, однако тестирование серьезных систем компьютерной математики принято производить в «экстремальных» условиях.
MathCAD — документ в среде MathCAD Prime 3.1 (рисунок), легко позволяет получить решение всеми упомянутыми способами. Для этого достаточно передвигать в документе соответствующие функции, так чтобы нужная из них оказалась непосредственно перед вычислением максимальной невязки max |Mv • a — У .
F(ж) := е~х • sin (20 • л;) Интерполируемая зависимость
а:=— 1 Ь:=1 Л7:=10 h:= h~a k:=O..N
Л7
■■—с + к' h у —F (х)
i •= О .. N J — O..M MV. . — х * Матрица Вандермонда
График интерполируемой зависимости
[MV| = 2.3990159 • ю-11 Определитель Вандермонда
conde (MV~) = 2.03 • 104 Число обусловленности
а •■= polycoeff (л;, у) ct'=MV~l • у а — rref (augment (iVfV~, у) )
а ■= geninv (MV) - у а •■= lsolve (ikTV~. у)
Лf-=MV -а —у Невязки
г/) = 2.13-10 Максимальная невязка и индекс
L6 J соответствующего узла
лг
■Р(ж) != У]а -хк Канонический полином
Процедура нахождения максимальной невязки и индекса соответствующего узла
Рисунок. Текст MathCAD-документа
3. Результаты
Таблица 1. Результаты тестирования при -1 < х. < 1 (в числителе — максимальная невязка, в знаменателе — индекс соответствующего ей узла)
Степень полинома N Определитель Вандермонда Число обусловленности Метод /Функция
Матричный geninv lsolve rref polycoeff
10 2.39-10-11 2ЯМ04 1.29-10-12 0 0.25 5 2.1310-13 6 7.4-10-13 10 2.69-10-13 0
50 я 0 4.48-1018 770.42 0 1.45 19 2.9-10-6 10 3.48 0 2.96-10-9 0
100 я 0 я да 400 0 1.45 38 5.64-10-8 91 2.36 7 4.68-1016 1
Таблица 2. Результаты тестирования при 0 — x. <
Степень полинома N Определитель Вандермонда Число обусловленности Метод/Функция
Матричный geninv lsolve rref polycoeff
10 6.66-1082 я да 1.46-10-13 10 2.09- 10-4 10 1.49^ 10-15 6 4.6-10-14 10 8.8-10-16 10
50 я да я да 1.73^ 108 50 107.47 50 0.05 35 8.73-108 50 1.3510125 50
100 я да я да 1.37^ 1012 100 0.6 95 1.211012 100
Результаты, приведенные в табл. 1, 2, показывают, что матрица Вандермонда для обоих диапазонов изменения аргумента является плохо обусловленной, а при —1 — x — 1 дополнительно к этому еще и близка к вырожденной (определитель стремится к нулю). Тем не менее при N = 10 все использованные способы (кроме функции geninv) выдали неплохие и сопоставимые результаты (табл. 1). При N = 50 начали проигрывать матричный способ и функция rref (реализующая метод Гаусса). Функция polycoeff при N = 100 неожиданно показала себя хуже всех. Очевидным лидером во всех случаях показала себя функция lsolve.
При 0 — x —100 (табл. 2) все способы, включая geninv, опять неплохо проявили себя при N = 10. А вот при N = 50 все, кроме lsolve, стремительно пошли «вразнос». При N = 100 geninv и polycoeff вообще отказались выдать результат.
Дополнительно следует отметить, что не во всех случаях максимальные невязки отмечаются в крайних точках отрезка интерполяции.
По утверждениям компании PTC программный код MathCAD Prime, включая вычислительные алгоритмы, написан с «нуля». Неоднозначность приведенных результатов побудила сравнить их с аналогичными результатами в среде предыдущих версий Mathsoft MathCAD. Результаты расчетов в среде Mathcad 14 приведены в табл. 3, 4.
Таблица 3. Результаты тестирования при -1 < xi < 1
Степень полинома N Определитель Вандермонда Число обусловленности Метод/Функция
Матричный geninv lsolve rref polycoeff
10 2.69-10-13 2.03-104 8.38-10"13 0 0.25 5 1.6510-13 9 7.40-10"13 10 2.69-10"13 0
50 я 0 4.48-1018 121.13 50 1.45 19 9.1910-8 42 3.48 0 2.96-10"9 0
100 Я 0 — 3.11103 100 1.45 38 4.29-10"7 3 2.36 7 2.88-1017 100
Таблица 4. Результаты тестирования при 0 < xi < 100
Степень полинома N Определитель Вандермонда Число обусловленности Метод/Функция
Матричный geninv lsolve rref polycoeff
10 6.66-1082 — 9.659-10"14 10 2.09-10-4 10 5.33-10"15 9 4.42-10-14 10 2.66-10-15 10
50 0 — 3.84-109 50 108.47 50 0.21 49 8.73-108 50 1.14-10125 50
100 0 — 1.011011 100 — 1.46 96 1.211012 100 —
Практическое совпадение в обеих версиях MathCAD показали только geninv и rref, что, скорее всего, говорит об одинаковости реализованных алгоритмов. Несколько ухудшились результаты стандартного матричного метода и в ряде случаев функции lsolve. Функция polycoeff ведет себя или одинаково, или несколько хуже. В этой связи хотелось бы порекомендовать компании PTC обратить внимание на приведенный в [3] эффективный алгоритм обращения матрицы Вандермонда.
Попутно следует отметить непонятные результаты вычисления в MathCAD 14 определителя Вандермонда при N = 50 и N = 100 (табл. 4). Вместо больших значений или сообщения о невозможности его вычисления получились строгие нули. Отмеченная особенность вряд ли имеет практическую ценность в связи с переходом (во многом, вынужденным) многих пользователей прежних версий Mathsoft MathCAD на MathCAD Prime.
В целом следует констатировать избыточность способов решения СЛАУ, заложенных в среде MathCAD, независимо от версий. В данном примере лучше всего себя проявила функция lsolve. Данный вывод трудно считать однозначным, поскольку сделан он на основе интерполяции всего лишь одной зависимости. Вполне возможно, что в иных случаях лидерство может и поменяться. Очевидно, при решении любой произвольной СЛАУ стоит перебором упомянутых способов выбрать тот, который дает минимальные невязки.
Литература
[1] Очков В. Ф., Богомолова Е. П. Интерполяция, экстраполяция, аппроксимация или «Ложь, наглая ложь и статистика» // Cloud of Science. 2015. T. 2. № 1. С. 61-88.
[2] Интерполяция функций интерполяционными полиномами [Электронный ресурс] URL: http://matlab. exponenta.ru/spline/book1/10 .php
[3] Кольвах В. Ф., Кольвах Д. В. Эффективный алгоритм обращения матрицы Вандермон-
да // Труды молодых ученых. Владикавказский научный центр РАН. 2003. № 4. С. 17-20.
Автор:
Александр Иванович Козлов — кандидат технических наук, доцент, доцент кафедры «Электроснабжение промышленных предприятий», Чувашский государственный университет им. И. Н. Ульянова
Comparison of the built-in methods of the decision systems of the linear algebraic equations in MathCAD
A. I. Kozlov
Chuvash state university n.a. I. N. Uliyanova Moscow avenue, 15, Cheboksary, Russia 428015
e-mail: [email protected]
Abstract. They are considered direct methods of the decision of the systems of the linear algebraic equations (SLAE) in systems computer mathematicians MathCAD Prime and MathCAD 14. Cite an instance use built-in function with estimation to their inaccuracy on example of the global interpolation by canonical multinomial.
Key words: MathCAD, systems of the linear algebraic equations, interpolation, matrix of Vandermond, testing.
References
[1] Ochkov V. F., Bogomolova E. P. (2015) Cloud of Science, 2(1):61-88. [In Rus]
[2] http://matlab.exponenta. ru/spline/book1/10 .php
[3] Kolivah V. F., Kolivah D. V. (2003) Works young scientist. Vladikavkazskiy scientific centre by RAS, 4:17-20. [In Rus]
