Введение к работе
, !
- '--'-дНуальность теми, Математическое моделирование динами .о -ских процессов находит применение как при научных исследованиях и проектировании, так и в процессах управления. С разработкой и совершенствованием высокопроизводит ьных ЭВМ, на современном этапе развития науки и —хники наблюдается появление новых типов задач и постоянное усложнение их математических моделей. Эти модели описываются сложными нелинейными системами обыкновенных дифференциальных уравнений или уравнениями в частных производных.
Многие современные численные методы решения краевых задач, в конечном счете, приводят к необходимости решения последовательности задач Коши для обыкновенных дифференциальных уравнений. Тем самым, а исследовании стационарных и динамических систем важное место занимают численные методы решения задачи Коши для систем обыкновенных дифференциальных уравнений Возникает вопрос, насколько полно ресена проблема численного исследования таких задач?
Еще сравнительно недавно,. несколько десятилетий назад, было принято считать, что, за исключением некоторых второстепенных вопросов, существующие на то время численные методы обеспечивали достаточную эффективность реюния задачи Коши для обыкновенных дифференциальных уравнений. Однако, уже в работах Кертиса и Гиршфельдера, а затем Далквиста- было начато обсуждение проблемы жесткости систем дифференциальных уравнений (СДУ), с которой существующие численные методы на могли справиться.
. С тох пор, вопросам разработки теории численных методоь, пригодных для решения жестких систем обыкновенных дифференциальных уравнений, уделяется усиленное внимание.
Наиболее интенсивно, развивались неявны численные методы. Были построены и теоретически.обоснованы различные типы таких методов как одношаговых, так и многошаговых, например, неярнью методы типа Рунге-Кутты, многошаговые кеп да Гира и Энрай э, -методы типа Розенброка и др. Все зги методы в процессе их реализации при решении нелинейных СДУ требуют многократного решения систем нелинейных алгебраических уравнений, что является одним ' из основных ограничений применения неявных методов Рунге-Кутты. При решении системы N ди$ферени»альньк уравнений с использованием S- этапного метода в общем случае необходимо решать нелинейную систему из N-S алгебраических уравнений .с большим
количесгвом вычислений правой части СДУ. Хотя теоретически, при определенных условиях, неявные методи Рунге-Кутты обладают свойством нелинейной устойчивости, при их практической реализации это свойство может не сохраняться, что требует соответствующих ограничений на шаг интегрирования в каждом отдельном случае. Ввиду приведенных замечаний, эти методы являются в ряде важных для практики случаев малоэффективными .
Более экономными, с точки зрения вычислительных затрат, являются методы типа Розенброка. К некоторой разновидности таких методов следует отнести и методы Е.Н.Новикова. Все эти методы не являются полностью неявными и на обладают таким же сильным свойством нелинейной устойчивости, как неявные методы. Построены методы только ди четвертого порядка согласованности, обладающие свойством А-устойчивости. Вывод формул А- устойчивых методов вызоких порядков довольно сложный, а их реализация требует существенного увеличения количества вычислений правой части СДУ. Кроме того, на каждом шаге интегрирования необходимо проводить вычисление матрицы Якоби правой части СДУ. Тем самым, значительные вычислительные затраты и невысокий порядок точности результатов ставит под сомнение целесообразность использования методов типа Розенброка для исследования жестких задач.
Наиболее широкое распространение получили неявные многошаговые методы' Гира. Они .как и другие многошаговые численные методы, наиболее экономны по количеству вычислений правой части системы уравнений. Для них существуют довольно простые оценки локальной погрешности. Однако, только методы не выие второго порядка согласованности обладают свойством А-устойчивости. Для методов более вызоких порядков введены понятия А(а)- устойчивости , жесткой устойчивости. , Ввиду специфики области устойчивости, многошаговые численные методы выше второго порядка не обепечиаают вычислительной устойчивости для систем с осциллирущими решениями. Кроме того, характеристические уравнения всех известных линейных многошаговых методов имеют, так называемые, "паразитические" корни, которые существенно влияет на область устойчивости и точность приближения решения, а тем самым, на величину шага интегрирования. Совр^мочнью многошаговые численные методы непосредственно предусматривают численное интегрирование систем с постоянным тагом. Дія повывония эффективности, в программном обеспечении катодов предусматривают определенные алгоритмы изменения
шага интегрирования. Эти алгоритмы обеспечивают приближенное определение ( метод, ::нгерполяц!!и . вектор Нордсика и др. ) дополнительных значений решения, и его производной. Они требует увеличения вычислительных затрат, а погрешности их результатов увеличивают погрешность искомого решения. Приведенные вше рассуждения позволяют прийти к вызоду, -.го и неявні» линейные многошаговые численные методі ,е являются достаточно эффективными.
Как указывалось выше, любые неявные численные методы в процессе реалізації требуют решения систем не; иейных алгебраических уравнений, что в свою очередь порождает проблему существования единственности ревения этих систем. Для рэпения нелинейных алгебраических систем наиболее широко используется итерационньй процесс Ньютона. Итерационный процесс, в свою очередь, требует сследования сходимости последовательности итераций и выбора начального приближения.
Все эти замечания дополнительно подтверждает, что линейны*, неявные численные методы не решают лрс. лемы жесткости дифференциальных уравнений.
Недостаточная эффективность неявных численных методоз привела к поиску новых направлений в построении численных методов решения жестких систем дифференциал :ых уравнений. С этой цель» в работах Вам5ека, Ламберта, В.В. Бобкова, П.И. Бод-нарчука, Ю.В. Ракитокого и др. были использованы дробно-рациональный преобразования Падэ. для построения, так -называемых, дробно-рациональных численных методов. В' этом направлении приводились и исследования автора. Эти методы для скалярного дифференциального уравнения при опроделент'ч условиях обладают А- или L- устойчивостью. Однако, при непосредственном распространении их на системы уравнений возникали непреодолимые сложности. Эти сложности пытались преодолеть путем раздельной .дискретизации по отдельным компонентам векторного' решения ' системы, т.е. путем, т.н., покомпонентной . реализации. Такой подход был очень привлекательным, ибо методы становились явными и их реализация требовала незначительного количества вычислительных затрат. Но в работах Соттаса было показано, v о 'методы Вамбека при покомпонентной реализации неустойчивы на жестких системах. В дальнейшем оказалось, что при покомпонентной ' реализации и другие 'дробно-рациональные мет^ы не обладают устойчивость».
Методы Ю.В.Ракитокого и В.В.Бобкова предназначены для ра-
шения линейных систем с постоянными коэффициентами. Они базируются на аппроксимации матричной экспоненты многочленами по степеням матрицы системы. Более перспективными и эффективными являются методы В.В.Бобкова, обладающие свойотвом адаптивности, что даот возможность гибко изменять шаг на различных этапах интегрирования жесткой системы.
Походя из приведенного выше краткого анализа современного состояния теории численных методов реиения жестких систем обыкновенных дифференциальных уравнений, следует, что ЕОПРОСЫ построения и исследования новых более эффективньк численных методов являются актуальными.
Целью диссертационной работы является: создание и развитие основ теории новых видов дробно-рациональных приближений решения жестких систем различных классов обыкновенных дифференциальных- уравкэний, обладающих требуемым типом устойчивости; построение и обоснование одношаговых и многошаговых дробно- рациональных численных методов; реализация разработанных алгоритмов в виде программных модулей; исследование эффективности предлагаемых численных методов.
Методика исследований. Для достижениия цели использовались общая теория алгебраических и дифференциальных уравнений, теория современных численных методов.
Научная ноэизна.; В процессе исследований разработаны основы общей теории новых эффективных дробно-рациональных безите-рационных численных методов решения задачи Коши для жестких систем обыкновенных дифференциальных уравнений. В частности:
- введена новая структура дробно-рациональных приближений (ДРП) решения задачи Копи для ОДУ в виде усреднения с дробно -рациональными весовыми коэффициентами тейлоровских приближений последовательных порядков точности; исследованы условия обеспечения порядка аппроксимации и устойчивости приближений, показано отсутствие "барьера Далквиста",т.е. возможность построения А-или I-устойчмвых ДРП любого конечного порядка согласованности;
доказано, что аппроксимации Падэ являются частными случаями ДРП показательной функции;
построены и обоснованы по точности и устойчивости дробно-рациональные численные методы решения линейных систем с постоянными коэффициентами, изучены вопросы их распространения на нолинейны системы с применением процесса кзазилинеаризации,
определены апостериорные оценки точности приближений дробно-рациональными методами линейных систем относительно нелине .ной системы;
- для жестких СДУ высших порядков впервые предложена
методика построения' приближений решений без приведения к
системам первого порядка, введе. , понятие и определены
необходимые и достаточные .ловия Ао-устойчивости приближений
для систем с осциллирующими решениями;
изложены и обоснованы новь» принципы пог-роения приближений режний неявных систем ОДУ первого порядка на базе их преобразования к явным системам второго порядка, что, в отличие от существующих методов» позволяет избежать необходимости решения систем алгебраических уравнений;
- разработана "-юрия устойчивых однопаговых дробно-рацио
нальных численных методов на базе новой универсальной,
отличающейся от общепринятой, методики построения вложенных
линейных однопаговых методов произвольного п ядка; с целью
уменьшения вычислительных затрат на пат интегрирования рас
смотрена структура блочных однош^'овых методов;
7 излагается теория многошаговых дробно-рациональных численных методов о переменны* шагом интегрировз' 'я, характеристические уравнения которых, в отличие от других известных многошаговых методов, не содержат посторонних корней; исследованы вопросы устойчивости и доказано, что на устойчивость не влияет переменный шаг интегрирования и многошаговость методов', определены оценки локальных гогрешостей методов;
- построены алгоритмы реализации численных методов о
выбором оптимального порядка и величины шага интегрирования,
разработаны новые тестовые задачи, на базе которых проведены
экспериментальные исследования качественных- и количественных
характеристик методов. * . .
Практическая, ценность. Полученные результаты могут быть
использованы при исследовании реальных динамических процессов,
математические модели которых описывается жесткими системами
дифференциальных уравнений. Программные модули реализации
алгоритмов включены в библиотеки кате атического. обеспечения
ЭВМ. ".". - - ' - :
Дппробация роботы. Результаты работы докладывались на: ежегодных научно-технических і ;Н$еренциях Львовского политехнического института.(1974-1990гг.); Объединенном семинаре по
матетатике и шханике Западного . научного центра АН УССР (1978- 1990гг.); республиканском научном семинара "Проблемы вычислительной математики" (г.Львов,1978 - 1990гг.); научных семинарах кафедры вычислительной математики Московского госуниверситета (1985г.), кафедры численных методов математической физики Киевского госуниверситета (1987,1990гг.); Всесоюзных школах-семинарах академика А.А.Самарского "Теоретические и прикладные проблемы-вычислительной математики" (г.Дрогобыч,1980г.; г.Рига,1982г.; г.Львов.1983г; г.Немала, 1985г.; Шусенокое, 1986г.; г.Одесса, 1987г.;г.Светлогорск,1988г.; г.Владивосток, 1989г.); Всесоюзном семинаре "Вопросы оптимизации вычислений" (г. Алушта, 1987г.); Всесоюзных конференциях "Новые подходы к ре-Еению дифференциальных уравнений" (г.Дрогобич, 1987, 1991гг.); Всесоюзном семінаре "Численные методы решения жестких систем" (г.Красноярск,1988г.); Боесоюзных школах-семинарах "Распараллеливание обработки информации(г.Львов,1987,1989гг.); Международной конференции "Математическое моделирование и прикладная математака"(г. Москва, 1990г.), Республиканской конференции "Математическое моделирование элементов и фрагментов БИС" (г.Рига, 1990г.)Международной конференции. "Информационные технологии в анализе изображений и распознавании образов"(г.Львов,1990г.)
Публикации. Го результатам диссертации опубликовано Є5 работ, основные из которых вклхнены в список литературы автореферата.
Структура работы.Работа состоит из введения, трех глав.за--клкчения.списка литературы из 166 наименований и раздела приложений.
СОДЕРЖАНИЕ РАБОта.
Во введении обоснована актуальность избранной темы, сформулированы цели диссертационной работы, описана методика исследований, определена научная новизна результатов, кратко изложены основные результаты работы.
Глага 1 "Дробно-рациональные приближения решений систем дифференциальных уравнений" посвяцена теории дробно-рациональных приближений (ДРП) решений различных классов жестких систем дифференциальных уравнений. Главное внимание уделяется изучению свойств ДРП, их устойчивости.
В 51.1 рассмотрены общие принципы построения ДРП решения задачи Кош для системы уравнений
Y' = f(x,Y). Y(x0) = Y0, (1)
Y R3, Г:Н8+1- R3, x « R,.
в виде g QhU .
Y[p] . J=o 3 p'J'n : .2)
где Y^ - приближение решения задачи (1) в сеточном узл ^+1 = xn+h, Тк - - тейлоровское приближение к-го порядка решения в том же узле относительно узла Гд, Cj (1=1 ..р)-. матричные коэффициенты, независимые от h, CQ=1 - единичная матрица. Под операцией деления в (25 следует понимать уіиножение обратной матрицы знаменател злева'на вектор числителя.
Формула (2) обладает свойством усреднения с дрсбно-рациона-дьнымп весовыми коэффициентами
С h
Е С<ПД 1=0 1
последовательности тейлоровских приближений.
Определение 1. Будем говорить, что приближение вида (2) является согласованна о порядком р, если *У(х)еСр+ ^^1 имеет несто' соотношение
Теорема 1. Пусть С0-1, тогда приближение (2) является сог-ласованнш с порядком р независимо от выбора значений коэффицг ентов С^ (1=1..р).
Введено понятие ошибки приближения, возникающей за счет ограничения порядка согласованности, которая определяется выражением
' У<Р-^1>
h Jo0;! тр=жтт
R0r p
E (Mr
1=0 i
В 1.2 обуславливается вьбор матричных коэффициентов С., в
виде . .
' Ct = (-1 iVjJ*. а=0..р). (3)
где Jn *= а ffi r=3n" значение матрицы Якобн правой
и части системы (1) в сеточном узле Хд, ttj,- скалярные параметри
(0^=1) и исследуются вопросы уотойчивооги ДРП , представленньк с учетом (3), соотношением
X (-1)V,J;| Ь*!^., „
„fnl i-O J п р-о»в
ylPJ „ Д=У (4)
i=0 in
на модельной линейной система
Г + АУ = О. (5)
где А - положительно определенная патрица о постоянными элементами.
С цель» упрощения процесса исследований показана правомочность замены системы (5) скалярным уравнением
Y' = ->.Y . . (6)
с комплексными значениями параметра X.
Приближение экспоненциального родения уравнения (6) формулой (4) определяется соотноиенияш
*J-V0)V. --(7)
3=0 1=0 1;> 1'1 х
В (а) = , о - Ал. (8)
Е of1
Теорема 2. Яхбой («С.ш> злзкенг таблицы Падэ содержится в 18).
Тем самим аппрокцикации Падэ являются, только частными случаями ДРП (8) для экспоненциальной функции.
Согласно известным определениям, на' устойчивость формулы (4) влияют параметры a.j (3=1..р) в соотношении (8).
Т<орема 3. Если в (4) значения параметров определяется соотношением aj=1/U. (1=1..р>, то операторная функция (8) обладает свойством Lq- допустимости.
Этот параграф состоит из трех подпунктов.
В 1.2.1 исследуются условия жесткой или А(а) устойчивости ДРП. Определена общая формула таких приближений произвольного конечного порядка *-.
-1T-
2H! TTJnrHn
In+1 p j hi ^
Приведены'графики областей уотоЯ'"івости для р = 1..7, анализ которьк показал, что пр'ближения первого и второго порядков согласованности являются 1 - устойчивыми, а приближения вьюших порядков обладают А (а).устойчивостью.
В 1.2.2 для обеспечения А- устойчивости ДРП (4) произволь
ного р-го порядка согласованности определены значения пара
метров ' '
' «і-<$%№.'. (l-l;.P).
где Ср - число сочетаний из р элементов* по І. _. Приведены примеры А-устойчивьк ДРП. от первого до восьмого порядков. Показана возиокность построения прчб." -кений различных порядков с общей операторной функцией устойчивости.
В 1.2.3 изложена- методика пос.роения 1-устойчивых ДРП. Определена общая формула для значений параметров
ч-оіЩІЇ1- <1=1..р>.
ДРП произвольного р-го порядка согласованности, при которьк обеспечивается L-усгойчивосгь приближений вида (4), Показано существование L-устойчивьк ДРП различных порядков с общей операторной функцией устойчивости. Приведены примеры 6 групп таких приближений от первого до седьмого порядков согласованности Построены и исследозаны графики их областей устойчивости.
В 1 -3 рассмотрены вопросы построения дробно-рациональных численных методов решения линейных систем с постоянными коэффициентами и их распространения на нелинейные системы с использованием процесса квазилинеаризации.' Он состоит .из 5. подпунктов.
В 1.3.1 на базе ДРП вида (4) поостроены дробно-рациональ-
. ные методы (ДРМ) решения задачи Копи для линейных однородных
систем (5). Метод произвольного р-го порядка описывается соот
ношением .
J=0 i=0
JT-^i-i
*i = -—p—п V (9)
На оснований общей методики, в качестве примера, определены значения параметров а± для ыэтодов до восьмого порядка, обеспечивающие 1-устойчивость для штодов нечетных порядков и А-устойчивооть для мэтодоч четных порядков.
Рассмотрено сходство ДРй рааекия линейных систем о системными ыетодзла Ю.В.Ракитского.
В 1.3.2 определены практические оценки точности ДРй репения
ЛИНеЙНЫХ СИСТеМ . n.-U-l
iya). р—^ Гп - (Ю)
Е аЛрк3 3=0 3
Исходя из значении параметров си, для катода заданного порядка р формула (І0) позволяет оценивать, величину локальной погрешности мэтода. Такие оценки даны для катодов, описанных в предыдущем подпункте.
В 1.3.3 приведена нетодика построения численных ДРМ решения задачи Кош для линейных неоднородньк систем первого порядка о постоянны*! коэффициентами
Г 4- AY » В , У(0) » Y0. (U)
где А - постоянная положительно определенная матрица размерности sxs, а В - s-мерныЯ постоянна вектор.
Для А(а)-устойчивых ДРМ определены связи мезду числителем и знаменателем дробно-рационального выражения методов
Р h* k „ где S * E -;]tH77 A " единичная матрица.
Приведены примеры А- и L-усгойчивых методов до восьмого порядка.
В 1.3.4 описан процесс преобразования с помощью квазилинеаризации задачи Кош для нелинейкой автономой системы
У = ГШ, Y(Iq) = Y0 (12)
к итерационной последовательности линейных систем
Vi(V = Yn« m = 0.1.2,... ' (13)
Yn+i и " n"T0G приближение решения исходной задачи в сеточном узле\+1.
В 1.3.5 исследовано влияние неточности решений линейных систем (13) в отдельных итерациях на сходимость итерационного процесса к реиению нелинейной системы, исходя из квадратичной сходимости точных реиений їд+1 п задачи (13) к решению Уп+1 нелинейной задачи (і?1 в coтoчнo^s узле з^+1, т.е.
где q « (OJ) - константа, не зависящая от п.
Пусть Yn+1 - приближенное ресоние в узле Jp+1 линейной задачи (13) на п-том отапе нтерс .лскного процесса, полученное некоторый ДРМ из семейства, приведенного в продвдудем подпункте. '
Тогда справедливо неравенство
. IVi|B1-Vl.nHa/ Vffl>i. (15)
где as JRpt^+^S. р- порядок пригоняемого ДРМ.
Пусть при реализации итерационного процесса, для некоторого пей вьтолняется условие
где 0 - достаточно малое заданной число. В процесое практической реализации последовательных итераций необходимо согласовать вьбор величин а и. Оп для достижения заданной точности в приближенного роптания Yn+1 m нелинейной системы, т.е. выполнение неравенства
От" ? W"*n+1 ІІБ' аП (,7)"
Теорема 4. Если для- итерационного процесса квазилинеари
зации (13) выполняются условия (14),(15)- hOHS 1/{9q), то
для достаточно большого ш справедлива оценка приближения реше
ния . ' ;
-u-
Из теоремы 4 следует, что в процессе квазилинеаризации при выполнении условий теоремы, влияние погрешности приближенных решений линейных систем на окончательный результат при больших m асиютготически не превышает удвоенной погрешности численного метода.
Используя результаты теоремы 4, получим Дщ ^ 4а.
Теорема 5. Если выполняется условие а 1/(24q), то при больших m для велччины й^ ассимптотически справедливы оценки
2
' а?п+1.т.1 -Vi.ral-Z (е-За). которое является критерием окончания итерационного процесса.
В 1.4 рассмотрены сложности реализации ДТП решений сильно -жестких систем. Реализация ДРП общей структуры (4) требует решения системы линейных алгебраических уравнений с матрицей системы в виде полинома от матрицы Якоби правой части системы (1). Кроме значительных вычислительных затрат при возведении матрицы Якоби в степени, соответственно отдельным членам полинома, такая полиномиальная матрица в случае жестких систем становится плохо обусловленной. Для преодоления возникающих' затруднений предложено специальную структуру ДРП в виде
X. (-1 )3^^ « w
, t = Д= _ „—, (їв)
(Е - арЫг/
где а - параметр приближения р-го порядка согласованности. Приближение М8) получаем из (4) при значениях о^ = Срар*
Реализация ДРП (18) требует решения р линейных алгебраических систем с общей матрицей системы С^ = Е - a fcJn. Если предварительно провести декомпозицию матрицы С^, то решение указанных р систем требует незначительных дополнительных вычислительных затрат.
Исследованы условия устойчивости ДРП вида (18), располагая, в данном случае, только одним параметром а Показано, что значения а для L-устойчивьк ДРП определяется уравнением
а для -устойчивых p-ro порядкаДРП из уравнений:
а) для четного порядка
^о<-1)іТр^ЬтСРаГ=* :
б) для нечетного ПО РЯД/ -.
Приведены значения параметра а для Ь- и А-устойчивых ДРП до шестого порядка.
С целью уменьшения . вычислительных . затрат при вычислении числителя соотношения (18), который содержит произведения степеней матрицы Як.би на вектор тейлоровских приближений, проведено разложение дроби (18) на сумму элементарных дробей, т.е.
-YtP' ,-Y Л *> * "*' (19)
1.1 № - aphJn,
Реализации'вычислений формулы 09) для "ДРП р-го порядка целесообразно приводить по схеме Горнера. которая представляется рекуррентными соотношениями линейных систем с общей матрицей С^
ABi-cJ"Vi-i + Bi-i'
' . во = « dj \Ei(-1)k^.n- -1 -X+ ВР
Определение правых частей в линейных, системах рекуррентного процесса (20) сводится к операциям умножения векторов на скалярные величины и их суммирования. .
Приведенные преобразования не влияют' на порядок согласованности и устойчивость ДРП.
В 1.5 излагается методика построения н& базе формулы О-Л
"ДРЛ решения задачи Коки для жестких систем дифференциальных
уравнений вьна первого порядка
Y{m) = тд-Д\...Д(ш~1)) . ' - (21).
с- начальными условиями
Y(x0)=Y0, S'(V«Y6«— Y(m"1)U0)=Y^m~1). rfleY^R3
без преобразования к системам первого порядка.
При преобразовании системы (21) к системе первого порядка и реализации вычислений рекуррентній процессом (20) возникала бы необходимость решения линейных систем размерности ВхШ.
Предлагаемые ДРП в процессе их реализации оперируют системами линейных алгебраических уравнений только размерности s, которая совпадает с размерностью исходной систеглы (21) и не зависит от порядхз системы.
1.5 состоит из трех подпунктов.
В 1.5.1 излагаются общи идеи построения ДРП для системы вида (21) без приведения к системам первого порядка на примере системы второго порядка
Г » f(Y;r). Y(z0;~Y0. Г (s0)=y6- Y* R- (22) При построении ДРП использован особьй случай матрицы Якоби расширенной системы первого порядка для системы (22)
f О Е І
где Гу»*Г(Уп.У)/*Г, Гу,=*Г(УпДД)/лУ*. которая применяется в процессе реализации рекуррентной процедуры (20) для систеш
Г Г = Z
\ (24)
I V = J(Y.Z)
с раздельный определением приближений Y^] и Y^+? n+i в св" точном узле X^f
В результате рекуррентный процесс определения приближений р-го порядка согласованности решния системы (22) описывается соотношениями:
Vie'" Di-1,2 + W-I.Y ' B«" Di-1.Y + V^is * (i=1"P)
yn?1=Yn+V' у,Й1"уА4Л«- (25)
гда Di* (й- h' (* ] ^"E"aphfY?" Ф2**
Для реализации процесса необходимо один раз на ваге ин-
тегрнрования провести декомпозиции матрицы С^ размерности з, а затем последовательно проводить р раз рекуррентные вычисление
В ..5.2 исследуется вопросы устойчивости ДРП вида (4) решений систем второго порядка на («дельном уравнении
Y" = -u2Y. ' (26)
Посла преобразования уравж ія (26) к виду
У = -AY. А - [ 2 ~1 } . (27)
операторная функция перехода (8) решения системы (27) принимает вид
D (h) = _Еі где р . н)Іе М )1-^ .
Е+ Е аЛі3А3 3~и
j=l >
Определенна 2. Еслл собственные значения к рицы D(п) при заданном h лежат на единичной окружности в комплексной плоскости, то ДРП решения называется и. - устойчивьм для указанного значения h.
Определение 3. Если ДРП являются и-устойчи^ыми при всех зна- чениях h>0, то такио приближения будем называть Аш-устойчивъми.
Свойство Aw - устойчивости представляет особый интерес при исследовании осциллирукщих задач.
Теорема 6. Для того, чтобы ДРП вида (4) обладали свойством Ш - устойчивости, необходимо и достаточно, чтобы кезду коэффициентами операторной функция D (h) существовала связь
Pi = (-1)1. (і - 1-Р)« Из теоремы следует, что А-усгойчивш ДРП р-го порядка является Аы-устойчивкми тогда и только тогда, если их операторная функция D(h) представила в виде *
Vh) = c(tt) ГД Vh) " Е + ? aih А
Среди ДРП вида (25), базирующихся на ДРП (18), только А-уо-
тойчивые приближения пэрвого и второго г^рядков согласованности
является Аш - устойчивыми.' '
-18-.
' В 1.5.3 методика построения ДРП для систем второго порядка обобщается на системы третьего порядка
У" = t(Y,Y',Y"), Y« Rs. (28)
1() = С,. Y'{xQ) = С2. У"(ї0) = 0Э.
В результате описан рекуррентный процесс реализации приближений, который требует решения линейных систем алгебраических уравнений размер: ;ооти s, совпадающей о размерностью исходной системы (28).
Приведенная мэтодика может быть использована для построения ДРП решения ситемы диїференциальньк уравнений произвольного m -го- порядка вида (21).
В 1.6 рассмотрены принципы построения ДРП решений отдельного класса неявных систем дифференциальньк уравнений первого порядка.
F(X.Y.Y') =0, YU0) - ї0, ' (29)
где х « R, ї є Б3-, Y* * Re, F(x,Y,Y') -непрерывная и дифференцируемая по воем аргументам вещественная е-- мерная вектор-функция и предполагается существование неособенной положительно определенной матрицы Fy.^F/tfY', а также ограниченность производной оУ/<Я.
Пусть, кроме торо, известно значение Y'U0) * X'Q, которое является корнем векторного уравнения ?(x0,Y0,Y') = 0.
Для построения ДРП система (29) была преобразована к виду.
Y- - -Py1(PyY« + Рх), Y(x0)-Y0. Y'(X0)=Y, (ЗО)
а затем была использована методика построения ДРП решения систем второго порядка раздельно для Y^?J и 1^ описанная в предыдущем параграфе.
В результате построены и обоснованы ДРП первого порядка
hPr Zn ' h(?y Zn+?_ )
n n n n
и второго порядка *І?і ' V^1FY2hZn-('(hFYVr<FYnV\)))
41l'V%*4Bvv vw (32)
n її n n
ОТНОСИТеЛЬНО решения Yn+1 И Yn+1=Z^Il+1
Проведены исследования устойчивости ДРП, в результате
которых показано, что при а.=1 приближения (31) являются L -
-і устойчивый, а при а.,= g - является А-устойчивыми.
Доказана
Теорема 7. Если параметр а,= g, то ДРП (32) обладает А-ус-
тойчивостью о общей операторной функцией
Е-^А
щи) = — .
Е + А
Ввиду сложности непосредст: иного вычисления производных высших порядков искомого решения Y(x), приближения более высоких порядков решения задачи (30) рекомендуется определять многошаговыми методами.
В 1-Т приведены дробно-рациональные численные методы решения линейных систем дифференциальных. уравнений с постоянными коэффициентами первого и произвольного k-го порядков, построенные на базе преобразованных ДРП по методико, описанной в 1.4. Они нэ требуют громоздких операций умножения матрицы на матрицу и матрицы на вектор. Параграф состоит из трех подпунктов.
В 1.7.1 изложен, методика построения методов решения задачи Коси для линейных систем о постоянными коэффициентами
У = AY + В, У(30) - Y0, їеН". A=RNzN, BsRN. (33)
В результате построений метод р-го порядка .описывается соотношениями
«p.,»- ^v1- 4Р)- <-1)P+k So »:*-* ^ (к=1 'р)
Определена оценка локальной погрешности метода произвольного р-го порядка точности
і »
где 7р(а,,)- |Дм)«7р^зттОІаЗ|.
^(^,) - точное решение задачи (33) в сеточном узле їц^. (| - некоторая норма в іг. В 1.7.2 построен алгоритм реализации методов из подпункта 1.7.1 о автомагическим изменением порядка метода и ката
интегриро' "іния. Для построения алгоритма изменения шага первоначально используются упрощенные оценки погрешностей
Rk = v4(ap) J іААіУ^'ц . (к>р-1 .p.p+1) (36)
Если Rp
При выполнении условия Вр<є, новьй шаг интегрирования определяется соотношениями
\юв - * . ' 1 (37)
где т = maxCtp^.ip.tp^), \ = (є/іу^. (k=p-1.p,p+1).
Порядок, метода на пооледущем шаге интегрирования соответствует порядку, определяющему максимальное значение коэффициента Т.
Для проведения экспериментальных исследований предлагаемых численных методов разработан программный модуль DM1S1 на алгоритмическом языке Бейсик. Эффективность методов проверяется нз трех тестовых задачах с наборами различных ' вариантов исходных данных , которые в достаточно полг ' мере отображают свойства и характер поведения решений линейных систем с постоянными коэффициентами. Они охватывают варианты устойчивых и неустойчивых лино]...ьх систем, жестких и нежестких систем, характеристические уравнения которых содержат действительны.», комплексные или кратнь» корни.
Результаты исследований сравнивается с результатами, полученными одной : - наиболее распространенных програми sriFP. Проведен анализ результатов, подтверждающий ripe, «уцества ДРМ как в обеспечении требуемой точности, так и времени решения.
Для примера, в таблице 1 приведены сравнительные результаты решений вариантов теста, заимствованного из. работы Араиунян 0."., Залеткин С.Ф., Захаров А.Ю.. Калиткин Н.Н. О тестировании программ решения обыкновенных дифференциальных уравнений: Препр./Ин-г прихл. матем. им. (.Келдыша'АН СССР/, 1983.-Г 139. -?0 с. Шесть вариантов теста охватывают различные типы осциллирующих задач: плох' обусловленную, мягкую, оызтроосциллирую-щую, жесткие и жестко-осциллирующую. . В. таблице приняты позначення: Е - максимально допустимая погрешность решения; КШ - обіцее количество тагов, интегрирования; КГ! - максимальная
погрекность приближенного решения на интервале интегрирования;
Таблица 1
КШП - количество иагов, о погрепностью, лровыгащей допуотииую; Т - время решения варианта задачи в секундах.
B 1.^.3 рассмотрены методы решения линейных систем с постоянными коэффициентами произвольного к -го порядка к-1
(38)
г(Ю
i-1 k_i
о начальными условиями /^ (0) = Y^*", . (3=0..k-1),
где Yelr, BeR , коэффициент A j - матрицы размерности NxN.
Эти методы описываются рекуррентными соотношениями
F' .Г V (D*.1+J5,Ak-1-mG0.i>'Fd.i=D3.i+apWd+1.l' 1^-^1>
4Р) - ^-ік - '$
ю=1
U-1
Р(р)
О, (1=1..р)
(39)
0р= Е - е ajj hhk_j. Gm>i - уКБ^^-Гс^д). (m=1..k-1)
"o.l
0.
u+1 n p
r(p>
где Si (3=1.-p) определяются по формуле из (34),
На каждом этапе рекуррентного процесса необходимо решать систему линейных алгебраических уравнений с общей для всего процесса матрицей CL размерности NxN.
Разработан алгоритм определения локальной погрешности методов и их. реализации с автоматическим изменением шага интегрирования : порядка метода. Составлен программный модуль EMLSK для экспериментальных исследований, проведенных на трех тестовых задачах, построенных автором и охватывающих различные типы поведения решения. Сравнительный анализ результатов исследований подтвердил преимущества ДРИ над другими современными численными иетодаьм.
Глава 2 "Одношаговыэ методы дробно-рационального вида",
поовящена методике построения и обоснования дробно-рациональных
одношаговых численных методов на базе теории ДРП из первой гла
вы. _ /
В 2.1 рассмотрены принципы построения дробно-рациональных олношагозых методов с использованием линейных явных методов Ру-нге-Кутты. Обоснованы сложности построения методов высших
С целью уменьшения вычислительных затрат в 2.2 разработана методика построения линейных явных однотаговых численных петодов решения автономных систем дифференциальных уравнений первого порядка с использованием матрицы Якоби J правой части системы. Приведены примеры таких методов до четвертого порядка точности. Построение по данной методике методов более высоких порядков также является сложным и громоздким.
В 2.3 изложены новые общие принципы построения вложенных явных линейных одношаговых численных методов высоких порядков, которые базируются на следующем утверждении.
Утверждение 1. Если функция Г(хЛ) системы (1) непрерывно дифференциируема не менее p+t раз и У (х^ + ал) -тейлоровское приближение р-го порядка решения системы (1) относительно сеточного узла х^, то Під + ah, Ур(хп + ah)) определяет разложение в ряд Тейлора р-го порядка производной решения относительно узла Хд, т.е.
Цх^ ah. Yp^+ah)) = Y^ + олУ +...+-^1^1^ 0{др+1).
Параграф состоит из двух подпунктов.
В 2.3.1 приведена методика построения вышеуказанных методов для неавтономных систем. Построение методов р-го порядка состоит в определении приближений решения в последовательности р+1 произвольных точек Хд + ctjh, (1=) ,р*1). В*едя обозначения
У*п +ah> = Yn?d ФІ = Тр.п эти приближения определяются соотношениями
Формулы (40).(4-1) определяют рекуррентные соотношения о набором свободных параметров a_i(p=1.2,... ,1=1, .р+1), а01=0. Для их реализации найдены формулы определения значении . коэффициентов Р у, которые зависят от параметров a^: ppl;j « Ppij^j)-
Если при Opi=1 ввести обозначение Ppi;jd )=Tpj. то соотношение (40) определяет приближение р -го порядка точности в сеточном узле xn+1=in+h,T.e.
Тем самим, соотношения (40-42) описывают рекуррентный процесс по следовательного построения многопараметрического семейства методов от пергэго до требуемого порядка.
Приведены примеры катодов построеных по указанной кот..;.::->.
В 2.3.2 использована методика из предыдущего подпункта для построения последовательности вложеных явных линейных одношаго-вых методов для <шгономных систем. Рассмотрены их особенности, которые позволяют уменьшить общоо количество вычислений правой ччсти системы дифференциальных уравнении.
С целью укэн мнения среднего количества вычислений правой части системы на саг интегрирования, приведенная выше методика применена в 2.4 к построению явных блочных однобитовых методо .. При общих количествах вычислений.правой частя, предполагаемых ранее на один шаг интегрирования, методом р-го порядка в блочных методах рекуррентный процессом определяются приближе. :я решения . непосредственно в р сеточных узлах . Последовательны» .этапы процесса увеличивают глубину блока сеточньк узлов и повышают порядок, приближений решений. Однако, в случае сильно жестких систем приближения решения лилейными блочнььи методами, за счет их неустойчивости, при расширении глубины блока, могут быть сильно искажены, что понижает" эффективность их дробно-рационального усреднения. Поэтому, блочные методы является наиболее зффектизньйя при их непосредственном использовании для исследования нежестких задач.
В 2,5 на осново ранее описанной методики исследуются вопросы построения двухиаговых явных числимых методов.
Пусть" Уда'--твЯя>роа>коо приближения решения р-го порядка в некоторой промежуточной точке Зд+оЖ, где J+1. Тогда справедливо следующее
Утверждение 2.Приближение *п+« согласуется с приближением ^g- с точностью Otti1»1), т.е. YJ&J* = *W' + 0(hp+1).
Указанное утверждение позволяет . при построении линейных ад' ленных методов (42) частично использовать коэффициенты K_j найденные ранее в пр.дыдущем узле, что существенно сокращает количество вычислений правой части, системы. .,
Изложена методика построения таких методов произвольного р .
-го порядка. Количество вычислений правой части на каздом иаго интегрирования равно порядку метода.
В заклкмешш ко второй главе приведен анализ эффективности использования дробно-рациональных одношаговых методов. Значительные вычислительные затраты при реализациии методов высоких порядков ставит под сомнение целесообразность их применения. Поэтому, их рекомендуется использовать для начального процесса реализации многошаговых методов.
Глава 3 "Многошаговые котоды др^оно-рационального вида" заключает в себе изложение основ теории дробно-раииональных многошаговых численных штодов. Она содержит важные прикладные аспекты теории ДРП.
В 3.1 приведена методика построения линейных иногопаговых численных кетодов с переменньа сагой интегрирования. Она базируется на принципах построения одношаговых методов, рассмотрены* во второй главе, о пыбором опорными раное определенных сеточных узлов резания. Пра этом учитшается порядок найденных ранее приближений ревония в предыдущих узлах шаговой
СОТКИ Зд. Зп_1. Хд.2
Используя обозначения
* - *n+rV hJ - Vi-rVr sl - k**y (Ы--P-1 > (43) линейньй многопаговьЯ ?.йтод р-го порядка представим в виде
Определены обцив раоч&тнкэ формулы коэффициентов метода произвольного р-го порядка из скстеиы согласованности кегода
SV1- ^3^-^. (M..p-t). (45)
реиенио которой описывается вмра*1»!ЖИ
Рй>я (-1)р^ l3g4 : . -Ot.-p-f). «б)
(k)' v*
где &A определяются рвкуррентнши соотношениями
-26-(k)_ (k-i)+ .(k-1 )._
f Sk_,, J >. k-1 *3 I Sk , J * k-1
Sot)= '» 0)= ' 1=1--P"2- 3=1.-P-1. k=2..p-l.
Приведены примеры таких методов до четвертого п^мдла. Следует заметить, что, если интегрирование проводить с постоянным шагом, т.е. 3^=1» то метод (44) совпадает с методом Адамоа-Банфорта того же порядка.
Указаны подходы к построению на базе методов (44) многоша-. овых дробно-рациональных методов.
В 3.2 излагается методика построения линейных явных многошаговых методов с переменным шагом интегрирования для автономных "четеи с использованием значений второй производной реше-
ГДЄ KCp-1] _ гіуІР-ІІ). КІР-23 - J(YfP-?]i-K[p_2)
V1 + V2 = p-2 .. ' . (49)
Параграф состоит из трех подпунктов.
В 3.2.1 для выделения из описанного выше множества многошаговых методов, путем введения дополнительного условия viat»2t приведена ..етодика построения методов четных порядков.
В 3.2.2, при общем условии (49), рассмотрено два варианта дополнительных условий: a) vz - V1-1 и б) v2= vi+1 при построении методов, нечетных порядков. Приведена расчетные формулы КО г офЗ'.ііу'.ентов мє" ">дов различных порядков.
В 3.2.3 описана методика построения отдельного класса многошаговых методов вида (48), использующих в ранее найденных сеточных узлач только вторую производную решения (г>1=0). Получены обцие формулы определения коэффициентов для . методов произвольного порядка.
В 3.3 определены оценки локальной погрешности линейных многошаговых методов общего вида (48).
Доказана следующая Теорема 8 Локальн'т погрешность є многоша. ового метода р-го порядка типа (48) определяется соотношением
V у!р+1} (a)hP+1 (50)
-27-rnf- r-v - Swh < г < Хд + h . v =max(vl,v2).
p^p^-S^P-1)]. (51)
Параграф состоит из двух подпунктов.
В 3.3.1 исследованы оценки локальной погрешности методов вида (44) для неавтономных систем, а в 3.3.2 оценки для методов, использующих только вторую производную (v1=0).
В 53.4 исследуются вопросы построения многошаговых дробно-рациональных методов (МДМ) о переменным шагом интегрирования. Он состоит из двух подпунктов.
В 3.4.1 рассмотрены МДМ вида
ЇЙ - -И; : : . (52)
Z (-1)\іг\Г*
где *n+i " приближение решения в сеточном узле zJ]+1 линейным многощаговым кетодом k-го порядка, а также их частные случаи на базе структуры (19). Исследуется влияние структуры линейных иногопаговых методов на устойчивость ИДМ. Предлагается использовать в МДМ для каждого метода отдельного р-го порядка (р=1,2,...) свою последовательность лилейных иногокаговых методов вида
J=0
(53)
Y&W* =5oaS-A-i **%^№№-з » (w..p-п
п+1 п. j_r> іР л~3
где a|y (1=1..3) - свободные коэффициенты. .
Т&<У2Эмя_9. Если коэффициенты dj?' методов (53)
определяется соотношения?*!
_(р) ш lc-1P 1р-К>1 „ „ . P;kn(p) /,,4 «_k,
(54)
aip-k ^ ' %-k= 1 "jE aiP-k (1-1.-Р-к)
где с^. параметры метода (52), то характеристическое уравнение МДМ порядка р имеет только один отличный от нуля корень, т.е. не содержит "паразитических" корней.
Исследованы условия устойчивости .МДМ и доказана
Теог^ма 10. Операторная функция D (а) МДМ (52) порядка р при условии (54) не зависит от значений коэффициентов используемых линейных шогошаговых методов (53) .
Из теоремы следует, что область устойчивости МДМ совпадает ' с областью устойчтосг базовых ДРП (4).
Установлены общие закономерности определения коэффициентов МДМ произвольных порядков.
В 3.4.2 определены оценки локальной шгрешости R_=0e В МДМ .
W^li^i-i^^Vi^ v
0).=
1 р
^п .h W
M>vi
їм *«
Е (-l)Ja,h^
j«o
Значення С' определяет
Тес та II.Локальная погрешность ладейных многооаговых м тодов вида (53) определяйся соотношениями
4-І - сДї(Р-Ь+1 Чг-)^1, г « (vVk11' Vh)-
р(Р> -. 1, ...Ги, с 1Р-к+1_Ргк"1(Р) (с ч \Р-*+1_
-«р-^^Й-к^р-^^]. <56)
Лооле преобразований, коэффициенты С* представиш в видя
(57) Рассмотрены вопросы практической реализации вычислений локальной погрешности.
В 3.5 описан алгоритм реализация МДМ соотнесениями
і; описі..ш экспериментальные иослпедования эффективности мето-
дов. Он состоит из трех подпунктов.
В 3.5.1 разработан алгоритм определения коэффициентов <х|' и р^' в методах (58).
В 3.5.2 построен алгоритм изменения шага интегрирования и порядка метода. При выполнении условия IL < &гххх, где єпах-ма-ксималым допустимая погрешность, прогнозируемая величина после дупгего тага п. определяется по формуле
^'^ Ч-тр . .
Алгоритм изменения порядка базируется на определении значения р, при которой .локальная погрешность минимальна.
В 3.5.3 описаны экспериментальный исследования эффективности МДМ. Для этой цели разработан програютьй модуль ВЇЇН1 и подобрано 25 тестовых задач, которые вклхмаот линейные и нелинейные системы ОДУ с различными свойствами: жестки», нежесткие,устойчивые и неустойчивые. Для исследования качества решений внутри интервала интегрирования построены 4 нелинейные тестовые задачи с известный! аналитическими точнюи решениями. Приведены результаты экспериментальных исследований с точностям є =10 , 10 ,1 (Г , которые систематизированы в отдельных таблицах. Таблицы содержат основные, наиболее ваянда, количественные и качественнее характеристики: общее чиоло иагов на интервал интегрирования, количество вычислений правой части системы, количество вычисления матрицы Якоби системы, максимальная глобальная погрешность на интервале интегрирования, количество патов о погрешностью, превышающей допустимую, погрешность в конце интервала интегрирования, время реиения варианта задачи в сек. Все показатели в таблицах приведены в сравнении с результатами решений задач известной программой STIFF.
Анализрезультатов тестирования показал, что ШМ1 по всем приведенным характеристикам существенно эффективнее программы STIFF. В отличие от программа STIFF, в тестовых задачах , для которых известны точные решения, DKM1 гарантирует на всем интервале интегрирования требуемую точноогь при меньшем количестве шагов, вычислений правой части и матрицы Якоби. Во всех тестовых задачах наблвдается значительный выигрыш по времени решения. Для отдельных задач в десятки и сотни раз.
В заключении обобщены основные результаты работы.
В приложении приведены программные модули: построение области устойчивости ДРП (ОСОП), реализации дробно-рациональных
методов гашения линейных систем первого порядка (DML51), решения линейных систем k-го порядка (DJQSK), реализации дробно-рациональных многошаговых методов (ШМ1). Также s приложении описаны тестовые задачи, заимствованные из других работ.
