Линейное уравнение на языке с

Решение линейного уравнения

Решение уравнения на Си.Программа на си должна решать линейное уравнение

Задается линейное уравнение ax=b. Задаются параметры уравнения a и b. Программа определяет, есть ли решение уранвения, если решение есть,то оно выводится на экран.При написании программы полезными будет уроки Ввод данных с клавиатуры в Си и Логические операции в Си Описание алгоритма решения линейного уравнения на Си. В начале программы с клавиатуры задаются значения параметров уравнения a и b.
В уравнении ax=b необходимо проанализировать все возможные варианты.
Если a=0 и b=0, то решение уравнения это любое действительное число
Если a=0 b и b≠0,т уравнения не имеет решений
Если a≠0,то решение уравнения рассчитывается по формуле x=b/a
Чтобы проверить выполнение двух условий одновременно, необходимо использовать логическую функцию И. В Си он задается с помощью &&.

Возможный вариант реализации программа на си решение линейного уравнения
#include
#include
#include

На заметку.
1.Чтобы вывести вещественное число на экран без лишних нулей используется %g
printf («Ищем решение уравнения %gx=%g «, a,b);
2.Чтобы задать условие равенства 0 используется ==. Например условие a=0 задается в Си a==0. Чтобы задать условие не равенства 0 используется логическое отрицание, которое задается знаком ! Например b≠0 задается в Си !(b==0)
Задания для самостоятельной работы:
1.Решение квадратного уравнения. Задается уравнение ax^2+bx+c=0. Задаются параметры a и b и c с клавиатуры. Определяется дискриминант. Если решения есть , то выводятся на экран, если решения нет то сообщается об этом.
В этой программе необходимо вычислить дискриминант d=b*b-4*a*c
Если d 0, то два решения x1=(-b+sqrt(d))/2a и x2 =(-b-sqrt(d))/2a

Возможный вариант реализации программа на си решение квадратного уравнения

#include
#include
#include
#include

main()
<
float a; // коэффициент при x^2
float b; // коэффицент при x
float c;// свободный член
float d;// переменная для дискрименанта

float x1;// переменная для первого корня
float x2;// перменная для второго корня

// ввод коэффициентов уравнения

printf («Введите коэффициент a\n»);
scanf(«%f», &a);
printf («Введите коэффициент b\n»);
scanf(«%f», &b);
printf («Введите свободный член с\n»);
scanf(«%f», &c);
printf («Ищем решение уравнения %gx^2+%gx+%g=0\n», a,b,c);
// расчет дискрименанта
d=b*b-4*a*c;

// анализ дискрименанта
if (d

if (d>0)
<
x1=(-b+sqrt(d))/(2*a);
x2=(-b-sqrt(d))/(2*a);
printf («Уравнение имеет два корня\n»);
printf(«x1=%g\n»,x1);
printf(«x2=%g\n»,x2);

Полезно почитать по теме условия в си

Курсовая работа (Д.664.2.1.12.02.002.09.ПЗ) по дисциплине «Программирование на языке высокого уровня С/С++» Тема: «Решение системы линейных уравнений методом LUP-разложения»

ФЕДЕРАЛЬНОЕ ГОСУДАРСТВЕННОЕ БЮДЖЕТНОЕ ОБРАЗОВАТЕЛЬНОЕ УЧРЕЖДЕНИЕ ВЫСШЕГО ПРОФЕССИОНАЛЬНОГО ОБРАЗОВАНИЯ

«ВОСТОЧНО-СИБИРСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ ТЕХНОЛОГИЙ И УПРАВЛЕНИЯ»

Кафедра систем информатики

по дисциплине «Программирование на языке высокого уровня С/С++» Тема: «Решение системы линейных уравнений методом

студент гр. д664

доцент кафедры СИ

Тема курсовой работы посвящена разработке алгоритма решение системы линейных уравнений с неизвестными вида методом LUP-разложения.

В ходе работы было определено понятие LUP — разложения, описан метод решения задачи, входные и выходные данные, и был составлен эффективный алгоритм, записанный в виде блок-схемы.

1 Теоретический материал 5

1.1 Словесная постановка задачи 5

1.2 Анализ предметной области 5

2 Проектный раздел 6

2.1 Формальная постановка задачи 6

2.2 Алгоритм решения задачи 8

3 Экспериментальный раздел 12

Список использованных источников 15

Введение

Системы линейных уравнений возникают в самых различных задачах. Такую систему можно записать в виде матричного уравнения . Если матрица не вырождена, то для нее найдется обратная матрица , и вектор является решением. Однако на практике этот подход очень часто сталкивается с вычислительной неустойчивостью (numerical instability): вещественные числа хранятся в памяти лишь приближённо и ошибки округления могут накапливаться. Другой подход – так называемое LUP — разложение – в меньшей степени подвержен этой опасности, к тому же, работает он в три раза быстрее.

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

Для достижения целей работы был поставлены и решены следующие задач:

изучить предметную область задачи; описать метод решения задачи; реализовать алгоритм в виде блок-схемы; провести тестирование алгоритма; сделать вывод по проделанной работе.

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

Теоретический материал Словесная постановка задачи

Разработать эффективный алгоритм решения системы линейных уравнений с неизвестными вида методом LUP-разложения.

Анализ предметной области

LUP-разложение (LUP-декомпозиция) – представление данной матрицы в виде произведения , где матрица является нижне-треугольной с единицами на главной диагонали, – верхне-треугольная общего вида, а – матрица перестановок, получаемая из единичной матрицы путём перестановки строк или столбцов. Такое разложение можно осуществить для любой невырожденной матрицы. LUP-разложение используется для вычисления обратной матрицы по компактной схеме, вычисления решения системы линейных уравнений. По сравнению с алгоритмом LU-разложения алгоритм LUP-разложения может обрабатывать любые невырожденные матрицы и при этом обладает более высокой устойчивостью[1].

Зачем нужно LUP-разложение? Зная его, мы сводим систему к двум системам с треугольными матрицами, а такие системы решаются просто. Вот как происходит сведение. Умножим обе части уравнения системы на , получаем уравнение, а согласно произведению , мы имеем . Выразим , и нам останется решить две системы: , где мы найдем , и следом , откуда мы и найдем конечный ответ.

Проектный раздел Формальная постановка задачи

Входные данные: – количество линейных уравнений и количество неизвестных в них;

– матрица коэффициентов при неизвестных для всех линейных уравнений;

– столбец свободных членов;

Выходные данные: — неизвестные, которые надо определить, являющиеся решением СЛАУ.

Для начала необходимо построить LUP – разложение невырожденной матрицы .

Опишем рекурсивный метод построения LUP – разложения. При (случай при – тривиален) разобьем матрицу на четыре части:

,

где – вектор столбец из элементов , – вектор-строка из элементов , – матрица размера .

Матрица состоящая из разности и называется дополнением Шура (Schur complement) элемента в матрице . Дополнение Шура – матрица размером , и мы можем найти для него LUP – разложение, таким образом, пройдясь по всей матрице рекурсивно.

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

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

После построения LUP – разложения нам остается решить две системы с треугольными матрицами. Сначала из нижнетреугольной системы мы найдем методом прямой подстановки, а затем из верхнетреугольной системы найдем методом обратной подстановки. Вектор будет искомым решением [2].

Прямая подстановка (forward substitution) позволяет решить нижне-треугольную систему за время , пройдясь от первого уравнения к последнему. Матрица является нижне-треугольной с единицами на главной диагонали, поэтому система имеет вид:

Первое уравнение сразу дает нам значение . Подставив его во второе уравнение, получим . Теперь, подставив и в третье уравнение, получим . И так далее: для получения подставим в — е уравнения и получим:

Обратная подстановка (back substation) для уравнения осуществляется аналогично, начиная с последнего уравнения, и так же требует времени Матрица – верхне-треугольная, и система имеет вид:

Мы последовательно вычисляем по формулам

или, в общем виде:

Алгоритм решения задачи

Первым делом необходимо считать входные данные. Для этого опишем ввод переменных. Сначала считываем – колличество уравнений и неизвестных в них. Затем в двойномом цикле от 0 до -1 считываем массив – массив коэффициентов при неизвестных. И наконец вектор – столбец свободных членов (рис. 1).

Рисунок 1 – Схема считывания входных данных

Напишем процедуру LUP-Decomposition (рис. 2), осуществляющую LUP – разложение матрицы . Вместо рекурсивного вызова функции будем использовать циклы. Матрицу перестановок представим массивом , считая, что в — й строке единица стоит на -м месте. Матрицы и будут вычисляться на месте матрицы , и после работы процедуры

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

Рисунок 2 – Схема LUP – разложения

Теперь решим систему с помощью прямой и затем обратной подстановкой. Для этого опишем процедуру LUP-Solve (рис. 3). В ней, по заданным , ищется , делая прямую и затем обратную подстановку. Сначала делаем прямую подстановку для определения , и затем обратную для . Внутри каждого из двух циклов присутствуют циклы суммирования, от этого сложность алгоритма равна

Рисунок 3 – Схема прямой и обратной подстановки

Рисунок 4 – Схема вывода ответа

Осталось вывести ответ СЛАУ, а именно . Для этого организуем цикл от первого до последнего элемента (рис. 4).

В качестве примера решим систему , где

,

а – вектор неизвестных [3]. Рассмотрим LUP – разложение матрицы :

(a) Изначально имеется исходная матрица и тождественная перестановка в столбце слева (одномерное представление матрицы ). На первом шаге ищем ведущий элемент, у нас он равен 5 и находится на третьей строчке. (б) Переставляем строки 1 и 3, одновременно меняя их в столбце перестановок. Выделенные столбец и строка представляют собой векторы и . (в) Вектор заменен на вектор , так же было вычислено дополнение Шура. Теперь матрица состоит из 3 частей: элементы сверху, элементы слева, и дополнение Шура (остальная часть матрицы). (г) Повторяем шаги а и б. (д) Ищем и дополнение Шура. На последнем шаге никаких изменений не происходит. (и) Разложение матрицы на .

Прямой подстановкой найдем из системы . Система имеет вид [1]:

.

Вычисляя сначала , затем , и, наконец, , получаем

.

Теперь обратной подстановкой найдем из уравнения , которое имеет вид

.

Вычисляя последовательно , , и , получим

.

Тестирование показало корректность работы алгоритма.

Заключение

Для решения СЛАУ методом LUP-разложения в курсовом проекте было поставлено и решено множество целей. Основными из них являются: разложение матрицы и прямая/обратная подстановка. Алгоритм имеет как достоинства, так и недостатки. Основным недостатком является сложность алгоритма равная и отражающаяся на времени решений больших матриц, но преимущество алгоритма заключается в высокой устойчивости при решении СЛАУ, что позволяет получить наиболее точный результат, чем при использовании других алгоритмов.

Достоверное и точное решение систем алгебраических и дифференциальных уравнений на языке Си Текст научной статьи по специальности « Компьютерные и информационные науки»

Аннотация научной статьи по компьютерным и информационным наукам, автор научной работы — Маничев В. Б., Жук Д. М., Глазкова В. В.

Приведены основные утверждения и научные результаты многолетних работ в области численного моделирования динамических систем, положенные в основу библиотеки математических программ SADEL (Sets of Algebraic and Differential Equations solvers Library) Си-библиотеки, предназначенной для решения алгебраических и дифференциальных уравнений с максимально возможной компьютерной точностью. Показаны результаты решения тестовых и практических задач, подтверждающие эффективность разработанных программ библиотеки SADEL

Похожие темы научных работ по компьютерным и информационным наукам , автор научной работы — Маничев В. Б., Жук Д. М., Глазкова В. В.

RELIABLE AND ACCURATE SOLVING OF SYSTEMS OF ALGEBRAIC AND DIFFERENTIAL EQUATIONS IN C LANGUAGE

Basic statements and scientific results of many-year works in the field of numerical simulation of dynamic systems are given that underlie the mathematical program library SADEL (Sets of Algebraic and Differential Equations solvers Library). SADEL is the C language library destined for solving algebraic and differential equations with maximum possible computer precision. Results of solving the test and practical problems are given that confirm the effectiveness of the SADEL programs

Текст научной работы на тему «Достоверное и точное решение систем алгебраических и дифференциальных уравнений на языке Си»

ИНФОРМАТИКА И ВЫЧИСЛИТЕЛЬНАЯ

ДОСТОВЕРНОЕ И ТОЧНОЕ РЕШЕНИЕ СИСТЕМ АЛГЕБРАИЧЕСКИХ И ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ НА ЯЗЫКЕ СИ

В.Б. Маничев1, Д.М. Жук1, В.В. Глазкова2

ХМГТУ им. Н.Э. Баумана, Москва, Российская Федерация e-mail: manichev@bmstu.ru; zhuk@bmstu.ru

2МГУ им. М.В. Ломоносова, Москва, Российская Федерация e-mail: glazkova@mlab.cs.msu.su

Приведены основные утверждения и научные результаты многолетних работ в области численного моделирования динамических систем, положенные в основу библиотеки математических программ SADEL (Sets of Algebraic and Differential Equations solvers Library) — Си-библиотеки, предназначенной для решения алгебраических и дифференциальных уравнений с максимально возможной компьютерной точностью. Показаны результаты решения тестовых и практических задач, подтверждающие эффективность разработанных программ библиотеки SADEL.

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

RELIABLE AND ACCURATE SOLVING OF SYSTEMS OF ALGEBRAIC AND DIFFERENTIAL EQUATIONS IN C LANGUAGE

V.B. Manichev1, D.M. Zhuk1, V.V. Glazkova2

1Bauman Moscow State Technical University, Moscow, Russian Federation e-mail: manichev@bmstu.ru; zhuk@bmstu.ru

2Lomonosov Moscow State University, Moscow, Russian Federation e-mail: glazkova@mlab.cs.msu.su

Basic statements and scientific results of many-year works in the field of numerical simulation of dynamic systems are given that underlie the mathematical program library SADEL (Sets of Algebraic and Differential Equations solvers Library). SADEL is the C language library destined for solving algebraic and differential equations with maximum possible computer precision. Results of solving the test and practical problems are given that confirm the effectiveness ofthe SADEL programs.

Keywords: simulation, dynamic systems, ordinary differential equations, differential algebraic equations, integration methods, linear algebraic equations.

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

из причин стали существенные погрешности компьютерных вычислений, результаты которых попали в конечную конструкторскую документацию. Главным недостатком вычислительных методов, алгоритмов и программ, используемых при проектировании и управлении и реализованных в лучших зарубежных программно-методических комплексах (ПМК) моделирования динамических систем, является возможная выдача ошибочных результатов компьютерных вычислений без предупреждения пользователей об их недостоверности. Подавляющее большинство универсальных и, особенно, специализированных программ численного моделирования разнообразных динамических систем и объектов разработано и разрабатывается с использованием языка Си (Си++). При этом пользователи ПМК моделирования динамических систем часто сами пишут программы решения систем обыкновенных дифференциальных и дифференциально-алгебраических уравнений (ОДУ-ДАУ), а также линейных алгебраических уравнений (ЛАУ), считая это простой задачей, а затем разрабатывают прекрасные интерфейсы и визуализацию неверных результатов решения сложных вычислительных задач. Библиотека SADEL (Sets of Algebraic and Differential Equations solvers Library) разработана как библиотека программ-решателей для жестких систем ОДУ-ДАУ и плохо обусловленных систем ЛАУ и предназначена для решения указанных проблем при численном моделировании динамических систем [1, 2]. В настоящей работе используется следующая терминология.

Математическая модель (mathematical model) — система ОДУ-ДАУ, отображающая функционирование моделируемой динамической системы с требуемой достоверностью и точностью.

Математическое моделирование (mathematical modeling) — процесс получения математических моделей динамических систем на основе фундаментальных законов физики, химии и др.

Компьютерная модель — алгоритмическая и программная реализация математической модели.

Численное (компьютерное) моделирование (далее для сокращения просто моделирование (simulation)) — процесс получения результатов математического моделирования динамических систем с помощью компьютерных моделей.

Жесткие системы ОДУ-ДАУ (stiff ODE-DAE systems (stiff set of ODE-DAE)) — системы ОДУ-ДАУ со степенью жесткости более 106.

Плохо обусловленные системы ЛАУ (badly conditioned (ill-conditioned) LAE systems (ill-conditioned set of LAE)) — системы ЛАУ с числом обусловленности более 106, к решению которых часто сводится решение жестких и сверхжестких систем ОДУ-ДАУ неявными методами интегрирования.

Достоверность численного моделирования (veracity of simulation) — отклонение результатов численного моделирования динамических систем от заведомо точных результатов не более чем на 20 % (точные результаты получают аналитически или экспериментально, как правило, оценивая относительную погрешность численного моделирования в процентах).

Точность численного моделирования (accuracy of simulation) — фактическая степень отклонения результатов достоверного численного моделирования от заведомо точных результатов.

Математическая точность (mathematical accuracy) решения систем ОДУ-ДАУ — заданная относительная погрешность нормы векторов решения систем ОДУ-ДАУ (параметр eps) по сравнению с заведомо точным аналитическим решением в относительных единицах (например, в ПМК MATLAB по умолчанию eps = 0,001).

Компьютерная точность (precision) решения систем ОДУ-ДАУ и ЛАУ — число полученных верных значащих цифр в каждом элементе векторов решения систем ОДУ-ДАУ и ЛАУ по сравнению с заведомо точным аналитическим (идеальным для систем ЛАУ) решением.

Решение систем ЛАУ с повышенной точностью вычислений (LAE systems solution with extra precision calculations) — применение методов получения решений с повышенной точностью вычислений (extra precision calculations), реализованных в соответствующих ПМК, а также в библиотеках математических программ и описанных в документации на соответствующие ПМК и библиотеки. Возможно также использование переменных типа Decimal на языке C# или класса java.math.BigDecimal на языке Java.

Далее приведены методы для получения решений с повышенной точностью вычислений, реализованные в наиболее известных ПМК: Maple (http://maplesoft.com — метод Software Floating Point (SFP метод)), MATLAB (http://mathworks.com — метод Variable Precision Arithmetic (VPA метод)), Mathematica (http://wolfram.com — метод Arbitrary Precision Arithmetic (APA метод));

библиотеки математических программ для получения решений с повышенной точностью вычислений: Intel для десятичных вычислений с повышенной точностью (стандарт IEEE 754-2008 Decimal Floating-Point for Intel§Architecture Processors, IEEE Standard for Float-ing-Point Arithmetic 7542008, IEEE, 2008:

http://ieeexplore.ieee.org/xpl/loginjsp?tp=&arnumber=5223319&url [3]), а также Си (Си++)-библиотеки для вычислений с повышенной точностью: http://en.wikipedia.org/wiki/Arbitrary-precision_arithmetic http://en.wikipedia.org/wiki/GNU_Multi-Precision_Library http://en.wikipedia.org/wiki/MPFR http://en.wikipedia.org/wiki/Class_Library_for_Numbers http://gmplib.org/

http://www. boost.org/doc/libs/1_53_0/libs/rational/rational.html http://speleotrove.com/decimal/dnusers.html#example7 http://www. ttmath.org/.

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

общем случае (без использования ПМК визуального математического моделирования типа ANSYS, NX-CAE (NASTRAN), SIMULIA (ABAQUS), MATLAB-SIMULINK, Maple-MapleSim, Mathematica-SystemModeler и др.), можно представить как последовательность четырех основных этапов.

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

2. Для выбранного алгоритмического языка определение численных методов и алгоритмов решения или выбор из соответствующей библиотеки готовой программы-решателя уравнений математической модели (в этом случае третий этап не нужен).

3. Разработка на выбранном алгоритмическом языке и отладка расчетной программы-решателя уравнений математической модели.

4. Ввод исходных данных для разработанной или выбранной программы-решателя уравнений математической модели, получение и обработка результатов численного моделирования.

Соответственно можно выделить четыре основных источника ошибок в процессе численного моделирования.

1. Погрешность исходной математической модели. Это очень важный источник ошибок, так как он является неустранимым. Например, второй закон Ньютона в дифференциальной форме не всегда дает корректные результаты математического моделирования, даже при отсутствии других источников ошибок математического моделирования. Анализ этих ошибок выходит за рамки настоящей статьи, поскольку относится к разным предметным областям математического моделирования (multi-physics, multi-discipline simulation).

2. Ошибки, связанные с неверным выбором численных методов и алгоритмов. Это устранимый источник ошибок, так как пользователям ПМК численного моделирования динамических систем не следует самому программировать решение вычислительных задач, а следует использовать готовые программы-решатели уравнений математических моделей из соответствующих ПМК или известных библиотек математических программ. Это устранимый источник ошибок математического моделирования.

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

4. Ошибки, связанные с вводом неверных исходных данных. Это устранимый источник ошибок математического моделирования.

Результатом многолетних исследований проблем численного моделирования динамических систем являются следующие утверждения.

Математические модели динамических технических систем необходимо получать на основе фундаментальных физических законов только в форме систем ДАУ, не разрешенных относительно производных, и решать эти системы без каких-либо преобразований и без получения «правой» части для производных в явном, аналитическом виде. При математическом моделировании динамических процессов в реальных технических системах и объектах состояние этих систем и объектов обычно рассматривают в пространстве дифференциальных переменных или в пространстве «переменных состояния» с получением в явном аналитическом виде «правой» части для производных этих переменных по времени, но этого недостаточно для достоверного и точного математического моделирования реальных технических систем и объектов, так как любой дифференциальной переменной состояния в технических системах обычно соответствует алгебраическая переменная. Например, в электрике и электронике дифференциальным переменным состояния — напряжениям на емкостях и токам через индуктивности — соответствуют алгебраические переменные — токи через емкости и напряжения на индуктив-ностях, в механике дифференциальным переменным состояния — скоростям тел определенной массы — соответствуют алгебраические переменные — силы инерции и другие параметры, поэтому состояние таких динамических систем и объектов следует рассматривать в пространстве дифференциально-алгебраических переменных. Основным физическим свойством алгебраических переменных является возможность идеального скачка значения этих переменных в бесконечно малый отрезок времени, в то время как для дифференциальных переменных состояния такие скачки невозможны, поэтому приведение систем ДАУ к системам ОДУ в нормальной форме Коши с преобразованием алгебраических переменных в дифференциальные переменные состояния физически неверно и полученная таким образом математическая модель динамических процессов в реальных технических системах и объектах будет физически недостоверна.

Эффективность численных методов решения систем ОДУ-ДАУ, НАУ и ЛАУ включает в себя два аспекта:

1) достоверность и точность получаемых решений данных уравнений;

2) затраты машинного времени и оперативной памяти на получение решений.

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

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

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

— при базовом, начальном математическом и компьютерном моделировании динамических технических систем не требуется высокая математическая точность выдаваемых пользователю результатов, так как параметры математических моделей этих систем (и соответствующие коэффициенты систем ОДУ-ДАУ и ЛАУ) получены, как правило, экспериментально с невысокой математической точностью. Поэтому заданная математическая точность конечных результатов численного моделирования по умолчанию может быть невысокой, но она должна быть гарантировано обеспечена.

Программа-решатель систем ДАУ в библиотеке 8АББЬ рассматривает системы ДАУ в еще более расширенном координатном пространстве переменных, в котором к дифференциально-алгебраическим переменным добавлены производные дифференциальных переменных [4]:

где X — вектор координатного базиса дифференциальных переменных состояния размерностью т, У — вектор алгебраических переменных полного координатного базиса дифференциально-алгебраических переменных размерностью к, РХ = йХ/йг — вектор производных дифференциальных переменных по времени для расширенного координатного пространства переменных размерностью т,Ь — независимая переменная (обычно время); О — вектор-функция математической модели динамической системы размерностью т + к. Заданы начальные условия Х0 = X(0) и отрезок интегрирования

Новые алгоритмы для решения приведенной системы основаны на совместном решении систем НАУ Ог(Х,У, РХ,гг) = 0 и систем из т алгебраических уравнений численного интегрирования Н^ (X, РХ, Н) = 0,

сформированных на соответствующих стадиях методов интегрирования, относительно дифференциально-алгебраических переменных Х,У и производных дифференциальных переменных РХ (DABC методы Маничева-Глазковой) [2]:

= Нп йгзРХз — аг3Х] — ЬгХп-1 — НпСгРХп-1 = 0;

Хп = Х3; РХп = РХ3, г = 1. ,з, = з — ; Нп = (гп — гп-\) — п-й шаг интегрирования; йгз, (щ, Ьг,сг — параметры метода.

Параметры алгебраических уравнений численного интегрирования И (Х,РХ,Н) = 0 были получены на основе анализа пяти базовых классов задач для фундаментальных решений линейных неоднородных систем ОДУ. Было показано, что для достоверного решения всех пяти классов фундаментальных решений при больших значениях шага интегрирования (т.е. при невысокой математической точности интегрирования) численные методы решения систем ОДУ-ДАУ должны быть ЛЬ-устойчивыми, т.е. абсолютно (А-) устойчивыми строго в левой (Ь-ЬеА) полуплоскости комплексной плоскости устойчивости методов численного решения систем ОДУ-ДАУ [5]. При невысокой математической точности интегрирования шаг интегрирования может стать таким большим, что не ЛЬ-устойчивые методы интегрирования могут потерять устойчивость для некоторых фундаментальных решений. Область абсолютной устойчивости ЛЬ-устойчивых методов интегрирования показана на рис. 1.

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

Численное решение систем ОДУ-ДАУ во многих численных методах их решения сводится к многократному численному решению соответствующих систем нелинейных алгебраических уравнений (НАУ) и ЛАУ [5]. Высокий порядок таких систем и необходимость их многократного решения в процессе математического моделирования предъявляют жесткие требования к эффективности применяемых численных методов решения систем ОДУ-ДАУ и особенно ЛАУ.

Для решения жестких систем ОДУ-ДАУ, которым часто соответствует плохо обусловленная матрица Якоби вдоль траектории решения, рекомендуется использовать неявные методы интегрирования, которые сводятся к решению систем НАУ на каждом шаге интегрирования. Для получения решений систем НАУ с гарантированной компьютерной точностью эти системы следует решать только такими итерационными методами, которые в «точке» решения системы НАУ сводятся к решению системы ЛАУ с матрицей Якоби в качестве матрицы коэффициентов такой системы ЛАУ. Если не получать решения систем НАУ с гарантированной компьютерной точностью, то траектории решения систем ОДУ-ДАУ могут быть ошибочными. В библиотеке 8ЛББЬ используется метод Ньютона-Рафсона и решение систем НАУ и ЛАУ с гарантированной удвоенной компьютерной точностью языка Си (Си++) на каждом шаге интегрирования, т.е. с гарантированной компьютерной точностью (15 верных значащих цифр для всех элементов вектора

Рис. 1. Область абсолютной устойчивости (заштрихована) AL-устойчивых методов интегрирования

решений систем ЛАУ при вычислениях с удвоенной точностью на языке Си). Аналитически данная проблема не имеет решения, например, теоретически умножение квадратной матрицы с постоянными коэффициентами на обратную равно единичной матрице, но для плохо обусловленных матриц при вычислении на компьютере единичную матрицу мы не получим. Только многолетние вычислительные эксперименты позволили решить эту проблему в программах-решателях систем ЛАУ библиотеки SADEL на основе перечисленных методов вычислений с повышенной точностью.

При разработке математических программ на языке Си необходимо также учитывать особенности программирования арифметических операций с вещественными числами с плавающей точкой, так как именно такие операции используют упомянутые ПМК и библиотеки математических программ при решении систем ОДУ-ДАУ и ЛАУ. Институт инженеров по электротехнике и электронике (Institute of Electrical and Electronic Engineers, IEEE) разработал международные стандарты, которые описывают представление вещественных чисел с плавающей точкой:

— стандарт ANSI/IEEE 754:1985 определяет требования к реализации двоичной плавающей арифметики;

— стандарт ANSI/IEEE 754:2008 обобщает стандарт ANSI/IEEE 754:1985, допуская дополнительно, кроме двоичного, десятичное основание представлений мантиссы и экспоненты и произвольную длину машинного слова.

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

Стандарты IEEE754 определяют следующие основные форматы представления вещественных чисел с плавающей точкой.

1. С одинарной точностью (соответствует типам SINGLE PRECISION (32 двоичных разряда) или REAL*4 в языке Фортран и float в языке Си).

2. С удвоенной точностью (соответствует типам DOUBLE PRECISION (64 двоичных разряда) или REAL*8 в языке Фортран и double в языке Си), этот тип представления принят по умолчанию во всех упомянутых ПМК и стандартных библиотеках математических программ.

3. С расширенной точностью (условно говоря, соответствует типам REAL*10 (80 двоичных разрядов) в языке Фортран и long double в языке Си).

Формат с учетверенной точностью QUAD PRECISION (128 двоичных разряда) или REAL*16 в языке Фортран стандарт IEEE754 отдельно не определяет, поэтому в языке Си следует создавать специальные форматы для данного представления вещественных чисел.

При реализации вычислительных методов решения жестких систем ОДУ-ДАУ и плохо обусловленных систем ЛАУ не следует использовать вычисления с одинарной точностью, причем при решении плохо обусловленных систем ЛАУ следует использовать методы вычислений с повышенной точностью, если вычисления с общепринятой удвоенной точностью не обеспечат компьютерной точности результатов. При программной реализации формул неявного интегрирования жестких систем ОДУ-ДАУ не должно быть также арифметических операций деления на шаг интегрирования, чтобы не было ограничений из-за «машинного нуля» на значение минимального шага интегрирования [6].

Результаты решения тестовых и практических задач. На основе приведенных утверждений и научных результатов была разработана библиотека SADEL — библиотека программ-решателей на языке Си для численного решения жестких систем ОДУ-ДАУ и плохо обусловленных систем ЛАУ с гарантированной достоверностью и компьютерной точностью. В результате при невысоких требованиях к математической точности конечных результатов численного моделирования получаются только корректные, достоверные и точные результаты решения данных систем по сравнению с соответствующими программами-решателями из зарубежных библиотек и пакетов математических программ: LAPACK (Linpack), MAGMA, NAG, IMSL, Intel MKL, MATLAB, Mathcad, Maple, Mathematica), что показано в работе [2]. Например, с помощью данных программ удалось решить с гарантированной компьютерной точностью ряд тестовых плохо обусловленных систем ЛАУ, которые не решаются с обычной удвоенной точностью языка Си другими известными программами-решателями систем ЛАУ из указанных библиотек и ПМК. В качестве примера рассмотрим тестовую задачу с матрицей Гильберта (Hilbert matrix) 10-го порядка в качестве матрицы коэффициентов системы ЛАУ. Все элементы матрицы коэффициентов и вектора правых частей, как суммы соответствующих строк матрицы коэффициентов, были вычислены с удвоенной точностью языка Си. Абсолютно точное решение такой задачи — единичный вектор. При решении с удвоенной точностью этой тестовой задачи с матрицей 10-го порядка с помощью программы-решателя систем ЛАУ из библиотеки MAGMA (LAPACK-Linpack) был получен следующий результат: 1.000000000464697 0.999999960244827 1.000000840984605 0.999992393461859 1.000036137385782 0.999900979283788 1.000162025520190 0.999843779478675 1.000081852296234 0.999982030574678

Курсивом выделены неверные значащие цифры в решении. Соответствующие программы-решатели из библиотеки SADEL получили единичный вектор для этой тестовой задачи с точностью в 15 верных значащих цифр для всех элементов вектора решений данной системы ЛАУ.

Сравнение программы-решателя систем ДАУ в библиотеке SADEL с зарубежными программами-решателями систем ОДУ-ДАУ было проведено

для параметров методов интегрирования программ-решателей систем ОДУ-ДАУ, рекомендуемых для этих программ-решателей по умолчанию. Все программы-решатели систем ОДУ-ДАУ сравнивались при заданной невысокой относительной точности интегрирования eps = 0,001 (в ПМК MATLAB это точность, задаваемая по умолчанию). В таблице приведены результаты сравнения только наиболее трудных для современных программ-решателей систем ОДУ-ДАУ жестких систем ОДУ с многопериодным решением. Если программа-решатель систем ОДУ-ДАУ не решит достоверно и точно хотя бы одну из приведенных в таблице тестовых задач, то эту программу-решатель не следует использовать для моделирования динамических систем, математическая модель которых представляет собой жесткую систему ОДУ-ДАУ с многопериодным решением [2].

В таблице сравнивались только методы интегрирования программ-решателей систем ОДУ-ДАУ, рекомендуемые в соответствующих ПМК и библиотеках математических программ для решения жестких систем ОДУ-ДАУ. Знак минус означает невозможность получения решения или (в большинстве случаев) качественно неверный результат без всякого предупреждения пользователя о возможных ошибках. Например, на рис. 2 приведен результат качественно неверного решения теста 4 с помощью программы-решателя систем ОДУ-ДАУ неявным методом BDF из библиотеки C-Library NAG (а), а также правильное решение этого теста с помощью методов М2 и М3 из библиотеки SADEL (б). Предупреждений о возможном недостоверном решении в программе-решателе библиотеки C-Library NAG не было.

Далее приведен пример решения практической задачи химической кинетики [7]. Данная задача является жесткой за счет большого разброса порядков в константах скоростей реакции и собственных значений матрицы Якоби, коэффициент жесткости более 106. Для решения этой задачи сначала был выбран пакет математических программ Mathematica 8.0. Для проверки достоверности и точности полученного решения была выбрана Си-библиотека программ-решателей жестких систем ОДУ-ДАУ и плохо обусловленных систем ЛАУ SADEL.

Как следует из рис.3,а, метод, выбранный в пакете Mathematica 8.0 по умолчанию выдал колебательный график решения. Неявный метод М3, реализованный в библиотеке SADEL, выдал правильный график решения — рис. 3, б. Следует отметить, что после настройки параметров методов интегрирования, реализованных в пакете Mathematica 8.0, был также получен правильный график решения, однако большинство пользователей математических программ не являются профессиональными специалистами в численных методах и программах для решения систем ОДУ-ДАУ, поэтому достоверность и требуемая точность решения систем ОДУ-ДАУ должна быть обеспечена для параметров программ-решателей этих уравнений, рекомендуемых для этих программ-решателей по умолчанию.

Заключение. Библиотека Си программ SADEL для решения жестких систем ОДУ-ДАУ и плохо обусловленных систем ЛАУ является практическим

Сравнение программ-решателей систем ОДУ-ДАУ для жестких систем ОДУ с многопериодным решением

Тесты Шрограммы-решатели систем ОДУ-ДАУ 8АБЕЬ-Си библиотека 2010. Методы М2,МЗ МАТЬАВ 2007. Метод Оёе15в Мар1е 2007. Метод ЯозепЬгоск МаИтетаИса 8 2011. Метод ВОБ МАО С-ЫЬгагу 2012. Метод е 1’>Г)К ШБЬ С-ЫЬгагу 2012. Метод ВОБ

ТЕСТ 1. Уравнения Ван дер Поля ми = 10е + + + +

ТЕСТ 2. Уравнения Ван дер Поля ми = ю9 +

ТЕСТ 3. Высокодобротный фильт Ы = 1, кг = 1, ки = 0,01 + + +

ТЕСТ 4. Высокодобротный фильтр ы = ю-104, ки = 1, Ы = 1 + +

ТЕСТ 5. Локально неустойчивая система ОДУ ми = 10е + +

ТЕСТ 6. Моделирование свечения лазера + + + +

Рис.2. Решение задачи ТЕСТ 4 при TOL = 0,001:

а — метод BDF из библиотеки C-Library NAG, б — метод М3 из Си библиотеки SADEL

результатом многолетних научно-исследовательских работ авторов в качестве математического ядра платформы математического и компьютерного моделирования разнородных (multi-physics, multi-discipline) динамических технических систем и объектов (программный комплекс ПА10 (SADEL-PA10)) [2], превосходящей подобные зарубежные программные продукты MATLAB-SIMULINK, Maple-MapleSim, C-Library NAG и C-Library IMSL, Intel MKL в части решения жестких систем ОДУ-ДАУ и плохо обусловленных систем ЛАУ со стандартной удвоенной точностью вычислений языка Си.

Разработка программ-решателей систем ОДУ-ДАУ и ЛАУ на языке Си для библиотеки SADEL позволяет также сделать следующие выводы.

1. Корректные математические модели динамических процессов в реальных технических системах и объектах необходимо получать в расширенном дифференциально-алгебраическом координатном базисе переменных на основе фундаментальных физических законов в форме неоднородных систем ДАУ, не разрешенных относительно производных, и решать эти системы без каких-либо эквивалентных математических преобразований и без получения «правой» части для производных в явном, аналитическом виде.

2. Для достоверного и точного решения систем ДАУ при невысоких требованиях к математической точности результатов следует использовать толь-

Рис.3. График зависимости дифференциальной переменной х8 системы ОДУ от времени при TOL = 10-6: а — МаЛешаНеа 8.0; б — 8ЛББЬ

ко ЛЬ-устойчивые методы интегрирования, при использовании других методов интегрирования необходимо гарантировать устойчивость численного решения при увеличении шага интегрирования в случае невысоких значений параметра ерв.

Новые научные результаты предполагается получить в направлениях разработки новых высокоэффективных неявных БЛБС методов интегрирования с полностью заполненными матрицами Б и Л в новых системах алгебраических уравнений интегрирования, разработки новых методов и алгоритмов автоматического выбора шага интегрирования для ЛЬ-устойчивых неявных методов интегрирования, а также эффективных алгоритмов и программ решения систем ЛАУ высокой и сверхвысокой размерности с повышенной точностью вычислений.

Работа выполнена при поддержке РФФИ, гранты № 12-07-00324/12 и № 12-07-31029.

1. Андронов А.В., ЖукД.М., Кожевников Д.Ю., Маничев В.Б. Библиотека математических программ-решателей на языке Си: SADEL // http://www.pa10.ru (дата обращения 20.10.2013).

2. Жук Д.М., Маничев В.Б., Сахаров М.К. SADEL — библиотека «сверхточных» решателей для программного комплекса ПА10 (SADEL-PA10) // Сб. науч. трудов МЭС-2012. М.: ИППМ РАН, 2012. С. 147-153.

3. IEEE Standard for Floating-Point Arithmetic 754-2008, IEEE, 2008 // http://iee-explore.ieee.org/xpl/login.jsp?tp=&arnumber=5223319&url=http%3A%2F%2Fiee-explore.ieee.org%2Fxpls%2Fabs_all.jsp%3Farnumber%3D5223319 (дата обращения 21.10.2013).

4. Норенков И.П., Трудоношин В.А., Федорук В.Г. Метод формирования математических моделей для адаптируемых программных комплексов анализа радиоэлектронных схем // Радиотехника, 1986. № 9. С. 67-72.

5. Жук Д.М., Маничев В.Б., Ильницкий А.О. Методы и алгоритмы решения дифференциально-алгебраических уравнений для моделирования систем и объектов во временной области // Информационные технологии. 2010. Ч. 1. № 7. С. 16-24. Ч. 2. № 8. С. 23-26.

6. Маничев В.Б. Новые алгоритмы для программ анализа динамики технических систем // Вестник МГТУ им. Н.Э. Баумана. Сер. Приборостроение. 1996. № 1. С. 48 -56.

7. Valerii N. Snytnikov, Tatyana I. Mishchenko, Vladimir N. Snytnikov, Sergei E. Malykhin, Vasilii I. Avdeev, Valentin N. Parmon. Autocatalytic gas-phase dehydrogenation of ethane // Research on Chemical Intermediates, March 2012. Vol. 38. Issue 3-5. Р. 1133-1147.

[1] Andronov A.V., Zhuk D.M., Kozhevnikov D.Yu., Manichev V.B. Biblioteka

matematicheskikh programm-reshateley na yazyke Si: SADEL [Library of mathematical problem solvers in C: SADEL] Available at: http://www.pa10.ru (accessed 14.01.14) (in Russ.).

[2] Zhuk D.M., Manichev V.B., Sakharov M.K. SADEL — library «extremely-precise»

problem solvers for complex software PA10 (SADEL-PA10). Sb. nauch. trudov VI Vserossiyskoy nauchno-tekhn. konf. «Problemy razrabotki perspektivnykh mikro- i nanoelektronnykh sistem MES-2012 [Call. Papers of All-Russ. Sc.-Techn. Со^. «Development problems of advanced micro-and nanoelectronic systems»], Moscow, IPPM RAN Publ., 2012. pp. 147-153 (in Russ.).

[3] IEEE Standard for Floating-Point Arithmetic 754-2008, IEEE, 2008. Available at:

http://ieeexplore.ieee.org/xpl/login.jsp?tp=&arnumber=5223319&url=http%3A% 2F%2Fieeexplore.ieee.org%2Fxpls%2Fabs_all.jsp%3Farnumber%3D5223319 (accessed 28.12.2013).

[4] Norenkov I.P., Trudonoshin V.A., Fedoruk V.G. Method of mathematical models

formation for adaptable software systems of electronic circuits analysis. Radiotekhnika [Radioengineering], 1986, no. 9, pp. 67-72 (in Russ.).

[5] Zhuk D.M., Manichev V.B., Il’nitskiy A.O. Methods and algorithms for solving

differential-algebraic equations for simulation systems and objects in the time domain. Informatsionnye tekhnologii [Information technologies], 2010, no. 7, pp. 16-24 (ch. 1); no. 8, pp. 23-26 (ch. 2) (in Russ.).

[6] Manichev V.B. New algorithms for program analysis of the dynamics of technical

systems. Vestn. Mosk. Gos. Tekh. Univ. im. N.E. Baumana, Priborostr. [Herald of the Bauman Moscow State Tech. Univ., Instrum. Eng.], 1996, no. 1, pp. 48-56 (in Russ.).

[7] Snytnikov V.N., Mishchenko T.I., Snytnikov V.N., Malykhin S.E., Avdeev V.I., Parmon V.N. Autocatalytic gas-phase dehydrogenation of ethane. Research on Chem. Int., 2012, vol. 38, iss. 3-5, pp. 1133-1147.

Статья поступила в редакцию 26.11.2012

Владимир Борисович Маничев — канд. техн. наук, доцент кафедры «Системы автоматизированного проектирования» МГТУ им. Н.Э. Баумана. Автор более 80 научных работ в области математического обеспечения систем автоматизированного проектирования.

МГТУ им. Н.Э. Баумана, Российская Федерация,105005, Москва, 2-я Бауманская ул., д. 5.

V.B. Manichev — Cand. Sci. (Eng.), assoc. professor of «Systems of Computer-Aided Design» department of the Bauman Moscow State Technical University. Author of more than 80 publications in the field of mathematical support of computer-aided design systems. Bauman Moscow State Technical University, Vtoraya Baumanskaya ul. 5, Moscow, 105005 Russian Federation.

Дмитрий Михайлович Жук — канд. техн. наук, доцент кафедры «Системы автоматизированного проектирования» МГТУ им. Н.Э. Баумана. Автор более 80 научных работ в области математического обеспечения систем автоматизированного проектирования.

МГТУ им. Н.Э. Баумана, Российская Федерация,105005, Москва, 2-я Бауманская ул., д. 5.

D.M. Zhuk — Cand. Sci. (Eng.), assoc. professor of «Systems of Computer-Aided Design» department of the Bauman Moscow State Technical University. Author of more than 80 publications in the field of mathematical support of computer-aided design systems. Bauman Moscow State Technical University, Vtoraya Baumanskaya ul. 5, Moscow, 105005 Russian Federation.

Валентина Владимировна Глазкова — канд. физ.-мат. наук, ассистент кафедры «Автоматизация систем вычислительных комплексов» факультета «Вычислительная математика и кибернетика» МГУ им. М.В. Ломоносова. Автор более 10 научных работ в области вычислительной математики, программирования и систем автоматической категоризации текстовой информации.

МГУ им. М.В. Ломоносова, Российская Федерация, 119991, Москва, ГСП-1, Ленинские горы, д. 1, стр. 52.

V.V. Glazkova — Cand. Sci. (Phys.-Math.), assistant lecturer of «Automation of Computing Complex Systems» department of «Computational Mathematics and Cybernetics» faculty of the Lomonosov Moscow State University. Author of more than 10 publications in the field of computational mathematics, programming, and systems of automatic categorization of textual information.

Lomonosov Moscow State University, Leninskie Gory 1, stroenie 52, Moscow, 119991 Russian Federation.