Всероссийская научно-техническая конференция "Проблемы разработки перспективных микро- и наноэлектронных систем" МЭС-2012 к. т. н. , доцент Маничев Владимир Борисович к. т. н. , доцент Жук Дмитрий Михайлович Студент Сахаров Максим Константинович Кафедра РК 6 (САПР) МГТУ им. Н. Э. Баумана http: //rk 6. bmstu. ru 2012 г.
ТЕМА ДОКЛАДА: SADEL (SETS OF DIFFERENTIAL AND ALGEBRAIC EQUATIONS SOLVERS LIBRARY): СИ библиотека «сверхточных» программрешателей (extra precision solvers) систем ОДУДАУ и ЛАУ для программного комплекса ПА 10 (SADEL-PA 10) Сайт проекта ПА 10 (SADEL-PA 10) www. рa 10. ru или www. па 10. рф
ПРОЕКТ «ПА 10 (SADEL-PA 10)» РАЗРАБОТКА ПРОГРАММНО-МЕТОДИЧЕСКОГО КОМПЛЕКСА ДЛЯ МАТЕМАТИЧЕСКОГО И КОМПЬЮТЕРНОГО МОДЕЛИРОВАНИЯ РАЗНОРОДНЫХ (MULTI-PHYSICS, MULTI-DISCIPLINE) ДИНАМИЧЕСКИХ СИСТЕМ С ОТКРЫТЫМ МАТЕМАТИЧЕСКИМ ЯДРОМ: СИ БИБЛИОТЕКА SADEL
ХРОНОЛОГИЯ РАЗРАБОТКИ ПРОГРАММ АНАЛИЗА (ПА) ДИНАМИЧЕСКИХ СИСТЕМ НА ОСНОВЕ РЕШЕНИЯ СИСТЕМ ОДУ-ДАУ Кафедра РК 6 (САПР) МГТУ им. Баумана ПА 1 (Норенков И. П. 1974 г. ) – (ПАЭС) Программа Анализа Электронных Схем (ассемблер) ПА 2 (Иванов С. Р. 1975 г. ) – анализ источников питания (фортран) ПА 3 (Хартов В. Я. 1976 г. ) – анализ цифровых блоков (фортран) ПА 4 (Жук Д. М. 1978 г. ) – анализ радиационной стойкости электронных схем (ассемблер) ПА 6 (Федорук В. Г. , Трудоношин В. А. 1983 г. ) – анализ разнородных (multi-physics, multi-discipline) объектов (ассемблер) ПА 7 (Уваров М. Ю. , Маничев В. Б. 1989 г. ) – анализ разнородных объектов (С - DOS) ПА 8 (Федорук В. Г 1994 г. ) – анализ разнородных объектов (C - Unix) ПА 9 (Норенков И. П. , Уваров М. Ю. , Трудоношин В. А. 1998 г. ) – анализ разнородных объектов (Java - Windows, Linux): Регистрационное свидетельство № 10634 от 4 июля 2007 г. ПА 10 -PAX (Алиев Е. А. , Баженов В. В. , Пузырев Е. А, Маничев В. Б. 2005 г. ) – анализ разнородных объектов (. NET - Windows) ПА 10 (SADEL-PA 10) – (C, C++, Java, . NET, Open. GL, Direct. X, NVIDIA®
Суть проекта “ПА 10 (SADEL-PA 10)” SC-Super. Computer, PC-Personal. Computer
Программно-методические комплексы для моделирования динамических систем Моделирование электронных схем: SPICE-SCEPTRE: Cadance Design Systems Inc. [www. cadence. com/us/], SPICE: Synopsys Inc. [www. synopsys. com], SPICE: Mentor Graphics Inc. [www. mentor. com] Моделирование механических систем: MSC. ADAMS: MSC Software Corp. [www. mscsoftware. com], EULER: Auto. Mechanics Inc. [www. euler. ru], LMS Imagine. Lab Amesim: LMS International Inc. , [www. lmsintl. com] Моделирование систем управления: МВТУ: МГТУ им. Н. Э. Баумана [http: //model. exponenta. ru/mvtu/20050615. html], MSC. EASY 5: MSC Software Inc. [www. mscsoftware. com] Моделирование разнородных (multi-physics, multi-discipline) динамических систем, описываемых системами ОДУ-ДАУ (прямые аналоги ПА 10 (SADEL-PA 10)): MATLAB-SIMULINK: Math. Soft Inc. , [www. mathworks. com], Maple-Maple. Sim: Maple. Soft Inc. [www. maplesoft. com]
Классическая постановка задач решения систем ОДУ с правой частью 1. Нормальная форма Коши в координатном базисе переменных состояния (явная форма): дифференциальных переменных состояния размерностью m - вектор-функция - независимая правых частей размерностью m переменная (обычно – время) Заданы начальные условия отрезок интегрирования t=[0, TK]. и
Классическая постановка задач решения систем ДАУ с правой частью 2. Дифференциально-алгебраическая форма разных индексов в полном координатном базисе (полуявная форма): - вектор алгебраических переменных размерностью k для полного базиса - вектор-функция размерностью k Заданы согласованные начальные условия:
Классическая постановка задач решения систем ДАУ в общем виде 3. Дифференциально-алгебраическая форма в полном координатном базисе (неявная форма): - вектор-функция размерностью m + k Заданы согласованные начальные условия и отрезок интегрирования. Задача (1) была поставлена в полном координатном базисе Л. Петзольд в 1982 г. и была разработана программа на основе метода Гира DASS – Differential Algebraic Systems Solver:
Постановка задачи решения систем ДАУ в расширенном координатном базисе - вектор переменных состояния размерностью m - вектор производных дифференциальных переменных по времени размерностью m - вектор алгебраических переменных размерностью k - независимая переменная (время) - вектор-функция размерностью m + k Заданы начальные условия и отрезок интегрирования. Найти с гарантированной достоверностью и точностью для всех известных тестовых задач. Задача (2) была поставлена в расширенном координатном базисе: Норенков И. П. , Трудоношин В. А. , Федорук В. Г. Радиотехника, № 9, 1986 г.
Алгоритм решения систем ДАУ в расширенном координатном базисе для неявного метода Эйлера На каждом шаге интегрирования hn для момента времени tn решается система нелинейных алгебраических уравнений (НАУ), размерностью 2 m + k - вектор-функция размерностью m, соответствующая формуле неявного метода интегрирования Эйлера относительно векторов общей размерности 2 m+k методом Ньютона-Рафсона, который сводится к решению соответствующих систем ЛАУ - шаг интегрирования, n – номер шага интегрирования
Схема алгоритма неявных 1 -но стадийных одношаговых методов Эйлера и трапеций (M 1 и M 2 число стадий S=1) tn-1 t 1=tn t hn
Схема алгоритма неявных 2 -х стадийных одношаговых методов (M 4 число стадий S=2) tn-1 t 1=(tn-tn-1)/2 t 2=tn t hn
Схема алгоритма неявных 3 -х стадийных одношаговых методов (M 6 число стадий S=3) tn-1 t 2 t 3=tn t hn
ПЯТЬ БАЗОВЫХ КЛАССОВ ЗАДАЧ РЕШЕНИЯ ЛИНЕЙНЫХ СИСТЕМ ОДУ Уравнение Далквиста (1963 г. ) для исследования устойчивости методов интегрирования - комплексный коэффициент - аналитическое решение - начальные условия Комплексная плоскость устойчивости методов интегрирования III II IV I V - постоянный шаг интегрирования 1. 2. 3. 0 II IV 4. 5. III Фундаментальные решения лежат в основном в области III Фундаментальные решения имеются в области IV Фундаментальные решения имеются в области V
Фундаментальные решения линейных систем ОДУ для классов задач 1, 2, 3, 4, 5 1. 4. 2. 5. 3.
Требование достоверности - основное требование к методам интегрирования Для достоверного решения 5 -ти классов задач Далквиста методы должны быть AL-устойчивыми – абсолютно A-устойчивы строго в левой (Left) полуплоскости комплексной плоскости устойчивости (Евстифеев Ю. А. , Маничев В. Б. 1986 г. ). Область абсолютной устойчивости AL-устойчивых методов интегрирования:
Формулы S-стадийных методов (1988 г. ) интегрирования Маничева-Глазковой Если матрица D единичная, то получим методы Рунге и Кутты, а матрица A и векторы B и C будут соответствовать таблице Д. Бутчера. Наши требования: методы должны быть AL- устойчивыми и должны иметь максимальный порядок точности p
Параметры неявных методов (1988 г. ) M 1(A-), M 2, M 4 и M 6 (AL-) Маничева-Глазковой (совпадают с параметрами методов Лобатто. IIIA)
Кардинальные проблемы реализации ALустойчивых S–стадийных методов интегрирования Маничева-Глазковой 1. Главный недостаток AL-устойчивых методов – наличие «ложных колебаний» (“ringing”), поэтому Хайрер в книгах (1996, 1999) для жестких систем ОДУ-ДАУ назвал такую устойчивость нежелательной (рис. 3. 2. стр. 58) и эти методы не нашли широкого применения на практике. Эта проблема была нами решена в 2006 году. 2. Проблема «сверхточного» (extra precision) решения систем ЛАУ. Эта проблема была нами решена в 2010 -2011 годах.
Алгоритм Маничева-Ильницкого (2006 г. ) реализации AL-устойчивых методов интегрирования без «ложных» колебаний Если в начале или в пределах текущего шага интегрирования по методу М 2 , М 4 , М 6 присутствуют моменты времени, соответствующие разрывам производных интегрируемых функций, то в конце текущего шага выполняется коррекция значений векторов PX и Y путем выполнения корректирующего шага интегрирования, много меньшего текущего шага интегрирования по методу М 1 – неявному методу Эйлера. Вектор X не изменяется!
Проверка алгоритма Маничева-Ильницкого на тестовой задаче Кокина С. А. Математическая модель тестовой задачи Кокина в форме исходной (полученной по законам физики) системы ДАУ вида (2): Аналитическое решение этой задачи для напряжения и тока на емкости С 2 при начальных условиях :
Корректное и точное решение для напряжения uc 2(t), рассчитанного методом М 2 - неявным методом трапеций
«Ложные колебания» для тока ic 2(t), рассчитанного методом М 2 - неявным методом трапеций
Ток ic 2(t) , рассчитанный методом М 2 , по алгоритму Маничева-Ильницкого
Постановка задачи «сверхточного» (extra precision) решения систем ЛАУ A =(aij) – матрица коэффициентов размером n×n, b – вектор-столбец правых частей, размерностью n, х - вектор неизвестных размерностью n. Необходимо вычислить все элементы вектора x с гарантированной точностью в 15 верных значащих цифр (точность double precision на языке Си)
Хронология разработки библиотеки SADEL 1. В 1995 г. была разработана программа-решатель DMAN для решения систем ДАУ в расширенном координатном базисе на языках ФОРТРАН и Си 2. В 2006 г. была решена проблема «ложных» колебаний 3. В 2009 г. было выполнено всестороннее тестирование программы DMAN на универсальном суперкомпьютере СКИФ МГУ ( «Чебышев» ) и программа была зарегистрирована в Роспатенте 4. В 2010 г. была решена проблема «сверхточного» решения систем ЛАУ и разработана 1 -я версия библиотеки SADEL в стандарте языка Си
Сравнение программ-решателей систем ЛАУ ТЕСТОВАЯ ЗАДАЧА ПРОГРАММА-РЕШАТЕЛЬ УИЛКИНСОНА (WILKINSON) 2 ВОХМИНЦЕВА 5 -ГО ПОРЯДКА СИСТЕМ ЛАУ -ГО ПОРЯДКА ОБУСЛОВЛЕННОСТЬ СИСТЕМЫ COND(A)=106 ОБУСЛОВЛЕННОСТЬ СИСТЕМЫ COND(A)=1012 ТЕСТОВАЯ ЗАДАЧА С МАТРИЦЕЙ ГИЛЬБЕРТА (HILBERT) 10 -ГО ПОРЯДКА ОБУСЛОВЛЕННОСТЬ СИСТЕМЫ COND(A)=1. 5*1013 LAE-SOLVER-01 (SADEL) + (15) LAE-SOLVER-02 (SADEL) + (15) MATLAB - (8) - (5) - (4) MATHCAD - (8) - (5) - (4) MATHEMATICA - (8) - (5) - (4) MAPLE - (8) - (5) - (н) NAG-LAPACK - (8) - (5) - (н) IMSL - (8) - (5) - (н) INTEL MKL - (8) - (5) - (н) MAGMA-LAPACK - (8) - (5) - (4) Знак (–) означает, что решатель СЛАУ не обеспечил точность в 15 верных значащих цифр для всех элементов вектора решений теста и не выдал никакого предупреждающего сообщения об этом, в скобках (m) указано полученное количество верных значащих цифр в решении теста, н – нет данных.
Тестовая задача с матрицей Гильберта (Hilbert matrix) 10 -го порядка (все элементы матрицы А и вектора В, как суммы соответствующих строк матрицы А, вычислены с удвоенной точностью) Вектор B – построчная сумма коэффициентов: i=1 Bi= 2. 928968253968250 i=2 Bi= 2. 019877344877340 i=3 Bi= 1. 603210678210670 i=4 Bi= 1. 346800421800420 i=5 Bi= 1. 168228993228990 i=6 Bi= 1. 034895659895660 i=7 Bi= 0. 930728993228990 i=8 Bi= 0. 846695379783620 i=9 Bi= 0. 777250935339180 i=10 Bi= 0. 718771403175440 Абсолютно точное решение этой тестовой задачи: Единичный вектор
Результаты решения с удвоенной точностью на языке Си тестовой задачи с матрицей Гильберта (Hilbert matrix) 10 -го порядка Библиотека MAGMA (LAPACKLinpack): 1. 00000464697 0. 999999960244827 1. 000000840984605 0. 999992393461859 1. 000036137385782 0. 999900979283788 1. 000162025520190 0. 999843779478675 1. 000081852296234 0. 999982030574678 Библиотека SADEL: i=1 Xi= i=2 Xi= i=3 Xi= i=4 Xi= i=5 Xi= i=6 Xi= i=7 Xi= i=8 Xi= i=9 Xi= i=10 Xi= 1. 000000000000000 1. 000000000000000
Сравнение программ-решателей систем ОДУ-ДАУ ТЕСТЫ Программа- ТЕСТ 1 ТЕСТ 2 ТЕСТ 3 ТЕСТ 4 решатель Жесткость Коэффициенты Жесткость Параметры MU=1 e 6 MU=1 e 9 kt=1, ki=1, ku=1, ki=1, MU=1 e 6 реального ku=0. 01 kt=1 e-104 систем ОДУ-ДАУ SADEL C-Library М 2, М 3, 2010 Mathcad Radau, 2007 MATLAB Ode 15 s, 2007 Maple Rosenbroсk, 2007 Mathematica BDF, 2010 лазера + + + + - - - + + + - - - - - + + Знак (-) означает недостоверное решение теста без предупреждения
Сравнение программ-решателей систем ОДУ-ДАУ ТЕСТЫ Программарешатель систем ОДУ-ДАУ ТЕСТ 1 ь MU=1 e 9 MU=1 e 6 SADEL C-Library М 2, М 3, 2010 NAG C-Library 2012 Метод BDF IMSL C-Library 2012 Метод BDF ТЕСТ 2 ТЕСТ 3 ТЕСТ 4 Коэффициенты Жесткость Параметры kt=1, ki=1, ku=1, ki=1, MU=1 e 6 реального ku=0. 01 Жесткость ТЕСТ 2 kt=1 e-104 лазера + + + + - - - Знак (-) означает недостоверное решение теста без предупреждения
ТЕСТ 1 Пример расчета жесткой системы ОДУ 2 -го порядка (тест Ван дер Поля) MU – параметр жесткости. Сравнение было проведено для часто встречающихся на практике параметров жесткости MU=106 и MU=109
Описания, НЕКОРРЕКТНЫЕ РЕШЕНИЯ ТЕСТА 1 в MATHCAD, MATLAB, MATHEMATICA, MAPLE и КОРРЕКТНОЕ РЕШЕНИЕ в SADEL
ТЕСТ 2. ЛИНЕЙНАЯ ЭЛЕКТРИЧЕСКАЯ СХЕМА С ИЗВЕСТНЫМ АНАЛИТИЧЕСКИМ МНОГОПЕРИОДНЫМ РЕШЕНИЕМ (тест Маничева) RLC: high Q filter
Система ОДУ для ТЕСТА 2
Описания, НЕКОРРЕКТНЫЕ РЕШЕНИЯ ТЕСТА 2 в MATHCAD, MATLAB, C-Library NAG и КОРРЕКТНОЕ РЕШЕНИЕ в Си библиотеке SADEL
ТЕСТ 3 Нелинейная жесткая система ОДУ, имеющая локально-неустойчивое решение (тест Скворцова). Сравнение программ проведено для часто встречающихся на практике значений жесткости MU=106.
Описания, НЕКОРРЕКТНЫЕ РЕШЕНИЯ ТЕСТА 3 в MATLAB, Maple, C-Library NAG и КОРРЕКТНОЕ РЕШЕНИЕ в SADEL
ТЕСТ 4 Нелинейная жесткая система ОДУ для моделирования процессов реального лазера (тест Евстифеева). Сравнение программ проведено для параметров реального работающего лазера
Описание, НЕКОРРЕКТНОЕ РЕШЕНИЕ ТЕСТА 4 в Maple и C -Library NAG, КОРРЕКТНОЕ РЕШЕНИЕ в SADEL
Задача химической кинетики – реакция пиролиз (описание задачи в пакете Mathematica 8. 0) • • • s=NDSolve[{x 1'[t]Š-k 1*x 1[t]-k 2*x 1[t]*x 2[t]-k 4*x 1[t]*x 6[t]+k 8*x 4[t]k 14*x 1[t]x 12[t], x 2'[t]Š 2 k 1*x 1[t]-k 2*x 1[t]x 2[t]-k 6*x 2[t]x 5[t]+k 7*x 8[t]k 10*x 2[t]x 5[t]k 11*x 2[t]x 9[t]+k 14*x 1[t]x 12[t], x 3'[t]Šk 2*x 1[t]x 2[t]+k 10*x 2[t]x 5[t]+k 11*x 2[t]x 9[t], x 4'[t]Šk 2*x 1[t]x 2[t]-k 3*x 4[t]+k 4*x 1[t]x 6[t]+k 5*x 5[t]x 6[t]2 k 8*x 4[t]+k 9*x 5[t]x 8[t], x 5'[t]Šk 3*x 4[t]-k 5*x 5[t]x 6[t]k 6*x 2[t]x 5[t]+k 7*x 8[t]+k 8*x 4[t]-k 9*x 5[t]x 8[t]-k 10*x 2[t]x 5[t]k 13*x 5[t], x 6'[t]Šk 3*x 4[t]-k 4*x 1[t]x 6[t]-k 5*x 5[t]x 6[t]k 12*x 9[t]x 6[t], x 7'[t]Šk 4*x 1[t]x 6[t]+k 12*x 9[t]x 6[t], x 8'[t]Šk 6*x 2[t]x 5[t]-k 7*x 8[t] -k 9*x 5[t]x 8[t]+k 14*x 1[t]x 12[t], x 9'[t]Šk 10*x 2[t]x 5[t]-k 11*x 2[t]x 9[t]k 12*x 9[t]x 6[t], x 10'[t]Šk 11*x 2[t]x 9[t]+k 12*x 9[t]x 6[t], x 11'[t]Šk 9*x 5[t]x 8[t], x 12 '[t]Šk 13*x 5[t]k 14*x 1[t]x 12[t], x 1[0]Š 0. 040746269, x 2[0]Š 0. 0, x 3[0]Š 0. 0, x 4[0]Š 0. 0, x 5[0]Š 0. 003955224, x 6[0]Š 0. 0, x 7[0]Š 0. 0, x 8[0]Š 0. 0, x 9[0]Š 0. 0, x 10[0]Š 0. 0, x 11[0] Š 0. 0, x 12[0]Š 0. 0}, {x 1, x 2, x 3, x 4, x 5, x 6, x 7, x 8, x 9, x 10, x 11, x 12}, {t, 0. 0, 14. 55}, Method "Stiffness. Switching"] Plot[Evaluate[{x 8[t]}/. s], {t, 0, 14. 55}, Plot. Style->Automatic] x 8[14. 55]/. s
Задача химической кинетики – реакция пиролиз (неверный график для переменной х8 в пакете Mathematica 8. 0)
Задача химической кинетики – реакция пиролиз (верный график для переменной х8 в SADEL )
Задача химической кинетики – реакция пиролиз (исправленное описание и верный график для переменной х8 в пакете Mathematica 8. 0)
Скачать презентацию Всероссийская научно-техническая конференция Проблемы разработки перспективных микро- и
b52011c6758c33a82ddf47f86a1c3020.ppt