Исследование линейных многошаговых методов тема автореферата и диссертации по математике, 01.01.07 ВАК РФ

Введение

Глава I. Методические вопросы построения линейных многошаговых методов (ЛММ)

§ I.I. Порядок, устойчивость и сходимость ЛММ.

§ 1.2. Методы для решения жестких задач.

1.2.1. Явление жесткости.

1.2.2. Специфика понятия устойчивости при решении жестких задач

1.2.3. Особенности ЛММ для жестких задач.

§ 1.3. Методы в форме Нордсика и их эквивалентность

Глава 2. Построение методов универсального характера, ориентированных на решение жестких и нежестких задач.

§ 2.1. Принципы построения.

§ 2.2. Построение к-шаговых методов порядка к

2.2.1. а -устойчивые методы 2-го порядка

2.2.2. Жестко устойчивые методы 3-го порядка

2.2.3. Жестко устойчивые методы 4, 5 и 6-го порядков

§ 2.3. Построение /с-шаговых методов порядка k-i

Глава 3. Алгоритмизация построенных ЛММ и их тестирование

§ 3.1. Программная реализация построенных методов с использованием алгоритма выбора шага и порядка

§ 3.2. Численное решение смешанных дифференциальноалгебраических систем вида F (у' с[, £) = О

5F дг с разреженными матрицами

3.2.1. Алгоритм решения дифференциально-алгебраических систем.

3.2.2. Особенности решения больших систем

F с разнеженными матрицами частных производных и -jp

§ 3.3. Принципы тестирования ЛММ . III

§ 3.4. Результаты численных экспериментов

Глава 4. Математическое моделирование нестационарной радиационной электропроводности полимеров.

§ 4.1. Физико-математическая модель Роуза-Фаулера

Вайсберга (РФВ)

§ 4.2. Численный алгоритм решения уравнений модели и его тестирование

§ 4.3. Анализ температурной зависимости радиационной электропроводности и сравнение с экспериментом

Настоящая работа посвящена исследованию линейных многошаговых методов (ЛММ) решения задачи Коши для обыкновенных дифференциальных уравнений (ОДУ).

В этот класс входят методы Адамса, которые, наряду с методами Рунге-Кутты, получили наиболее широкое распространение для численного интегрирования ОДУ с использованием ЭВМ. Однако в 50-х годах было отмечено [59], что эти традиционные методы пригодны не для всех задач. Для некоторых типов ОДУ, названных жесткими, они требуют шага, величина которого очень мала по сравнению с длиной отрезка интегрирования, несмотря на достаточно медленное изменение искомых функций. Такие уравнения описывают физические модели с сильно различающимися временными постоянными. Источники их возникновения и свойства подробно описаны в [35].

С расширением множества решаемых задач жесткие системы ОДУ все чаще стали встречаться в научных исследованиях, что потребовало разработки для них эффективных методов, удовлетворяющих сформулированным в [62,75,132] специальным условиям устойчивости. Однако, как оказалось, наилучшие в этом смысле методы значительно уступают непригодным для жестких задач традиционным методам по тем овойствам точности и устойчивости, которые определяют эффективность метода при решении нежестких задач, а часто и по вычислительным затратам на один шаг интегрирования (что видно из сопоставления одношаговых методов для жестких задач - неявных методов Рунге-Кутты [51] или семейства методов Розенброка [123], представители которого получены, например, в [2,24,29,55,123], с явными методами Рунге-Кутты для нежестких задач). Поэтому для решения жестких и нежестких задач обычно применяются различные группы методов, чем и объясняется тот факт, что современные программы численного интегрирования ОДУ либо содержат одну из этих групп и предназначены для задач одного типа - жестких или нежестких (см., например, основанные на неявных одношаговых методах алгоритмы для жестких задач в [3,4,6,7,19,20,32]), либо являются универсальными за счет включения разных групп методов, но, по сути, реализуют два автономных алгоритма, каждый из которых ориентирован на решение задач одного типа.

Ярким примером описанной ситуации служит также класс ЛШ. Наиболее эффективными для нежестких задач являются методы Адам-са, непригодные в жестком случае. С другой стороны, специальным требованиям устойчивости при решении жестких задач наилучшим образом удовлетворяют формулы дифференцирования назад 1-6 порядков (для 1-го и 2-го порядков они впервые предложены в [59]), недостаточно эффективные в нежестком случае. Поэтому программы, основанные на методах Адамса, например, DVDQ [98], Ш5[124 J)E [125], предназначены только для нежестких задач, на формулах дифференцирования назад, например, программа SOLVER [22], учитывающая разреженность якобиана, только для жестких задач, а универсальные программы, например, DIFSUB [77], GEAR [84], STIFF [94,26], EPISODE [b2], включают методы Адамса и формулы дифференцирования назад (называемые также методами Гира).

Алгоритмы всех указанных программ базируются на теоретическом анализе ЛШ общего вида с постоянным шагом [9,18,38,44, 60,61,76,83,106], некоторые аспекты которого будут подробно рассмотрены в первой главе, а также на многочисленных исследованиях (например, [28,40,48,73,74,79-82,92,97,100,102,107] ), посвященных вопросам реализации ЛШ и касающихся преимущественно методов Адамса и формул дифференцирования назад.

В ряде работ изучены разные варианты предсказывающе-ис-правляющего процесса, которые используются для вычислений на каждом шаге интегрирования и включают явный метод, как правило, метод Адамса-Бэшфорта, в качестве предиктора и неявный метод, чаще всего метод Адамеа-Мултона или формулу дифференцирования назад, в качестве корректора. Существуют два подхода к реализации этого процесса, что связано с необходимостью применять итерационный метод для решения корректора, нелинейного относительно значения искомой функции в текущей точке. Первый подход, называемый исправлениями до сходимости, состоит в том, что сходимость достигается за счет проведения нужного числа итераций или за счет уменьшения шага интегрирования при ограниченном числе итераций. Другой подход, возможный при условии, что корректор решается методом простой итерации, заключается в применении фиксированного числа итераций, что тактически приводит к явному алгоритму. Для изложения этого подхода обычно пользуются введенными в [88] обозначениями: Р - применение предиктора, С - одна итерация корректора, Е - вычисление правой части системы ОДУ.

В [76] исследованы оба подхода для ЛММ общего вида и показано, что при исправлениях до сходимости устойчивость определяется корректором и не зависит от предиктора, а при втором подходе на абсолютную устойчивость, т.е. устойчивость при фиксированном шаге h для модельной линейной задачи у'- J(/ а - комплексная постоянная), влияют и предиктор и корректор. В [104] изучена устойчивость разных алгоритмов, реализующих второй подход, и доказано, что для заданного корректора максимальная область абсолютной устойчивости (совокупность точек плоскости их , в которых метод абсолютно устойчив) для алгоритмов Р(ЕС)Мт может быть больше, чем для Р^С)Н£, при любых предикторах и любом числе итераций М • Устойчивости разных вариантов предсказывающе-исправляющего алгоритма, основанного на методах Адамса, посвящена работа [82]. При сравнении их областей абсолютной устойчивости автор признает наилучшим алгоритм РкЕСк+1Е ( к означает порядок метода), который при к>5 имеет лучшие свойства устойчивости, чем исправления до сходимости [82]. Однако, т.к. при £ разница в устойчивости в этих двух случаях незначительна, то для практических вычислений рекомендуются исправления до сходимости по причине их экономичности.

Выло предпринято несколько попыток увеличить за счет выбора предиктора область абсолютной устойчивости алгоритмов с корректором Адамеа-Мултона, реализующих второй подход. В [58, 96 ] получены предикторы для алгоритма РЕ СЕ, а в [94] - для РЕ С, В [127] предложено расширить область абсолютной устойчивости за счет последнего вычисления производной как линейной комбинации ее значений, полученных на предыдущих итерациях. Однако вое эти усовершенствования позволили добиться лишь незначительного улучшения.

Поскольку независимо от конкретного варианта реализации и вида предиктора и корректора процесс, организованный в рамках второго подхода, является явным, он непригоден для решения жестких задач так же, как и исправления до сходимости простыми итерациями, сильно ограничивающими величину шага интегрирования. В жестком случае для решения корректора обычно применяется метод Ньютона или его модификации, использующие якобиан системы ОДУ, что впервые было предложено в [74]. Получаемая при этом погрешность изучена в [109].

Трудности, вызванные неявностью ЛММ для жестких задач, привели к поиску путей повышения устойчивости явных методов с помощью якобиана. На его использовании основаны обобщенные многошаговые методы [Ю5], которые являются "линейно неявными" в отличие от полностью неявных. Подмножество этих методов изучено в [87]. К обобщенным многошаговым методам относятся методы, описанные в (49,IOS], а также модификации явных методов Адамса-Бэшфорта [93,118], пригодные для решения жестких задач (отметим, что значительно раньше якобиан был введен Розенброком [123] в методы Рунге-Кутты).

Другую модификацию ЛММ для жестких задач представляют методы, использующие вторую производную искомой функции. Одно- и двухшаговые методы этого типа предложены соответственно в [112] и [91]. Они входят в семейство разработанных в [бб] к -шаговых методов порядка к + 2 , к^ /2,7 , которые хотя и являются более точными и устойчивыми, чем формулы дифференцирования назад, но требуют большого количества вычислений якобиана, необходимых для определения второй производной искомой функции, что делает их слишком трудоемкими.

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

В [Ю0] обсуждаются 10 различных способов изменения шага, наиболее эффективными из которых признаны интерполяционный вариант, основанный на методе с постоянным шагом и на получении информации в новых точках при изменении шага с помощью интерполяции через известные значения, и вариант с переменным шагом, использующий методы, выведенные непосредственно для неравномерной сетки. Было отмечено (например, в [48]), что эти варианты не эквивалентны.

Интерполяционный вариант для методов Адамса при сохранении информации в виде разностей назад исследован в [I00J, где показано, что уменьшение шага вдвое при высоком порядке метода может вызвать увеличение погрешности. Это является недостатком интерполяционного варианта и не зависит от вида сохраняемой информации. В [ЮО] предложен также эффективный алгоритм удвоения шага и деления его пополам с использованием в последнем случае сглаживания разностей.

Упростить процесс интерполяции при произвольном изменении величины шага позволяют методы, предложенные Нордсиком [П7]. Они эквивалентны методам Адамеа-Мултона [120], от которых отличаются только видом сохраняемой информации. Методы в представлении Нордсика, эквивалентные формулам дифференцирования назад, получены Гиром [73] и использованы в программах DIFS^S [77], QEAR [84], ЭТТ/Т[94,26] . Реализованный в этих программах алгоритм с автоматическим выбором величины шага и порядка метода допускает естественное обобщение (описанное Гиром в [78]) на дифференциально-алгебраические системы, неразрешенные относительно производных. Вопросы существования и единственности решения таких систем обсуждаются в [1,10,72,114,115] : в [114,115] формулируются достаточные для этого условия, в [1,10] рассматриваются частные случаи сингулярных систем, а в [72] поведение решения изучается на конкретных примерах. В [114,115] и в [ill] предлагаются различные способы обобщения ЛММ на эти системы.

Вариант с переменным шагом потребовал обобщения ЛММ на неравномерную сетку. Методы с переменным шагом, эквивалентные в случае постоянного шага методам Адамса, разработаны в [39], а также в [121], где доказана их устойчивость и сходимость при легко выполняемых требованиях на стратегию выбора шага. Дальнейшее обобщение методов Адамса [*97,102] привело к созданию алгоритма РЕ СЕ для решения дифференциальных уравнений, имеющих порядок выше первого. Методы с переменным шагом при сохранении информации в виде разностей назад, эквивалентные при постоянном шаге формулам дифференцирования назад, получены в [48]. В [52] предложено обобщение на неравномерную сетку эквивалентных формулам дифференцирования назад методов в представлении Нордоика и дано описание созданной на этой основе программы EPISODE. Ее модификация, позволяющая уменьшить количество вычислений якобиана системы ОДУ и LU -факторизаций матриц, представлена в [92]. В [23,47] проведен анализ абсолютной устойчивости формул дифференцирования назад 2-го и 3-го порядков на неравномерной сетке, шаги которой образуют геометрическую прогрессию.

Отметим, что обобщение на неравномерную сетку проведено не для всех ЛММ, а только для методов, явно использующих идею интерполяции, и, в частности, для методов Адамса и формул дифференцирования назад. В [ИЗ] сделана попытка построить на неравномерной сетке многошаговые методы с применением интерполяции на основе обобщенных многочленов. При этом шаблон, т.е. множество точек, в которых задаются значения искомой функции и ее производной, являетоя фиксированным и определяет порядок аппроксимации. В этот класс попадают методы Адамса и формулы дифференцирования назад на неравномерной сетке, но при постоянном шаге он является более узким, чем класс ЛШ общего вида.

Оба способа изменения шага - интерполяционный и с переменным шагом, а также изменение порядка исследуются в [79-81], где для предсказывающе-исправляющих алгоритмов с фиксированным числом исправлений и с исправлениями до сходимости, использующих методы Адамеа, получены следующие результаты: вариант с переменным шагом сходится при любой схеме выбора шага и порядка; интерполяционный вариант сходится, еоли схема выбора шага и порядка гарантирует после каждого изменения их постоянство по крайней мере на к шагах, где к - порядок метода. В этих работах рассматриваются также алгоритмы на основе устойчивых аппроксимирующих ЛММ общего вида, для которых доказано, что при малой скорости изменения шага оба варианта дают устойчивый алгоритм, а если предусмотрена смена методов в процессе счета, то алгоритм устойчив при не зависящих от ОДУ и величины шага ограничениях на частоту изменений.

Кроме изучения теоретических свойств разных алгоритмов численного интегрирования ОДУ с помощью ЛММ проводилось их широкое тестирование с целью сравнения между собой и с алгоритмами, базирующимися на других методах. При этом в разных работах выдвигаются различные принципы тестирования, соблюдение которых, по мнению авторов, должно обеспечить правильность выводов о качестве проверяемых алгоритмов. В [99,101] , например, перечисляются требования к получению и представлению результатов тестирования, цель которого состоит в том, чтобы испытать алгоритм в различных предельных ситуациях и изучить его характеристики как при этих, так и при промежуточных условиях.

Этот подход демонстрируется на проверке программы с переменным шагом и порядком, использующей методы Адамса. Другая идея реализована в [67,89] на примере сравнения методов для нежестких задач, автоматизированного в результате разработки программы DETEST, которая решает набор систем ОДУ" при разных условиях точности и для каждого варианта выдает значения нескольких параметров, используемых для оценки надежности метода. Программа включает алгоритм выбора шага на основе оценки приращения глобальной погрешности на шаге, отнесенной к единичному шагу интегрирования, что обеспечивает при достаточно высокой заданной точности решение задач с почти одинаковой глобальной погрешностью при использовании различных методов, требующих разного объема вычислительной работы, который и входит в оценку качества метода. Результаты тестирования по различным методикам (ом. также [69,126] ) показали, что при решении нежестких задач с достаточно высокой заданной точностью алгоритмы с переменным шагом и порядком, основанные на методах Адамса, являются более эффективными, чем алгоритмы, базирующиеся на других методах, особенно в тех случаях, ковда вычисление правой части системы ОДУ требует больших затрат.

В [64,68-71] предприняты попытки применить подход, описанный в [89], к сравнению методов для жестких задач. В [70, 71] подбор и классификация тестов были направлены на выявление зависимости результатов от математических свойств ОДУ. В [68] в качестве тестов взяты жеоткие задачи из области химии. Проведенные расчеты позволили сделать вывод о хорошем поведении программ, использующих формулы дифференцирования назад.

Наличие нескольких таких программ с разными алгоритмами изменения шага потребовало их сравнения, которое и было выполнено в [68], где рассматривались программы DIF$UB, GEAR, EPI&0DE, а также в [53,54,85] , где проверялись только две последние из них, т.к. DIFSVB уступает1 GEAR , являющейся ее удачной модификацией и реализующей тот же интерполяционный вариант изменения шага, но с более мягким контролем погрешности. Тестирование показало, что для большинства жестких задач средней сложности программа GEAR по эффективности превосходит EPISODE, т.к. последняя требует большого количества вычислений якобиана и LU-факторизаций матрицы. Однако, т.к. алгоритм, реализованный в EPISODE, обладает лучшими свойствами устойчивости, то он может быть применен к некоторым задачам, с которыми GEAR по справляется.

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

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

Для методов общего назначения, ориентированных на решение жестких и нежестких задач, в работе сформулирован критерий оптимизации, в котором требование аппроксимации сочетается с условиями устойчивости при малых и больших шагах интегрирования. Получена группа отвечающих этому критерию ЛММ 2-6 порядков. Для анализа свойств ЛММ введена новая характеристика -область аппроксимации (на одном шаге интегрирования) экспоненты - решения модельной задачи y'-Jty в комплексной плоскости hA . Размер и форма таких областей для разных уровней точности относятся к важнейшим характеристикам ЛММ. Данный подход позволил уточнить определение жесткой устойчивости [75], также включающее требование точности при конечном шаге, которое однако в [75] не конкретизировано.

Для ЛММ общего вида доказано взаимно однозначное соответствие k -шаговых методов порядка к и методов Нордсика того же порядка и установлена связь между их характеристиками устойчивости. Получена модификация представления Нордсика для к -шаговых методов порядка к- 1 (поскольку к ним относится часть оптимизированных методов). Отмечена экономичность по памяти представления Нордсика для ^-шаговых методов общего вида.

На основе программы STIFF [94,26] разработаны две программы: программа для решения систем ОДУ, разрешенных относительно производных (будем называть ее модифицированным вариантом $TIFF), и программа STIFSP, которая реализует обобщенный на ЛММ произвольного вида алгоритм Гира [78] и предназначена для решения дифференциально-алгебраических систем большой размерности, неразрешенных относительно производных, с учетом разреженности якобиана. Обе программы включают три группы методов в представлении Нордсика - Адамса для нежестких задач, Гира и построенные оптимизированные методы для задач, жестких на всем отрезке интегрирования или на отдельных его участках (в том числе, для систем, содержащих алгебраические уравнения и представляющих предельный случай жестких задач). В программе STIF2P для решения возникающих в процессе счета систем линейных алгебраических уравнений с разреженной матрицей применяется адаптация программы А/шт.

Сравнение новых методов с методами Адамса и Гйра проведено в соответствии с обоснованными в работе общими принципами тестирования, касающимися, в частности, выбора в качестве тестов задач с известным точным решением и заданными математическими свойствами. Это позволило определить ограничения, возникающие при решении жестких задач на ЭВМ типа БЭСМ-6.

С помощью программы STIFSP (а именно той ее части, которая реализует построенные оптимизированные методы) проведено полномасштабное параметрическое численное моделирование нестационарной радиационной электропроводности полимеров в рамках модели Роуза-Фаулера-Вайсберга [п] . В результате расчетов определены: эффективная энергия активации в широком диапазоне температур; условия перехода от аномального переноса заряда к нормальному диффузионному; эффективная дрейфовая подвижность электронов в условиях однородной (по объему) генерации в широком временном интервале; нелинейная зависимость установившейся подвижности от мощности генерации; поправочные множители, позволяющие значительно улучшить согласие аналитических формул с результатами экспериментов.

Изложим далее содержание диссертации по главам.

В первой главе рассматриваются методические вопрооы построения ЛММ с постоянным шагом Н : предназначенных для решения задачи Коши: y'-mt), t4To,T,) у(То)-у0^ где У,У'>Уо,/ - векторы размерности s .

В § I.I описаны важнейшие характеристики ЛММ, используемые при их теоретическом анализе, - порядок, устойчивость, сходимость. Обсуждаются разные аспекты устойчивости и их связь с многочленом устойчивости Р (f) f hA с (f) , где к , к J pCf)= £ ^ f , fif ~L . На основе модельной задачи вводится критерий точности ЛММ при конечном h . Размер и форма областей абсолютной устойчивости и точности (при различных заданных ее уровнях) на комплексной плоскости hH вместе с величиной корней многочлена устойчивости при разных hA входят в число основных характеристик метода.

В § 1.2 обсуждается понятие жесткости и отвечающие ему определения устойчивости [62,75,132]. Затем рассматриваются методы Гира, наиболее широко используемые для решения жестких задач.

В § 1.3 исследуется представление Нордсика для ЛММ, которое превращает многошаговый метод фактически в одношаговый с запоминанием (в случае скалярной задачи) вектора информации С у», hyn/ tf(/n"/2!f. v h купк)/Ь /) r, где y(nL) - величины, приближающие производные решения yil)(tnJ. Доказывается взаимно однозначное соответствие к -шаговых методов порядка к и методов Нордсика того же порядка. Устанавливается связь между матрицей перехода да вектора вариаций (<Гу„/ф")г .{(^'))т которая используется для исследования устойчивости методов Нордсика, и многочленом устойчивости эквивалентных им ЛММ. Предлагается модификация представления Нордсика для к -шаговых методов порядка к- 1 . Обсуждаются преимущества и недостатки представления Нордеика при его использовании в алгоритмах с переменным шагом и порядком.

Вторая глава посвящена построению методов, ориентированных на решение жестких и нежестких задач.

В § 2.1 формулируются принципы их построения, объединяющие аппроксимацию при малых hA с устойчивостью при малых и больших hA для модельной задачи у/=Лу . Требование аппроксимации состоит в том, чтобы один из корней многочлена устойчивости, называемый главным, приближал с заданным порядком по h точное решение модельной задачи G ^ . Требование устойчивости заключается в минимизации всех корней многочлена устойчивости при Re т.е. корней &(f) , и его посторонних корней (всех кроме главного) при hA-*0 , т.е. корней многочлена (f)- J>(f) /(1-f'). При этом рассматривается две возможности оптимизации: минимизировать максимум модулей корней за счет понижения порядка на один в методах Адамса-Мултона и минимизировать в к -шаговых методах порядка к максимум модулей корней <JffJ и fi(f)* Дополнительную минимизацию предлагается провести за счет дальнейшего понижения порядка, что является оправданным, поскольку во многих практических расчетах, особенно при невысокой заданной точности, условия устойчивости накладывают более сильные ограничения на максимальный шаг интегрирования, чем условие точности, и поэтому методы высоких порядков не требуютоя (и не используются, даже если они включены в программу с автоматическим выбором порядка). В таких случаях методы пониженного порядка, но с улучшенными свойствами устойчивости позволят повысить эффективность расчетов.

В § 2.2 описано получение к -шаговых методов порядка к .

Методы 2-го и 3-го порядков исследуются аналитически, а для наховдения методов более высоких порядков применяется численный алгоритм направленного случайного поиска. При реализации первой возможности получены методы типа Адамса: жестко устойчивые 2-го и 3-го порядков, а также 4-го порядка, который лишь условно можно отнести к жестко устойчивым; показано, что при к>4 не существует -шаговых методов типа Адамса порядка к , пригодных для решения жестких задач. При втором способе минимизации построена группа методов 2-6 порядков, условно названных "методами с равными корнями", поскольку в них максимумы модулей корней 6"Yfj и р1 (()оказались равными. Проводится сравнительный анализ свойств устойчивости и точности полученных методов, методов Адамса и Гира.

В § 2.3 на примере методов Гира и построенных методов типа Адамса 2-го и 3-го порядков изучаются возможности улучшения свойств устойчивости за счет расширения шаблона на одну точку (что эквивалентно понижению порядка на один при фиксированном шаблоне).

В третьей главе излагаются вопросы, связанные с алгоритмизацией ЛММ и их тестированием.

В § 3.1 описывается программная реализация построенных методов с использованием интерполяционного варианта изменения шага на основе программы STIFF [94,26] , в которой выбор величины шага интегрирования и порядка метода осуществляется по алгоритму, предложенному в [77] .

В § 3.2 алгоритм [78], использующий методы Гира для численного решения дифференциально-алгебраических систем

0 t обобщен на ЛШ произвольного вида. Рассматривается специфика решения таких систем большой размерности, с разреженными матрицами частных производных jpr > jp •

Она учтена в предназначенной для этого программе STIFSP, которая реализует три семейства методов в рамках обобщенного алгоритма (78] и основана на описанном в § 3.1 модифицированном варианте программы STIFF.

В § 3.3 формулируются принципы тестирования ЛММ, обосновывается выбор в качестве тестов задач не по предметным областям, а по математическим свойствам, причем требуется знание их точного решения.

В § 3.4 приводятся результаты проведенных в соответствии с предложенными принципами численных экспериментов, которые подтверждают эффективность построенных методов при решении задач, жестких как на всем отрезке интегрирования, так и на отдельных его участках. Указываются ограничения на достижимую глобальную точность и дополнительные затраты на вычисления, возникающие при решении жестких задач на ЭВМ типа БЭСМ-6 и выявленные с помощью тестирования задач с известным точным решением.

Четвертая глава посвящена численному моделированию нестационарной радиационной электропроводности (НРЭ) полимеров с использованием программы

StifSp.

В § 4.1 дана физико-математическая модель процесса, учитывающая рекомбинацию носителей заряда при наличии их захвата на энергетические уровни - ловушки (модель Роуза-Фаулера-Вай-сберга), уравнения которой преобразованы к виду, удобному для проведения расчетов.

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

В § 4.3 на основе численных расчетов проведен анализ температурной зависимости НРЭ в широком диапазоне температур для различных распределений ловушек по энергии, а также нелинейных эффектов, связанных с рекомбинацией и заполнением ловушек.

Основные результаты диссертации опубликованы в работах [12 - 17] .

ЗАКЛЮЧЕНИЕ

В заключение сформулируем основные результаты работы.

1. Предложен новый подход к построению ЛММ, в рамках которого разработаны две группы оптимизированных методов 2-6 порядков точности. Они обеспечивают эффективное решение систем ОДУ без предварительного анализа на жесткость и реализованы в виде двух стандартных программ на языке ФОРТРАН для решения дифференциально-алгебраических систем, разрешенных и неразрешенных относительно производных.

2. Сформулированы принципы тестирования ЛММ, на основе которых проведено сравнение построенных методов с известными ранее ЛММ.

3. Проведено многопараметрическое численное исследование температурной зависимости нестационарной радиационной электропроводности полимеров на основе полной модели Роуза-Фаулератт 7

Вайсберга на временном интервале 10 хх - 10' с. Показана неприменимость квазиравновесного приближения на этапе спада электропроводности. Математическое моделирование позволило согласовать асимптотические закономерности с экспериментальными данными.

1. Абрамов А.А., Балла К., Конюхова Н.Б. Перенос граничных условий из особых точек для систем обыкновенных дифференциальных уравнений. - М., ВЦ АН СССР, Сообщ. по вычислительной математике, 1981.

2. Артемьев С.С., Демидов Г.В. Численные методы высокого порядка точности для решения жестких систем. В кн.: Математические проблемы химии. Новосибирск, ВЦ СО АН СССР, 1975.

3. Артемьев С.С., Демидов Г.В. Алгоритм переменного порядка и шага для численного решения жестких систем обыкновенных дифференциальных уравнений. ДАН СССР, 238, В 3, 1978.

4. Артемьев С.С., Демидов Г.В. Алгоритм переменного порядка и шага для численного решения жестких систем обыкновенных дифференциальных уравнений. Новосибирск, ВЦ СО АН СССР, 1978, препринт № 45.

5. Архипов В.И., Колесников В.А., Руденко А.И. Дисперсионный транспорт носителей заряда в поликристаллических слоях пен-тацена. Изв. АН Латв. ССР, Сер. физ. и техн. наук, 1981, JS 6, с.10.

6. Банчев ВД. Программа численного интегрирования задачи Коши для жестких систем обыкновенных дифференциальных уравнений. Дубна, Сообщ. ОИЯИ, PII-9677, 1976.

7. Банчев В.Ц. Алгоритм решения задачи Коши для жестких систем обыкновенных дифференциальных уравнений. Дубна, Сообщ. ОШИ, PII-9678, 1976.

8. Банчев БД. Об обобщенных методах интегрирования обыкновенных дифференциальных уравнений. Дубна, Сообщ. ОИЯИ, P5-I2I74, 1979.

9. Бахвалов Н.С. Численные методы. М.: Наука, 1975.

10. Боярияцев Ю.Е. О структуре общего решения краевой задачи для сингулярных систем ОДУ. Новосибирск, ВЦ 00 АН СССР, 1978, препринт № 44.

11. Вайсберг С.Э. Обратимые радиационные эффекты в полимерах. -В кн.: Радиационная химия полимеров. Под ред. Каргина В.А., стр.376-443, М.: Наука, 1973.

12. Герасимов Б.П., Кульчицкая И.А. Об одном семействе линейных многошаговых методов для ОДУ. Препринт Ин. прикл. матем. им. М.В.Келдыша АН СССР, 1982, № 51.

13. Герасимов Б.П., Кульчицкая И.А. Программа интегрирования систем ОДУ, неразрешенных относительно производных. Отчет Ин. прикл. матем. им. М.В.Келдыша АН СССР, 1983>t Ш СИ098.

14. Герасимов Б.П., Кульчицкая И.А. Об одном подходе к построению линейных многошаговых методов. Препринт Ин. прикл. матем. им. М.В.Келдыша АН СССР, 1984, № 76.

15. Герасимов Б.П., Кульчицкая И.A. STIR5P- пакет программ интегрирования дифференциально-алгебраических систем большой размерности. Препринт Ин. прикл. матем. им. М.В.Келдыша АН СССР, 1984, № 103.

16. Герасимов Б.П., Кульчицкая И.А., Мингалеев Г.С., Тютнев А.П. Температурная зависимость нестационарной радиационной электропроводности полимеров: теоретический анализ и численные расчеты. Препринт Ин. прикл. матем. им. М.В.Келдыша АН СССР, 1984 , В Г46.

17. Горбунов А.Д. Разностные методы решения задачи Коши для системы обыкновенных дифференциальных уравнений. М.: Изд-во МГУ, 1973.

18. Демидов Г.В., Лаевский Ю.М. Численно-аналитический метод решения жестких систем обыкновенных дифференциальных уравнений. Новосибирск, ВЦ СО АН СССР, 1979, препринт № 184.

19. Демидов Г.В., Новиков В.А. Программа STEK (модификация программы МERS0N). Новосибирск, Щ СО АН СССР, 1981, № 313.

20. Демидов Г.В., Юматова Л.А. Исследование точности неявных од-ношаговых методов. Новосибирск, Щ СО АН СССР, 1976, препринт № 25.

21. Еремин А.Ю., Марьяшкин Н.Я. Пакет программ SOLVER системы нелинейных функциональных и обыкновенных дифференциальных уравнений с разреженными якобиевыми матрицами. М., ВД АН СССР, Сообщ. по программному обеспечению ЭВМ, 1980.

22. Заворин А.Н. Влияние изменения шага на устойчивость некоторых формул дифференцирования назад. IBM и МФ, 21, № 6,1981.

23. Заворин А.Н., Хесина И.Я. О некоторых численных методах решения жестких систем обыкновенных дифференциальных уравнений. IBM и МФ, 13, № I, 1973.

24. Захаров А.Ю. Некоторые результаты сравнения эффективности решения систем обыкновенных дифференциальных уравнений. -Препринт Ин. прикл. матем. им. М.В.Келдыша АН СССР, 1979, № 125.

25. Захаров А.Ю., Турчанинов В.И. S7TFF программа для решения жестких систем обыкновенных дифференциальных уравнений (адаптация для ЭВМ БЭСМ-6). - Инструкция Ин. прикл. матем. им. М.В.Келдыша АН СССР, 1977.

26. Звягин И.П. Кинетические явления в неупорядоченных полупроводниках. М.: Изд-во МГУ, 1984.

27. Ионеску Д.В. Остаточный член в формуле численного интегрирования Адамса. ЖВМ и МФ, 2, № I, 1962.

28. Калиткин Н.Н., Кузьмина Л.В. Интегрирование жестких систем дифференциальных уравнений. Препринт Ин. прикл. матем. им. М.В.Келдыша АН СССР, 1981, № 80.

29. Корн Г., Корн Г. Справочник по математике для инженеров и научных работников. М.: Наука, 1977.

30. Курош А.Г. Курс высшей алгебры. М.: Наука, 1965.

31. Новиков Е.А. Программа для решения жестких систем обыкновенных дифференциальных уравнений. Алгоритмы и программы, ВНГИ центр, » 3 (47), 1982.

32. Пикаев А.К., Кабакчи С.А., Макаров И.Е., Ершов Б.Г. Импульсный радиолиз и его применение. М.: Атомиздат, 1979.

33. Постников М.М. Устойчивые многочлены. М.: Наука, 1981.

34. Ракитский Ю.В., Устинов С.М., Черноруцкий И.Г. Численные методы решения жестких систем. М.: Наука, 1979.

35. Самарский А.А. Теория разностных схем. М.: Наука, 1983.

36. Соболь И.М., Левитан Ю.Л. Получение точек, равномерно расположенных в многомерном кубе. Препринт Ия. прикл. матем. им. М.В.Келдыша АН СССР, 1976, № 40.

37. Современные численные методы решения обыкновенных дифференциальных уравнений. М.: Мир, 1979.

38. Стукало А.С. Общий метод Адамса с переменным шагом интегрирования. В сб.: Вычислительная математика, вып.2, изд-во Киевского Ун-та, 1966.

39. Тихонов А.Н., Горбунов А.Д. Асимптотические разложения погрешности разностного метода решения задачи Коши для систем обыкновенных дифференциальных уравнений. ЖВМ и МФ, 2, № 4, 1962.

40. Тютнев А.П., Саенко B.C., Пожидаев Е.Д. Радиационная электропроводность полимеров при импульсном облучении. ДАН СССР, 266, Ш I, 168-172, 1982.

41. Тютнев А.П., Саенко B.C., Дунаев А.Ф., Сичкарь В.П., Ванников А.В. Температурная зависимость нестационарной радиационной электропроводности полимеров. ДАН СССР, 277, $ 2, 424-429, 1984.

42. Широбоков Н.В. К определению жестких дифференциальных задач. ЖВМ и МФ, 24, Л 4, 1984.

43. Штеттер X. Анализ методов дискретизации для обыкновенных дифференциальных уравнений. М.: Мир, 1978.

44. Axelsson 0. A note on a class of strongly A-stable methods. BIT, 12, 1-4, 1972.

45. Brayton R.K., Gustavson E.G., Hachtel G.D. A new efficient algorithm for solving differential algebraic systems using implicit backward differentiation formulas. Proc. IEEE, 60, 98-108, 1972.

46. Brunner H. Marginal stability and stabilization in thenumerical integration of ODE. Math, of сотр., 24, N 111, 1970.

47. Brunner H. A class of A-stable two-step methods based on Schur polynomials. BIT, 12, 468-474, 1972.

48. Butcher J.C. Implicit Runge-Kutta processes. Math. Сотр., 18, 50-64.

49. Byrne G.D., Hindmarsh A.C. A polyalgorithm for the numerical solution of ODE. ACM Transactions on mathematical software, 1, N 1, 71-96, 1975.53» Byrne G.D., Hindmarsh A.C., Jackson K.R., Brown H.G.

50. Comparative test results for two ODE solvers EPISODE and GEAR. - ANL-77-19, Argonne national lab., Argonne, Illinois, 1977.

51. Dahlquist G. Convergence and stability in the numerical integration of ordinary differential equations. Math. Scand., 4, 33-53, 1956.

52. Dahlquist G. Stability and error bounds in the numerical integration of ordinary differential equations. Trans, of the Royal Inst, of Techn. Stokholm, 130, 1959.

53. Dahlquist G. A spacial stability problem for linear multi-step methods. BIT, 3, 27-43,, 1963.

54. Dill C., Gear C.W. A graphical search for stiffly stable methods for ODE. J.ACM, 18, 75-79, 1971.

55. Ehle B.L. A comparison of numerical methods for solving certain stiff ordinary differential equations. Tech. Rept. N 70, Dept. of Math. Univ. of Victoria, 1972.

56. Eisenstat S.C., Gursky M.C., Schultz M.H., Sherman A.H. Yale Sparse, Matrix Package II. The Nonsymmetric Codes. Res. Rept. N 114, Dept. of Comput. Sci., Yale Univ.

57. Enright W.H. Second derivative multistep methods for stiff ordinary differential equations. SIAM J. Numer. Anal., 11, 321-331, 1974.

58. Enright W.H., Bedet R., Parkas I., Hull Т.Е. Test resultson initial value methods for non-stiff ordinary differential equations. Tech. Rept. N 68, Dept. of Comput. Sci., Univ. Toronto, 1973.

59. Enright W.H., Hull Т.Е. Comparing numerical methods for the solution of stiff systems of ODE's arising in chemistry. -In: Numerical methods for Differential Systems (L.Lapidus, W.E. Schiesser ed.), Academic Press, N-Y., 45-66, 1976.

60. Enright W.H., Hull Т.Е. Test results on initial value methods for nonstiff ordinary differential equations. -SIAM J. Numer. Anal., 13, 944-961, 1976.

61. Enright W.H., Hull Т.Е., Lindberg B. Comparing numerical methods for stiff systems of O.D.E.'s Tech. Rept. N 69, Dept. of Comput. Sci., Univ. of Toronto, 1974.

62. Enright W.H., Hull Т.Е., Lindberg B. Comparing numerical methods for stiff systems of ODE's. BIT, 15, N 1, 10-48, 1975.

63. Pox L., Mayers D.F. On the numerical solution of implicit ordinary differential equations. IMA J. of numerical analysis, 1, 377-401, 1981.

64. Gear C.W. The numerical integration of Ordinary Differential Equations. Math. Сотр., 21, 146, 1967.

65. Gear C.W. Numerical initial value problems in ordinary differential equations. Englewood Cliffs: Prentice Hall, 1971•

66. Gear C.W. Algorithm 407, DIFSTJB for solution of ordinary differential equations. Communs. Ass. comput. Mach., 14, 185-190, 1971.

67. Gear C.W., Tu K.W. Watanabe D.S. The stability of automatic programs for numerical problems. In: Stiff differential systems (R.A. Willoughby ed.) Plenum Press, N.-Y., 1974» 111-121.

68. Gear C.W., Watanabe D.S. Stability and convergence of variable order multistep methods. SIAM J.Numer. Anal., 11, 1044-1058, 1974.

69. Hall G. Stability analysis of predictor-corrector algorithms of Adams type. SIAM J.Numer. Anal., 494-505, 1974.

70. Henrici P. Discrete variable methods in ordinary differential equations. Wiley, 1962.

71. Hindmarsh A.C. GEAR: ordinary differential equations system solver, 1972, Lawrence Livermore Lab., Rept. UCID-30001, Rev. 3, Univ. of California, Livermore.

72. Hull Т.Е., Creemer A.L. Efficiency of predictor-corrector schemes. J. Ass. comput. Mach., 10, 291-301, 1963.

73. Hull Т.Е., Enright W.H., Fellen B.M., Sedgwick A.E. Comparing numerical methods for ordinary differential equations. SIAM J. Numer. Anal., 9, 603-637, 1972.

74. Isaacson E., Keller H.B. Analysis of numerical methods. -Wiley, 1966.

75. Jackson L.W., Kenue S.K. A fourth order exponentially fitted method. SIAM J.Numer. Anal., 11, 965-978, 1974.

76. Jackson K.R., Sacks-Davis R. An alternative implementationof variable stepsize multistep formulas for stiff ODE's.

77. Tech. Rept. N 121 (1978), Dept. of Computer Sci., Univ. of Toronto.93« Jain R.K. Some A-stable methods for stiff ordinary differential equations. Math. Сотр., 26, 71-77, 1972.

78. Kahaner D., Sutherland C.D. STIFF- a subroutine package for solving ordinary differential equations with prescribed initial conditions. Los-Alamos Scientific Lab., D-205, 1975»

79. Klopfenstein R.W., Millman R.S. Numerical stability of a one-evaluation predictor-corrector algorithm for numerical solution of ordinary differential equations. Math. Сотр., 22. 557-564, 1968.

80. Krogh P.T. VODQ/SVDQ/DVDQ variable order integrators for the numerical solution of ordinary differential equations. - TTJ Doc. N CP-2308, NPO-11643, Jet Propulsion Lab., Pasadena, Calif., 1969.

81. Krogh P.T. Opinions on matters connected with the evaluation of programs and methods for integrating ordinary differential equations. SIGNUM Newsletter, 7, 1972.

82. Krogh P.T. Algorithms for changing the stepsize. SIAM J. Numer.Anal., 10, 949-965, 1973.

83. Krogh P.T. On testing a subroutine for the numerical integration of ordinary differential equations. J.Ass. Comput. Mach., 20, 545-562, 1973.

84. Krogh P.T. Changing stepsize in the integration of differential equations using modified divided differences. In: Lecture Notes in Mathematics, 362, Ni-Y.: Springer, 22-71, 1974.

85. Lambert J.D. Linear multistep methods with mildly varying coefficients. Math, of сотр., 24, N 109, 1970.104» Lambert J.D. Predictor-corrector methods with identical regions of stability. SIAM J. Numer. Anal., 8, 337-344, 1971.

86. Lambert J.D., Sigurdsson S.T. Multistep methods with variable matrix coefficients. SIAM J.Numer. Anal., 9, 715-733, 1972.

87. Lambert J.D. Computational methods in ordinary differential equations. Wiley, 1973.

88. Liniger W. Global accuracy and A-stability of one and two-step integration formulae for stiff ordinary differential equations. In: Proc. Conf. on Numerical Solution of Differential Equations, Dundee 1968, Berlin, Springer, 188-193.

89. Liniger W. A criterion for A-stability of linear multistep integration formulae. Computing, 3, 280-285, 1968.

90. Liniger W. A stopping criteria for the Newton-Raphson method in implicit multistep integration algortithms for nonlinear systems of ordinary differential equations. -Communs. Ass. Comput. Mach., 14, 600-601, 1971.

91. W.Liniger. Multistep and one-leg methods for implicit mixed differential algebraic systems. IEEE Transactions on circuit and systems, CAS-26, N 9, 755-762, 1979.

92. Liniger W., Willoughby R.A. Efficient integration methods for stiff systems of ordinary differential equations. -SIAM J. Numer. Anal., 7, 47-66, 1970.

93. Ntfrsett S.P. An A-stable modification of the Adams-Bashforth methods. In: Proc. Conf. on Numerical Solution of Differential Equations, Dundee 1968. - Berlin, Springer, 214-219.

94. N^rsett S.P. A criterion for A( «^-stability of linear multi-step methods. BIT, 9, N 3, 259-263, 1969.

95. Osborne M.R. On Nordsieck's method for the numerical solution of ordinary differential equations. BIT, 6, 51-57, 1966.

96. Sedgwick A. An effective variable order variable step Adams method. Tech. Rept. N 53, Dept. of Comput. Sci., Univ. of Toronto, 1973.

97. Shampine L.F., Gordon M.K. Computer solution of ordinary-differential equations. W.H. Freeman and Co., 1975.

98. Tyutnev A.P. et.al. Radiation-induced conductivity in polymers under contineous irradiation. Phisica Status Solidi (a), 83, 365-373, 1984.

99. Tyutnev A.P., Saeriko V.S., Pozhidaev E.D., Akkerman A.F. Transient radiation conductivity in polymers. Physica status solidi (a), 73, N 1, 81-83, 1982.

100. Tyutnev A.P., Saenko V.S., Vannikov A.V. Mingaleev G.S. Radiation-induced conductivity in polyethelene. Phisica status solidi (a), 78, N 2, 689, 1983.

101. Tyutnev A.P., Saenko V.S., Dunaev A.F. et al. Temperature dependence of transient radiation induced conductivity in polymers. Physica status solidi (a), 85, N 1, 679, 1984.

102. Widlund O.B. A note on -unconditionally stable linear multi-step methods. BIT, 65-70, 1967.

103. Yakoviev b.s., Novikov g.f. Кинетика рекомбинации электронов и ионов, генерированных ионизирующим излучением в твердых углеводородах. Int. J.Radiat. Phys. Chem., 7, N 6, p.679, 1975.