Содержание к диссертации
Введение
1 Методы второго порядка точности 12
1.1 Основные понятия и определения 12
1.2 Численная схема второго порядка для автономных систем 15
1.2.1 Условия второго порядка точности 15
1.2.2 Исследование устойчивости 16
1.2.3 Связь схемы второго порядка с явным методом Рунге - Кутта и
(га, к) - методом 18
1.2.4 Контроль точности вычислений 19
1.2.5 Контроль устойчивости 20
1.3 О применимости к неавтономным задачам схем. построенных для автономных задач 22
1.4 Численная схема второго порядка для неавтономных систем 25
1.5 Анализ результатов расчетов 27
2 Методы третьего порядка точности для автономных систем 28
2.1 Численная схема 1 29
2.1.1 Условия третьего порядка точности 29
2.1.2 Исследование устойчивости 32
2.1.3 Контроль точности и устойчивости 38
2.1.4 Связь схемы 1 с явным методом Рунге - Кутта и (т, к) - методом . 41
2.2 Численная схема 2 42
2.2.1 Условия третьего порядка точности 42
2.2.2 Исследование устойчивости 45
2.2.3 Контроль точности и устойчивости 54
2.3 Численная схема 3 57
2.3.1 Исследование устойчивости 58
2.3.2 Контроль точности и устойчивости 60
2.4 Анализ результатов расчетов 63
3 Методы третьего порядка точности для неавтономных систем 64
3.1 Численная схема 1 G5
3.1.1 Условия третьего порядка точности 65
3.1.2 Исследование устойчивости 70
3.1.3 Контроль точности и устойчивости 75
3.2 Численная схема 2 78
3.2.1 Условия третьего порядка точности 78
3.2.2 Исследование устойчивости 83
3.2.3 Контроль точности и устойчивости 87
3.3 Анализ результатов расчетов 89
Заключение 90
Список использованных источников
- Численная схема второго порядка для автономных систем
- О применимости к неавтономным задачам схем. построенных для автономных задач
- Контроль точности и устойчивости
- Условия третьего порядка точности
Введение к работе
При решении ряда задач, таких как проектирование радиоэлектронных схем, моделирование кинетики химических реакций, расчет динамики механических систем и других возникает необходимость численного решения задачи Коши для жестких систем обыкновенных дифференциальных уравнений j/ = f(t,y), y(to) = Vo, t0 t tk, (0.0.1) где у и / - вещественные N - мерные вектор - функции, уо - начальное условие, t -независимая переменная, которая меняется на заданном конечном интервале. Учет все большего числа факторов при построении математических моделей физических процессов, имеющих компоненты с сильно различающимися временными константами, приводит к жестким системам обыкновенных дифференциальных уравнений все большей размерности. Кроме того, для некоторых задач правая часть может быть разрывной.
Несмотря на рост быстродействия ЭВМ, сложность задач, возникающих на практике, опережает развитие вычислительной техники, что повышает требования к вычислительным алгоритмам. Поэтому проблема создания эффективных численных методов решения задачи Коши для жестких систем большой размерности является актуальной.
Для построения эффективных алгоритмов интегрирования необходимо 1) выбрать методы, соответствующие классу решаемых задач, 2) решить ряд вопросов, связанных со способом управления величиной шага и оценкой точности получаемых результатов.
Каждый из предложенных в данной работе методов основывается на двух численных формулах - явной и неявной. Из них строится метод, называемый аддитивным [1-4,50-55]. Для применения аддитивного метода правая часть / исходной дифференциальной задачи (0.0.1) разбивается на две части f{t,y) = p{t,ij) + g(t,y), которые будем называть "нежесткой" и "жесткой", соответственно. Аддитивный метод строится так, что явный метод используется для решения задачи с "нежесткой" частью p(t,y), а неявный - с "жесткой" g(t,y). Поэтому, говоря о методе, "жесткую" часть g{t,y) будем называть неявный частью, а "нежесткую" часть tp(t, у) - явной. В качестве неявной составляющей аддитивного метода используются (т, к) - методы [5-20], а в качестве явной - явные методы типа Рунге - Кутта. Таким образом, в основу алгоритмов интегрирования положены одноша-говые безытерационные методы вида Уп+і = Уп + h pf(tn, уп, h), п = 0,1,2,..., где начальное условие у0 задано, п - текущая точка интегрирования, h - шаг интегрирования, (р/ - заданная вектор - функция, зависящая от правой части исходной задачи. В данной форме можно представить не только безытерационные методы, но и неявные итерационные методы типа Рунге - Кутта с фиксированным итерационным процессом, в котором число итераций не зависит от номера шага интегрирования. Одношаговые безытерационные методы имеют следующие преимущества перед многошаговыми методами:
• многошаговые методы усредняют решение ("срезают экстремумы"), что для некоторых задач неприемлемо,
• многошаговые методы малоэффективны для задач с разрывной правой частью,
• безытерационность метода позволяет оценить вычислительные затраты на шаг интегрирования до проведения расчетов,
• отсутствие итерационного процесса значительно упрощает программную реализацию алгоритма интегрирования.
Достаточно полный обзор работ по решению (0.0.1) многошаговыми методами содержится в [25,26, 29, 30,148-162]. Поэтому далее на многошаговых методах останавливаться не будем.
Многие современные способы управления величиной шага основаны на контроле точности численной схемы. Обзор способов оценки ошибки приведен в [21,25-27,32-48]. Здесь для оценки глобальной ошибки метода р - го порядка точности используется вложенный метод (р — 1) - го порядка, а неравенство для контроля точности имеет вид \Уп,р Уп,р—1\ .. max —г-—І е. i i N \у ПіР\ + г где уп р и уп,р \ - соответственно, приближения к решению в точке tn методами р - го и (р — 1) - го порядков, г - положительный параметр. Если \у1п\ г, то контролируется относительная ошибка є, в противном случае контролируется абсолютная ошибка є-г. Указанный способ хорошо зарекомендовал себя в [29,30,32-37,44-48] и будет использоваться ниже.
Поскольку в аддитивные методы входят явные численные формулы, возникают определенные проблемы с устойчивостью численной схемы. При решении жестких систем явными методами на участках установления, которые, как правило, составляют большую часть интервала интегрирования, шаг ограничен не требованием точности, а требованием устойчивости численной схемы. При выборе величины шага исходя только из требования точности на этих участках возникает неустойчивость, приводящая к раскачиванию шага, что в лучшем случае снижает эффективность алгоритма интегрирования. Этот недостаток можно устранить дополнительным контролем устойчивости численной формулы [21]. Поскольку неравенство для контроля устойчивости грубое, то оно используется как некоторый ограничитель на рост шага. В результате прогнозируемый шаг hn+i по точности с ограничением но устойчивости вычисляется по формуле hn+i = max{hn, min{/iacc,/ist}}, где hacc - шаг, выбранный исходя из требования точности, hst - шаг, выбранный исходя из требования устойчивости, hn - предыдущий шаг.
Обычно для решения умеренно жестких задач применяются явные методы с расширенными областями устойчивости [81-135]. Однако явные методы имеют ограниченную область устойчивости, следовательно, при решении задач высокой степени жесткости требование устойчивости накладывает слишком обременительное ограничение на размер шага интегрирования [25]. Поэтому, явные методы применяются для решения задач средней жесткости [21].
Для решения жестких задач обычно используют L - устойчивые методы [5-20,56-76]. В случае большой размерности системы дифференциальных уравнений в L - устойчивых методах общие вычислительные затраты фактически полностью определяются временем вычисления и обращения матрицы Якоби. Затраты можно значительно уменьшить, если использовать одну и ту же матрицу Якоби на нескольких шагах интегрирования, то есть применять алгоритмы с замораживанием матрицы Якоби. В итерационных методах замораживание матрицы не влияет на порядок точности численной схемы, а определяет скорость сходимости итерационного процесса. Поэтому данный подход широко используется при реализации методов такого типа [26]. В безытерационных методах, к которым относятся методы типа Розенброка [77] и их различные модификации (например, [9], [32]), матрица Якоби включена непосредственно в численную формулу. Поэтому ее аппроксимация может приводить к потере точности численной схемы, что приводит к определенным проблемам с ее замораживанием. Методы типа Розенброка для автономной системы У = /Ы, У (to) = 2/о, k t tk, имеют вид т i-l Уп+i =Уп + Y Piku Dni = E- aihf (yn + Y lijkj), i=l j=\ i-l Dnik = hf(yn + pijkj) i=\ где E - единичная матрица, / = df/ду - матрица Якоби правой части, Pi,ai)Aj 7ij параметры. Наиболее эффективные реализации методов такого типа возникают тогда, когда все щ равны между собой и 7ij = 0. В этом случае на каждом шаге требуется вычисление и обращение только одной матрицы Dn = Е — ahf (yn).
Проблема замораживания матрицы Якоби значительно проще решается в рамках (т, к) - методов [5-20]. Заметим, что (гп, к) - методы так же просты в реализации, как и методы типа Розенброка, а свойства точности и устойчивости у них лучше [16]. Класс (т, к) - методов вводится следующим образом. Зададимся целыми числами т и к, к т, и рассмотрим следующие множества A/m = {l,2,...,m}, Мк — {ті е Л/т1 = Ш\ ГТІ2 • гпк т}, Ji = {rrij — 1 Є Mm\j 1, rrij Є Л/jt,raj г} 2 i m. Тогда (m, A;) - схемы имеют вид Уп+і - Vn + Piki + р2к2 + - + Pmkm, Dn = E- ahf n, t-i Dnki = hf{yn + /%) + ay ,-, і Є Л/ , (0.0.2) Dnki = ki-i + J afjfcj, і Є Л/т \ АД.
Здесь Е1 - единичная матрица, f n = д/(уп)/ду - матрица Якоби функции f,h- шаг интегрирования, &,-, 1 і m, - стадии метода, а,рі,(3ц,ац вещественные константы, определяющие свойства точности и устойчивости (0.0.2). Для описания традиционных одношаговых методов достаточно одной константы m - числа стадий. В данных схемах для описания затрат на шаг необходимо введение двух чисел m и к. Вычислительные затраты на шаг интегрирования в методах (0.0.2) следующие - один раз вычисляется матрица Якоби и осуществляется декомпозиция Dn, к раз вычисляется функция /, тп раз осуществляется обратный ход метода Гаусса.
Отметим, что при к = m и aij = 0 данные методы совпадают со схемами типа Розен-брока. В этом случае (к, к) - схему (0.0.2) можно поставить в соответствие к - стадийной полуявной численной формуле типа Рунге - Кутта, при реализации которой на каждом шаге используется одна матрица размерности N. Относительно таких формул известно, что нельзя построить к - стадийную схему выше (к + 1) - го порядка точности. Отметим, что схема максимального порядка может быть только Л - устойчивой, что недостаточно для построения эффективных алгоритмов интегрирования. Поэтому в настоящее время применяются методы типа Розенброка, в которых порядок точности совпадает с числом вычислений правой части дифференциальной задачи. Однако если рассматривать схемы (0.0.2) при m к, то при к равном 1 можно построить L - устойчивый (2, 1) - метод второго порядка, а при к равном 2 и 3 можно построить L - устойчивые методы (к + 2) - го порядка точности, соответственно.
Независимость порядка точности (т, к) - метода от числа шагов с замороженной матрицей достигается за счет того, что при получении условий порядка предполагается, что при реализации метода вместо матрицы Якоби f n правой части системы дифференциальных уравнений применяется матрица Л„, представимая в виде Ап = f n + hlin, где Rn -некоторая матрица, не зависящая от шага интегрирования. Это естественное требование в случае применения одной и той же матрицы Якоби на нескольких шагах интегрирования. Отметим, что при замораживании матрицы Якоби максимальный порядок точности L -устойчивой численной схемы с к вычислениями правой части дифференциальной задачи равен (к + 1).
Попытка повысить эффективность алгоритма интегрирования отделением жесткой и нежесткой компонент решения предпринята в проекционных методах. Идея проекционных методов состоит в разделении жесткой системы у — f(y) на жесткую систему небольшой размерности у а = /а(Уа,Уь) и нежесткую систему большой размерности у ь = /ь(Уа,Уь) с целью применения к жесткой системе неявного метода, а к нежесткой - явного. Трудность этого подхода заключается в проблеме разделения на жесткую и нежесткую системы.
Построенные в данной работе аддитивные методы предусматривают разбиение правой части f(t, у) на жесткую g(t,ij) и нежесткую p(t,y) части t/ = 4 (t,y) + 9(t,y), у{к) = Уо, k t tk. (0.0.3)
Предполагается, что большая часть жесткости сосредоточена в функции g{t,y). Методы конструируются так, чтобы неявный метод использовался для решения задачи с жесткой частью, а явный - с нежесткой. В качестве неявной составляющей аддитивного метода используется (т, к) - методы [5,13-16,18], а в качестве явной - явные методы типа Рунге -Кутта. При разбиении правой части задачи (0.0.1) на жесткую и нежесткую части можно выделить два случая
1) Естественное разбиение. При численном решении ряда задач, описываемых диф ференциальными уравнениями в частных производных, после дискретизации по пространственным переменным (методом конечных разностей или конечных элементов) возникает необходимость решения задачи Коши для системы обыкновенных дифференциальных уравнений вида (0.0.3), в которой правая часть естественным образом разбивается на жесткую и нежесткую составляющие. При этом разбиении выполняется предположение о том, что большая часть жесткости сосредоточена в функции g(t, у). Например, в некоторых задачах механики сплошной среды жесткая часть появляется после дискретизации дифференциального оператора второго порядка, а нежесткая - после дискретизации дифференциального оператора первого порядка. Кроме того, в некоторых задачах вида (0.0.3) матрица Якоби функции g(t, у) имеет специальный вид (например, в задачах механики сплошной среды она симметрична).
2) Искусственное разбиение. Полагая g(t,у) = By, p(t,y) = f(t,y) — By, где В — df(t,y)/dy, можно считать, что вся жесткость сосредоточена в функции g(t, у) [1-4, 50-53]. Тогда система (0.0.1) примет вид y = [f{t,y)-By} + By, (0.0.4)
В данной работе и в [4,50-53] в качестве В может быть использована произвольная аппроксимация матрицы Якоби df(t,y)/dy. При "разумной" аппроксимации этой матрицы можно считать, что основная жесткость задачи сосредоточена в функции g{t,y). Однако считать, что вся жесткость сосредоточена в функции g(t,y), а ip(t,y) есть нежесткая часть, вообще говоря, нельзя. Поэтому все построенные в данной работе методы оснащены неравенством для контроля устойчивости явной части численной схемы, которое в ряде случаев повышает эффективность алгоритма интегрирования. Контроль устойчивости явной части может оказаться полезным и для случая естественного разбиения правой части. Если есть уверенность, что вся жесткость сосредоточена в функции g(t, у), то контроль устойчивости следует отключить, потому что в некоторых случаях он может приводить к понижению эффективности алгоритма интегрирования.
Условия порядка для всех предложенных здесь схем выводились без каких-либо предположений относительно вида функций p(t, у) и д(t, у). В результате, для построенных аддитивных методов несущественно, естественным образом или искусственным была разбита правая часть, а имеет значение лишь то, как распределена жесткость задачи относительно функций p(t,y) и g(t,y). В идеале вся жесткость должна быть сосредоточена в функции g(t,y). Возможность произвольной аппроксимации матрицы Якоби исходной задачи, без снижения порядка точности, в аддитивных методах достигается за счет произвольности представления задачи (0.0.1) в виде (0.0.4). В частности, метод допускает замораживание матрицы Якоби.
Хотя совершенно произвольная аппроксимация матрицы Якоби не снижает порядка точности метода, она может понизить его эффективность.
1. Например, для задачи (30) из [26, с. 167] матрица В — diag(df(y)dtj) слишком грубо показывает направления изменения решения из-за особенности этой задачи - когда компонента решения у2 случайно становится отрицательной, она стремится к —со. В результате, для данной задачи алгоритм интегрирования с диагональной аппроксимацией матрицы Якоби менее эффективен, чем с полной матрицей Якоби, несмотря на экономию вычислительных затрат при ее обращении.
2. Несмотря на контроль устойчивости явной части численной схемы, алгоритм интегрирования будет достаточно эффективен, когда основная часть жесткости сосредоточена в функции g(t,y), соответствующей неявной части схемы. Например, теоретически возможным, но практически бесполезным, было бы положить 5 = 0. Тогда g(t,y) = 0 и вся нагрузка легла бы на явный метод Рунге - Кутта, который, по указанным выше причинам, малоэффективен для решения жестких задач.
3. Диагональную аппроксимацию целесообразно использовать для матрицы Якоби с диагональным преобладанием.
Поэтому, говоря о возможности произвольной аппроксимации, предполагается ее "разумность", то есть соответствие классу решаемых задач.
Цель настоящей работы состоит в разработке эффективных алгоритмов численного решения задачи Коши для жестких систем обыкновенных дифференциальных уравнений большой размерности. Построенные алгоритмы призваны расширить область применения, занимаемую явными методам (возможно, обладающими расширенной областью устойчивости и оснащенных неравенством для контроля устойчивости [21]), в сторону задач большей степени жесткости, а не конкурировать с методами, использующими полную аппроксимацию матрицы Якоби. При этом их вычислительные затраты на шаг интегрирования должны быть сопоставимыми с явными методами.
В основе построенных алгоритмов интегрирования переменного шага лежат аддитивные схемы, которые, в свою очередь, конструируются из двух схем - явной и неявной. В данной работе в качестве явной схемы используется явный метод тина Рунге - Кутта, а в качестве неявной - (т, к) - метод. Поскольку оба метода являются одношаговыми безытерационными методами, то и аддитивные методы можно отнести к одношаговым безытерационным.
В исходной задаче Коши правая часть может быть представлена как в виде (0.0.1), так и (0.0.3). В первом случае задача приводится к виду (0.0.4), то есть производится искусственное разбиение правой части. Во втором случае правая часть уже разбита исходя из некоторых соображений предметной области, согласно которым большая часть жесткости сосредоточена в функции g(t,y), то есть разбиение было произведено естественным образом.
В случае естественного разбиения правой части эффективность алгоритма интегрирования достигается за счет того, что отделение более жесткой части g(t,y) от менее жесткой p{t,y) позволяет использовать неявную схему, входящую в аддитивный метод, преимущественно для решения задачи у — g(t,y). Следовательно, неявная схема требует вычисление и обращение матрицы Якоби только функции g(t,y), а не всей правой части исходной задачи. Во многих прикладных задачах матрица dg(t, у)/ду, в отличие от d( f(t,y) + g(t,y))/dy, имеет специальный вид (например, в некоторых задачах механики сплошной среды она симметрична). В этом случае эффективность можно повысить за счет выбора специального метода решения линейных систем алгебраических уравнений, учитывающего информацию о виде матрицы Якоби (например, метода квадратного корня для симметричных матриц [31]), и, возможно, использующего параллельный алгоритм вычислений [136-147]).
В случае искусственного разбиения правой части затраты на обращение матрицы В можно значительно уменьшить за счет использования в качестве В не полной матрицы df(t,y)/dy, а некоторой ее аппроксимации специального вида и применения соответствующего этой аппроксимации метода решения линейных систем алгебраических уравнений. Применение аппроксимации матрицы Якоби специального вида не только уменьшит затраты на ее вычисление, но и позволит применить соответствующий метод решения линейных систем алгебраических уравнений, что может весьма значительно уменьшить затраты на ее обращение. В данной работе использовалась диагональная аппроксимация матрицы Якоби, поскольку именно в этом случае затраты на ее обращение наиболее низки и аддитивный метод по затратам на шаг интегрирования становится сопоставимым с явными методами, а по свойствам устойчивости, относительно неявной части, значительно превосходит их.
Возможно сочетание описанных способов, например, использование трехдиагональ-ной аппроксимации матрицы Якоби, применение метода прогонки для решения систем с этой матрицей и ее замораживание. При решении жестких задач на эффективность алгоритма интегрирования существенно влияют свойства устойчивости не только основной, но и промежуточных численных формул [21,26,77]. В данной работе этому также уделено немало внимания. Для многих задач эффективность алгоритма интегрирования можно повысить за счет контроля устойчивости явной части численной схемы.
Диссертация состоит из введения, трех глав, заключения и списка цитируемой литературы. Во введении обоснована актуальность темы, дан обзор основных работ и приведено краткое описание диссертации по разделам.
В первой главе на примере двух численных схем второго порядка точности демонстрируется общий подход к построению алгоритмов интегрирования переменного шага, в основе которых лежат аддитивные методы.
В первом параграфе вводятся основные определения, которые используются в дальнейшем. Второй параграф посвящен построению аддитивного метода второго порядка точности для автономных систем. Достаточно простой вид схемы позволяет пояснить некоторые особенности этих методов, характерные и для схем более высокого порядка. В пункте 1.2.2 вводятся некоторые определения, отражающие свойства только аддитивных методов. Пункт 1.2.3 раскрывает смысл понятия "аддитивный метод" и демонстрирует возможность "декомпозиции" аддитивной схемы на явную и неявную (с сохранением условий порядка точности). В пункте 1.2.4 предлагается способ получения неравенства для контроля точности, используемое для аддитивных методов второго порядка. В пункте 1.2.5 предлагается подход к построению неравенства для контроля устойчивости и стратегия управления величиной шага интегрирования, используемые для всех построенных в данной работе методов. В параграфе 1.3 обосновывается необходимость вывода численных схем отдельно для автономных и неавтономных систем обыкновенных дифференциальных уравнений. В параграфе 1.4 построен метод второго порядка для неавтономных систем. В параграфе 1.5 приведены результаты расчетов, подтверждающие эффективность построенных на основе аддитивных методов второго порядка алгоритмов.
Во второй главе исследуются три семейства численных методов решения задачи Копій для жестких аддитивных автономных систем. Построены алгоритмы интегрирования переменного шага с контролем точности вычислений. Приведены результаты расчетов, подтверждающие работоспособность и эффективность построенных алгоритмов.
Третья глава посвящена исследованию двух семейств численных методов решения задачи Коши для жестких аддитивных неавтономных систем. Построены алгоритмы интегрирования переменного шага с контролем точности вычислений. Приведены результаты численных экспериментов.
В заключении сформулированы основные результаты. В приложении приведены тестовые примеры.
Основные результаты диссертации опубликованы в работах [50-57].
В работе используется следующая система ссылок на формулы. Тройная нумерация означает номер главы, номер параграфа и номер формулы. Некоторые параграфы разбиты па пункты, нумерация формул внутри каждого параграфа сквозная.
Работа инициирована член-корреспондентом РАН Шайдуровым В.В., выполнена при финансовой поддержке грантов НШ-3428.2006.9 и РФФИ №05-01-00579, №06-08-920
Численная схема второго порядка для автономных систем
Рассмотрим задачу Коши для неавтономной системы обыкновенных дифференциальных уравнений y = f{t,y), y(t0) = y0, t0 t tk, (1.1.1) где у и / - гладкие вещественные N - мерные вектор-функции, t - независимая переменная.
Определение 1. Задача Коши (1.1.1) называется жесткой на интервале / С [to, ], осли для V t 6 I 1) Re(A2) 0, 1 і N, max Re(-Av) _ч Ki N mm Re(—A;) где Aj - собственные числа матрицы Якоби df(t, y(t))/dy. Далее в основном будем рассматривать жесткие задачи. Для численного решения (1.1.1) рассмотрим класс одношаговых безытерационных методов
Упн = Уп + h pf(tn, уп, /І), п = 0,1,2,..., (1.1.2) с заданным начальным условием /о, где h - шаг интегрирования, п - текущая точка интегрирования, pf - гладкая веіцественная N - мерная вектор - функция, зависящая от правой части / дифференциальной задачи. В форме (1.1.2) можно представить не только явные методы типа Рунге - Кутта, но и неявные и полуявные методы типа Рунге - Кутта с фиксированным итерационным процессом, в котором число итераций не зависит от номера шага интегрирования.
Обозначим через y(t) -точное решение (1.1.1), через уп - приближенное решение (1.1.1) в точке „, вычисленное но формуле (1.1.2). Определение 2. Величина еп = y(tn) - уп называется полной или глобальной погрешностью метода (1.1.2) в точке t = tn.
Определение 3. Метод вида (1.1.2) сходится, если для каждой задачи (1.1.1) имеет место тах0 „ л/ \\єп\\ - 0 при h - 0, и сходится с порядком р, если еп = 0(hp), где Mh = tk t0, - некоторая норма в R .
Определение 4. Погрешностью аппроксимации Sn+i схемы (1.1.2) в точке tn+\ Є [(о, h] называется величина, определяемая по формуле
Wi = y(tn+i) - y(Q - h pf(tn, y{tn), h).
Если yn = y{tn), то 6n+i — єп+\. Следовательно, погрешность аппроксимации есть ошибка, которая получается за один шаг интегрирования, поэтому во многих работах ее называют локальной ошибкой. Этот термин и будет использован далее.
Определение 5. Говорят, что метод (1.1.2) аппроксимирует (1.1.1), если при h -» 0 выполняется h l тах0 п м n -» 0 и имеет порядок аппроксимации р, если Sn = 0(hp+l). Перейдем теперь к определению устойчивости. Понятие устойчивости вводятся на линейном скалярном уравнении у = Ху, 1/(0) = 2/0, t 0, Re(A) 0. (1.1.3)
Будем предполагать, что если функция / линейна по у, то и (р/ линейна по у. Данное предположение не является обременительным, поскольку все практические методы обладают данным свойством. Применяя метод (1.1.2) к задаче (1.1.3), получим Vn+i = QWVn-Обозначим z = A/t. Функция Q(z) называется функцией устойчивости метода (1.1.2).
Определение 6. Метод (1.1.2) называется устойчивым для данного z Є С, если \Q{z)\ 1. Область R комплексной плоскости С называется областью устойчивости метода, если он устойчив при всех z Є R. Пересечение области устойчивости с вещественной осью называется интервалом устойчивости.
Определение 7. Одношаговый численный метод называется А - устойчивым, если его область устойчивости включает всю полуплоскость {z Є С : Re 2 0}.
Несмотря на свою простоту уравнение (1.1.3) успешно применяется для предсказания устойчивости практических численных методов решения нелинейных жестких систем вида (1.1.1), при этом параметр А трактуется как произвольное собственное число матрицы Якоби функции /. Кроме того, для уравнения (1.1.3) требование в Определении А - устойчивости является естественным. Действительно, при Re 2 0 для точного решения y(t) = eXty0 задачи (1.1.3) требование ехр(А) 1 выполняется всегда, то есть модуль точного решения является невозрастающей функцией, поэтому и от приближенного решения естественно потребовать выполнения того же свойства.
Свойства А - устойчивости оказывается недостаточно для эффективности алгоритма интегрирования, поскольку для многих А - устойчивых одношаговых методов \Q{z)\ — 1 при Re(2) - —со. В результате приближения к быстрозатухающим фундаментальным решениям могут затухать очень медленно. Определение 8. Одношаговый численный метод называется L - устойчивым, если он А - устойчив и \Q{z)\ - 0 при Re(z) -» —сю.
При решении задач с малым параметром эффективный метод интегрирования должен обладать не только L, но и D - устойчивостью. Определение D - устойчивости [28] вводится на некотором классе (обозначим его & ) линейных жестких задач с переменными коэффициентами вида y = A(t)y, у(0)=у0, 0 t T, (1.1.4) где вид матрицы A(t) определяется классом решаемых задач & и зависит от малого параметра е. Применяя метод (1.1.2) к (1.1.4), получим уп+} = M(tn,h)yn, где матрица M(tn,h) зависит от tn,h, A(t) и коэффициентов метода. Определение 9. Метод (1.1.2) называется D(y) - устойчивым, если (3 Л/)(ЗЛ)(У задачи вида (1.1.4) из класса &){V tn Є [0, T])(V/i Є (0, h])
О применимости к неавтономным задачам схем. построенных для автономных задач
В разделах 1.2 и 1.2.3 численная схема (1.2.4) строится для автономной задачи (1.2.1). Что делать, если решаемая задача имеет неавтономный вид y = p(t,y)+g(t,y), y(t0) = y0, t0 t tk? (1.3.1) Можно пойти двумя путями 1) построить численную схему для задачи (1.3.1), 2) привести задачу (1.3.1) к автономному виду и воспользоваться построенной ранее схемой (1.2.4).
Рассмотрим второй способ немного подробнее. Считая t зависимой переменной, добавим ос к вектору решений: у = Q, тогда к системе (1.3.1) следует добавить уравнение «- = 1. Положим тогда задачу (1.3.1) можно переписать в автономной форме У = Ф{У)+Ш, Ш = Уо, t0 t tk. (1.3.3)
Отметим, что численная схема, выведенная непосредственно для (1.3.1) будет существенно отличаться от схемы, полученной в результате применения (1.2.4) к (1.3.3). Отличие кроется в матрице Якоби, которая используется при построении всех рассматриваемых в данной работе методов. В первом случае матрица Якоби имеет вид п і? ,dg{tn,yn) ,дд(у) дд ип — an , а во втором ип — an
Следовательно, во втором случае вектор dg/dt будет участвовать в вычислениях, а в первом случае - нет.
Второй путь на первый взгляд кажется более привлекательным, поскольку выводить численные схемы для автономных систем легче, чем для неавтономных. Однако вхождение dg/dt в расчетные формулы приводит к ухудшению свойств D - устойчивости [28], что делает второй подход малопригодным для задач с малым параметром. Фактически все жесткие задачи являются задачами с малым параметром, при этом малый параметр играет роль параметра жесткости. В одних случаях он явно входит в правую часть, а в других - неявно [28]. В данной работе для решения неавтономных систем выбран первый подход.
Покажем на примере задачи из класса ,У (см. определение 10), что свойства D -устойчивости метода (1.2.4) второго порядка точности ухудшаются при выборе второго подхода к решению неавтономных систем. Так как ip никак не влияет на различие в схемах, полученных при первом и втором подходах (эти схемы существенно различаются только за счет dg/dt), то для простоты положим р = 0. Тогда, как отмечено в разделе 1.2.3, данная схема перейдет в (2, 1) - метод, а задача (1.3.1) перейдет в (1.3.7) . Переобозначив стадии hi, получим Уп+і = Уп+Рі кі+іЖ, Dnkx = hg(yn), (1.3.5) Dnk2 = kv Итак, пусть имеется схема (1.3.5) для автономной системы у = у{у), у(к)=Уо, h t tk. (1.3.6) и требуется численно найти решение неавтономной системы y = (j{t,y), y{to)=Vo, tQ t tk. (1.3.7) 1. Схема вида (1.3.5), но выведенная непосредственно для неавтономной системы (1.3.7), записывается так уп+1 = уп + Piki + р2к2, Dnki = hg{tn,yn), Dnk? = к\, где Dn = Е — ah[dg(yn)/dy]. Перепишем ее в более удобном для выкладок виде Уп+1 = Уп + Piki + р2к2, ki = hg(tn,yn) + ahd9 ,yn)ku (1.3.8) h = hg(tn, уп) + акЩ - (кх + к2). ду 2. Приведем задачу (1.3.7) к автономному виду У = 9(У), y(to) = Vo, t0 t tk. и применим к ней схему (1.3.5), получим Уп+i =Уп + Рікі+р2к2, h = hgijjn) + uh—{yn)ku ду к2 = hg(yn) + ah—(yn) Учитывая (1.3.2) и (1.3.4), перепишем полученную схему в старых переменных Уп+1 = Vn + V\h + lhk2, h = hg{tn,yn) +ah—{tn,yn)k1 Л-ah2 — {tn,yn), (1.3.9) Од 9g, Ot dy кг = hg{tn, yn) + ah—{tn, xjn) {hi + k2) + 2ah -r-(i„,yn). Исследуем свойства D(y) - устойчивости схем (1.3.8), (1.3.9) на модельной задаче из класса У , которая задается О 1 0 DW=I I - )=1-.(0 , Она имеет вид (1.1.4), где A[t) = є х I , Q j . 1) Для (1.3.8) матрица M(tn,h) (см. определение 9) имеет вид M{tn,h) = (є2 +{2а-рх - p2)h + а(а - pv)h2 \ (е + ah)2 a{tn){e{py+p2) + apih)h (є + ah)2
Учитывая, что а ф 0, легко видеть, что все элементы матрицы M(tn, h) ограничены при є - 0. Поэтому, хотя бы на данном семействе задач из классы И(У) оператор M(tn, h) равномерно ограничен по є и t„. 2) Для (1.3.9) имеем
Рассмотрим задачу Коши для неавтономной системы обыкновенных дифференциальных уравнений у = v(t,у) + g(t,y), у(к) = Уо, t0 t tk, (1.4.1) где у, р и д - гладкие вещественные N - мерные вектор - функции, t - независимая переменная. Будем полагать, что вся жесткость сосредоточена в функции g(t,y), a p(t,y) есть нежесткая часть.
При выводе формул потребуется разложение точного решения y{tn+\) в ряд Тейлора в точке tn h2 y{tn+i) = y{t-n) + h(y + gj + Y\(Pt + 9t + (Py P + (Py9 + 9y P + 9y9) + 0(1?), (1.4.2) где элементарные дифференциалы вычислены на точном решении y(tn). Для численного решения (1.4.1) рассмотрим четырехстадийную численную схему [4]
Контроль точности и устойчивости
На каждом шаге интегрирования вложенный метод требует дополнительно один обратный ход в методе Гаусса и не требует дополнительных вычислений правой части и матрицы Якоби. Поскольку на каждом шаге интегрирования в основном методе третьего порядка производится одна декомпозиция матрицы Якоби и четыре обратных хода в методе Гаусса, то для задач большой размерности увеличение вычислительных затрат, связанных с использованием вложенного метода, будет незначительным. Введем обозначение еп = гпах! 1 \угп — Уп,2І/(І2/пІ + г)- Тогда с использованием (2.1.17) точность вычислений метода (2.1.1) можно контролировать с помощью неравенства єп є, где г - положительный параметр. Если \угп\ г, то контролируется относительная ошибка є, в противном случае контролируется абсолютная ошибка Коррекцию асимптотического поведения (1.2.18) построенной оценки проводить не нужно, поскольку численная схема второго порядка точности (2.1.17) является L - устойчивой относительно неявной части.
Контроль устойчивости будем проводить точно так же, как описано в пункте 1.2.5, поскольку неравенство для контроля устойчивости довольно грубое и при вычислении прогнозируемого шага по формуле (1.2.20) используется только как ограничитель на рост величины шага интегрирования.
Ниже алгоритм интегрирования на основе схемы (2.1.1) с коэффициентами (2.1.15) будем называть ADD3 Требование, чтобы схема (2.1.19) имела третий порядок точности, приводит к четырем уравнениям, совпадающим со вторым, шестым, пятнадцатым и двадцатым уравнениями системы (2.1.3).
Подставляя () = 0 в (2.1.1) и учитывая, что в этом случае Dn — Е, получим б Схема (2.1.20) с точностью до переобозначений является явным трехстадийным методом Рунге - Кутта. Требование, чтобы схема (2.1.20) имела третий порядок точности, приводит к четырем уравнениям, совпадающим с первым, третьим, седьмым и десятым уравнениями системы (2.1.3).
Таким образом, чтобы схема (2.1.1) имела третий порядок точности, необходимо, чтобы образующих ее "подсхемы"(2.1.19) и (2.1.20) тоже имели третий порядок точности.
Попутно заметим, что схема (2.1.1) при tp = 0 является L - устойчивой, поскольку она А - устойчива (при р = 0 переходит в (т, к) - метод, который является А - устойчивым [81]) и удовлетворяет необходимому условию L - устойчивости относительно неявной части (1.2.9) в смысле определения 11, а следовательно и необходимому условию L - устойчивости при р = 0 в обычном смысле.
Для задач небольшой жесткости можно воспользоваться схемой, похожей на (2.1.1), но не содержащей Ді в четвертой стадии. Это слегка уменьшает вычислительные затраты на шаг интегрирования, но несколько ухудшает свойства устойчивости. Если потребовать, чтобы промежуточные численные формулы в четвертой стадии были L - устойчивыми, то естественно ожидать, что свойства устойчивости полученной схемы ухудшатся незначительно.
Для численного решения (1.2.1) рассмотрим шестистадийиую численную схему i=i где k{, 1 і 6, - стадии метода, а,р,-,сц,-,Ду,/?б./,7 числовые коэффициенты, определяющие свойства точности и устойчивости схемы (2.2.1). Запишем условия третьего порядка точности схемы (2.2.1). Для этого разложим стадии кі, 1 г 6, по степеням h до членов с h3 включительно и учтем (2.0.2), получим где элементарныИтак, на каждом шаге интегрирования вложенный метод требует дополнительно только один обратный ход в методе Гаусса и не требует дополнительных вычислений правой части и матрицы Якоби. Поскольку на каждом шаге интегрирования в основном методе третьего порядка производится одна декомпозиция матрицы Якоби (порядка 0(N3) арифметических операций) и четыре обратных хода в методе Гаусса (порядка 0(N2) арифметических операций), то для задач большой размерности увеличение вычислительных затрат, связанных с использованием вложенного метода, будет незначительным. Контроль точности (с использованием метода (2.3.10)), контроль устойчивости, а также выбор прогнозируемого шага производится аналогично пункту 2.1.3.е дифференциалы вычислены на приближенном решении Исследуются два семейства численных методов решения задачи Копій для жестких аддитивных неавтономных систем. При построении методов учтены условия согласованности. Построены алгоритмы интегрирования переменного шага с контролем точности вычислений. Приведены результаты расчетов, подтверждающие работоспособность и эффективность построенных алгоритмов.
Причины, по которым рассматриваемые в данной работе численные схемы выводятся отдельно для автономных и неавтономных задач, указаны в параграфе 1.3. Численная схема, выведенная непосредственно для (1.3.1) будет существенно отличаться от схемы, полученной в результате применения (1.2.4) к (1.3.3). Отличие кроется в матрице Якоби, которая используется при построении всех рассматриваемых в данной работе методов. В численной формуле, полученной для автономных систем, вектор dg/dt будет участвовать в вычислениях, а в разработанной непосредственно для неавтономных задач - нет. Очевидно, что строить методы для автономных задач существенно легче, чем для неавтономных. Однако вхождение dg/dt в расчетные формулы приводит к ухудшению свойств D - устойчивости [28], что делает данный подход малопригодным для задач с малым параметром. Достаточно много жестких задач являются задачами с малым параметром, при этом малый параметр играет роль параметра жесткости. В одних случаях он явно входит в правую часть, а в других - неявно [28].
Условия третьего порядка точности
При решении ряда задач, таких как проектирование радиоэлектронных схем, моделирование кинетики химических реакций, расчет динамики механических систем и других возникает необходимость численного решения задачи Коши для жестких систем обыкновенных дифференциальных уравнений j/ = f(t,y), y(to) = Vo, t0 t tk, (0.0.1) где у и / - вещественные N - мерные вектор - функции, уо - начальное условие, t -независимая переменная, которая меняется на заданном конечном интервале. Учет все большего числа факторов при построении математических моделей физических процессов, имеющих компоненты с сильно различающимися временными константами, приводит к жестким системам обыкновенных дифференциальных уравнений все большей размерности. Кроме того, для некоторых задач правая часть может быть разрывной.
Несмотря на рост быстродействия ЭВМ, сложность задач, возникающих на практике, опережает развитие вычислительной техники, что повышает требования к вычислительным алгоритмам. Поэтому проблема создания эффективных численных методов решения задачи Коши для жестких систем большой размерности является актуальной.
Для построения эффективных алгоритмов интегрирования необходимо 1) выбрать методы, соответствующие классу решаемых задач, 2) решить ряд вопросов, связанных со способом управления величиной шага и оценкой точности получаемых результатов.
Каждый из предложенных в данной работе методов основывается на двух численных формулах - явной и неявной. Из них строится метод, называемый аддитивным [1-4,50-55]. Для применения аддитивного метода правая часть / исходной дифференциальной задачи (0.0.1) разбивается на две части f{t,y) = p{t,ij) + g(t,y), которые будем называть "нежесткой" и "жесткой", соответственно. Аддитивный метод строится так, что явный метод используется для решения задачи с "нежесткой" частью p(t,y), а неявный - с "жесткой" g(t,y). Поэтому, говоря о методе, "жесткую" часть g{t,y) будем называть неявный частью, а "нежесткую" часть tp(t, у) - явной. В качестве неявной составляющей аддитивного метода используются (т, к) - методы [5-20], а в качестве явной - явные методы типа Рунге - Кутта. Таким образом, в основу алгоритмов интегрирования положены одноша-говые безытерационные методы вида Уп+і = Уп + h pf(tn, уп, h), п = 0,1,2,..., где начальное условие у0 задано, п - текущая точка интегрирования, h - шаг интегрирования, (р/ - заданная вектор - функция, зависящая от правой части исходной задачи. В данной форме можно представить не только безытерационные методы, но и неявные итерационные методы типа Рунге - Кутта с фиксированным итерационным процессом, в котором число итераций не зависит от номера шага интегрирования. Одношаговые безытерационные методы имеют следующие преимущества перед многошаговыми методами: многошаговые методы усредняют решение ("срезают экстремумы"), что для некоторых задач неприемлемо, многошаговые методы малоэффективны для задач с разрывной правой частью, безытерационность метода позволяет оценить вычислительные затраты на шаг интегрирования до проведения расчетов, отсутствие итерационного процесса значительно упрощает программную реализацию алгоритма интегрирования.
Достаточно полный обзор работ по решению (0.0.1) многошаговыми методами содержится в [25,26, 29, 30,148-162]. Поэтому далее на многошаговых методах останавливаться не будем.
Многие современные способы управления величиной шага основаны на контроле точности численной схемы. Обзор способов оценки ошибки приведен в [21,25-27,32-48]. Здесь для оценки глобальной ошибки метода р - го порядка точности используется вложенный метод (р — 1) - го порядка, а неравенство для контроля точности имеет вид где уп р и уп,р \ - соответственно, приближения к решению в точке tn методами р - го и (р — 1) - го порядков, г - положительный параметр. Если \у1п\ г, то контролируется относительная ошибка є, в противном случае контролируется абсолютная ошибка є-г. Указанный способ хорошо зарекомендовал себя в [29,30,32-37,44-48] и будет использоваться ниже.
Поскольку в аддитивные методы входят явные численные формулы, возникают определенные проблемы с устойчивостью численной схемы. При решении жестких систем явными методами на участках установления, которые, как правило, составляют большую часть интервала интегрирования, шаг ограничен не требованием точности, а требованием устойчивости численной схемы. При выборе величины шага исходя только из требования точности на этих участках возникает неустойчивость, приводящая к раскачиванию шага, что в лучшем случае снижает эффективность алгоритма интегрирования. Этот недостаток можно устранить дополнительным контролем устойчивости численной формулы [21]. Поскольку неравенство для контроля устойчивости грубое, то оно используется как некоторый ограничитель на рост шага. В результате прогнозируемый шаг hn+i по точности с ограничением но устойчивости вычисляется по формуле hn+i = max{hn, min{/iacc,/ist}}, где hacc - шаг, выбранный исходя из требования точности, hst - шаг, выбранный исходя из требования устойчивости, hn - предыдущий шаг.
Обычно для решения умеренно жестких задач применяются явные методы с расширенными областями устойчивости [81-135]. Однако явные методы имеют ограниченную область устойчивости, следовательно, при решении задач высокой степени жесткости требование устойчивости накладывает слишком обременительное ограничение на размер шага интегрирования [25]. Поэтому, явные методы применяются для решения задач средней жесткости [21].
Для решения жестких задач обычно используют L - устойчивые методы [5-20,56-76]. В случае большой размерности системы дифференциальных уравнений в L - устойчивых методах общие вычислительные затраты фактически полностью определяются временем вычисления и обращения матрицы Якоби. Затраты можно значительно уменьшить, если использовать одну и ту же матрицу Якоби на нескольких шагах интегрирования, то есть применять алгоритмы с замораживанием матрицы Якоби. В итерационных методах замораживание матрицы не влияет на порядок точности численной схемы, а определяет скорость сходимости итерационного процесса. Поэтому данный подход широко используется при реализации методов такого типа [26]. В безытерационных методах, к которым относятся методы типа Розенброка [77] и их различные модификации (например, [9], [32]), матрица Якоби включена непосредственно в численную формулу. Поэтому ее аппроксимация может приводить к потере точности численной схемы, что приводит к определенным проблемам с ее замораживанием.
