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

Р Г о 00

- 1 ПАП 1093

РОССИЙСКАЯ АКАДЕМИЯ НАУК СИБИРСКОЕ ОТДЕЛЕНИЕ ИНСТИТУТ МАТЕМАТИКИ

на правах рукописи

Фадеев Станислав Иванович

УДК 519.615: 519.624

ОРГАНИЗАЦИЯ ЧИСЛЕННОГО ЭКСПЕРИМЕНТА ПРИ ИССЛЕДОВАНИИ НЕЛИНЕЙНЫХ КРАЕВЫХ ЗАДАЧ МЕТОДОМ ПРОДОЛЖЕНИЯ РЕШЕНИЯ ПО ПАРАМЕТРУ

Специальность 05.13.16 - применение вычислительной техники, математического моделирования и математических методов в научных исследованиях

Диссертация на соискание ученой степени доктора физико-математических наук (в форме научного доклада)

Новосибирск 1992

Работа выполнена в Институте математики СО РЯН.

Официальные оппоненты : доктор Физино-математичееких наук,

профессор В.И.Быков доктор Физико-математических наук, профессор В.П.Ильин доктор технических наук, профессор Кирилов В. О.

Ведущая организация 8 Факультет вычислительной математики V кибернетики Московского государственного университета, г. Москва.

Зашита состоится £3 марта 1ЭЭЗ г. в 15 часов на заседании специализированного совета Д 002.10.ОН по защите диссертаций на соискание ученой степени доктора наук при Вычислительном центре СО РЯН/6300Э0, г.Новосибирск , пр-кт академика Лаврентьева, 6.

С диссертацией можно ознакомиться в читальном зале библиотеки Вычислительного центра СО РЙН/6300Э0, г.Новосибирск, пр-кт академика Лаврентьева, 6.

Автореферат разослан "_Ч_ 1993 г.

Ученый секретарь

специализированного ученого совета

Г.И.Забиняко

СОДЕРЖАНИЕ

стр.

Введение.................................................1

Продолжение решения по параметру.............................8

Характерный припер нелинейной краевой задачи................13

Параметризация. ...........................................1В

□6 использовании ортогональных прогонок при квазилинеаризации "параметризованной" краевой задачи..... 23 Численное исследование дискретной подели нелинейной краевой задачи. ..................................25

Дискретные модели, основанные на сплайн-интерполяции. ■■■■■■■■■■■..•■■.■■.•.■•■■«■■■»•>■■■-■-..*.... 32 Программная реализация. Комплексы программ

"SYSTEM", "BPR-Q" и "INTENS"..................................................35

Заключение........l ........................................39

Публикации............................................42

ВВЕДЕНИЕ

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

Ритуальность проблемы.

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

- £ -

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

Но во многих приложениях изучение зависимости решения от одного из параметров математической модели Формулируется как цель исследования. Примером могут служить математические модели стационарных процессов макрокатализа, часть» «следования которых являются такие вопросы как определение чувствительности найденного решения к изменению параметров модгли, расчет критических значений параметров, и т.д.

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

Продолжение решения по параметру и ветвление решений типа "поворот". Параметризация.

Предполагается, что пользователь программы имеет представление о решении краевой задачи при некотором "стартовом" значении параметра, не совпадающим с заданным. И тогда предпринимается попытка продолжить найденное решение краевой задачи при "стартовом" значении параметра до значения параметра, при котором решение краевой задачи представляет интерес. Здесь метод Ньютона-Канторо-вича-РаФсона (метод квазилинеариаации) выступает как Эффективное средство численного построения серии решений краевым задач, отвечающей изменению ваг за шагом значений параметра в ходе продолжения решения. Начиная с первого вага по параметру, имеется возможность эффективного задания начального приближения га счет выбора величины шага, а быстро сходящийся итерационный процесс завершается определением производных решения по параметру, которые затем используются, при задании очередного начального приближения решения, и так далее. В такой трактовке метод продолжения решения по параметру по существу является способом задания начального прибпи-жения решения краевой задачи при интересующем значении параметра и в определенном смысле пополняет информацию, которую при других, бопее благоприятных обстоятельствах, мы могли бы получить из качественного анализа пробпемы.

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

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

Вариант параметризации, обсуждаемый в докладе, имеет простой Физический смысл. Рассмотрим систему, состояние которой изучается в зависимости от параметра. Подвергнув возмущению параметр <ипи компоненту системы), будем наблюдать отклик системы на возмущение. И тогда в качестве "текучего" параметра системы для продолжения на один ваг решения, описывающего состояние системы, выберем компоненту (или параметр системы), наиболее сильно реагирующую на возмущение.

Геометрическая интерпретация.

Отметим; что графиком решения краевой задачи, в которой возможны только ветвления типа "поворот", является гладкая поверхность в пространстве компонент вектора неизвестных и скалярного аргумента, целиком состоящая из пространственных кривых, задаваемых решением краевой задачи при непрерывном изменении параметра. Это не противоречит ситуации, когда одному и тому же значению параметра будет отвечать несколько пространственных кривых, принадлежащей одной гладкой поверхности, или, что краевая задача описывает несколько таких гладких поверхностей. Аналогичная ситуация имеет место при численном анализе систем трансцендентных уравнений с параметром, когда графиком решения системы является гладкая пространственная кривая в пространстве компонент вектора неизвестных и параметра. Тем не менее возможно, что некоторые точки пространственной кривой характеризуются равенством нулю определителя матрицы Якоби.

Предложенная программная реализация продолжения реаения в сочетании с параметризацией на основе метода Ньютона—Канторовича-

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

К истории вопроса. Математические модели пленочной электромеханики.

Идея параметризации возникла в связи с математическим моделированием в пленочной электромеханике и, начиная с 1Э6Э года использовалась при решении нелинейных краевых задач в различных приложениях. (Список основных публикаций по теме доклада приводится.) Отметим, что в ряде случаев зависимость решения от параметра удавалось численно построить методом простой итерации, сходимостц которой имела место только за счет параметризации. Примером может служить расчет параметров гистерезиса пленочного электростатического реле в зависимости от жесткости гибкого электрода. В дальнейшем ( в 1Э85 году ) был предложен вариант параметризации в сочетании с методом продолжения решения по параметру дискретной модели краевой задачи, представленной системой трансцендентных уравнений с параметром. К этому же времени относится разработка первых версий комплекса программ "BPR—Q" ( численное исисследование дискретной' модели нелинейной краевой задачи) и комплекса программ"SYSTEM" (численное исследование систем трансцендентных уравнений как самостоятельная проблема ) на алгоритмических языках BASIC и FORTRAN. Тексты "BPR-Q" и "БУБТЕМ" < BflSIC-лрогранмы ) приведены в работах [1S, 13,17,24,£63. Краткие характеристики версий содержатся в докладе.

Дискретные модели краевой задачи.

Предлагаемые диснретные модели краевой задачи определяется заданием сетки по скалярному аргументу, а затем приближенным представлением интегральных равенств, вытекающих из системы дифференциальных уравнений, на каждом из отрезков сеточного разбиения с максимальным нагом Нгпах . При этом используются приближенные выражения интегралов по Формулам Котеса алгебраической точности 1, 3 и 5. что соответствует аппроксимации дискретной модель» сеточных

9 L

значений искомой вектор-Фуннции с погрешность» порядна Нмах ,Нтах^ и Нтах^ на решениях.(Насколько известно автору, такой подход к построении дискретной нодели краевой задачи применял H.Keller с привлечением линейной интерполяции решения на каждом из отрезков). Полученная этим способом система трансцендентных уравнений характери-

зуется блочно—двухдиагонапьной структурой матрицы Якоби. Важно отметить, что параметризация в предлагаемой варианте позволяет сохранить разреженность матриц систем линейных алгебраических уравнений на итерациях. Поэтому применение метода Гаусса, учитывающего разреженность матриц, обеспечивает экономичность процесса продолжения решения. Зту же цепь преследует процедура адаптации шага по"текуще— му"параметру. Наконец, при Формировании дискретной модели необходимо предусмотреть адаптацию сетки ( как расположение узлов,так и их число ) в процессе продолжения решения с тем, чтобы предусмотреть возможность появления решений с большими градиентами. В "BPR-Q" при адаптации узлов сетки используется правило, заимствованное из задачи интерполяции гладких функций сплайнами. При этом узлы сетки требуется выбирать таким образом, чтобы распределение погрешности интерполяции было равномерным по отрезкам разбиения.

Комплексы программ. Новизна.

Из публикаций известно о существовании нескольких комплексов <пакетов) программ < авторы Е. Doedel, M.Kubicek как примеры) в связи с рассматриваемой проблемой, в которых используются иные способы построения дискретных моделей < применение коллокационных методов) и продолжения решения. Как и предлагаемые в докладе комплексы они являются рабочим инструментом исследования. Однако сопоставление различных вариантов комплексов по эффективности затруднительно, в частности, из-за не полной доступности к текстам программ, входящих в комплекс. В связи с этим можно указать ряд математических моделей, исследованных при помощи "BPR—Q,r или "SYSTEM", которые могут служить тестами для сопоставления различных подходов. Это краевые задачи с множественностью решений, с образованием внутренних пограничных слоев, перемещающихся в процессе продолжения решения по параметру, с проблемами задания начального приближения, выбора стартовой позиции, и т.д. Подобного рода особенности моделей могут быть препятствием для применения компленсов программ, в которых не предусмотрена параметризация или адаптация сетки. Как показывает практика, численное исследование нелинейных проблем наиболее эффективно в режиме диалога, который позволяет организовывать вычисления с учетом информации об объекте исследования. В этих условиях программные средства, доступные пользователю в виде "черного ящика", настроенные на задачу определенного типа, очевидно, имеют ограниченное применение. Разработка "BPR-Q" и "SYSTEM" явилась реакцией на поток задач, с тем чтобы иметь возможность оперативно адаптировать комплесы программ к различным приложениям. Обратная связь проявляется в том, что практика общения с конкретными математическими моделями подсказывает различные модификации компленсов программ. Например, прокладывание маршрутов по двум параметрам модели (каждый из этих параметров последовательно играет ропь сналяр-ного параметра задачи) для выхода на начальное приближение решения

- в -

при.заданной совокупности параметров модели, в частности, при изучении изолированным решений.

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

С разной степень» краткости мы остановились во введении на содержании модулей цепочки, входящей в понятие численный эксперимент: объект—математическая модель—дискретная недель—метод—комплекс программ—расчет. Конкретные приложения комплексов программ приводятся в докладе.

Личный внлая автора.

К личному вкладу автора относятся 1) разработка дискретных моделей; £) разработка метода продолжения решения в связи с предложенным вариантом параметризации; 3) разработка на алгоритмическом языке BASIC версии "BPR-Q" и "SYSTEM", а также комплекса программ "INTENS", описание которых даны в докладе; 4> научное руководство программной реализацией других версий этих комплексов и их приложений; 5> численное исследование конкретных математических моделей, отраженное в списке публикаций. При этом в совместных публикациях автору принадлежат Формулировки дискретных моделей, в частности, основанных на сплайн-интерполяции, и результаты численных исследований с применением BASIC-программ.

Схема доклада.

В первом параграфе дается описание метода продолжения решения краевой задачи по параметру. Однако' смысловая нагрузка параграфа не ограничивается изложением известных Фактов. Для дальнейшего изложения было важно обратить внимание на то обсоятельство, что для достаточно хорошего начального приближения в методе квазилинеаризации линейные краевые задачи на итерациях столь же хорошо обусловлены, как и линейная краевая задача, определяющая производные решения по параметру. Это следует из близости по норме соответствующих Функций Грина. Кроме того, здесь рассматривается процедура адаптации шага по параметру задачи, применяемая и дпя продолжения решения по "текущим" параметрам, которая используется в комплексах программ "BPR-Q","SYSTEM" и "INTENS".

Во втором параграфе приводится пример нелинейной краевой за—

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

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

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

Пятый параграф содержит описание дискретных моделей нелинейной краевой задачи, определяемых на сетке с использованием Формул Котеса, и процедуру адаптации сетки. Здесь же дан вариант метода продолжения решения в сочетании с параметризацией для систем трансцендентных уравнений.

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

В седьмом параграфе приводятся краткие описания комплексов программ "PPR-Q","SYSTEM" и "INTENS". Разработка "INTENS" связана с исследованием параметрической чувствитепьностьи решений задачи Коши для автономной системы уравнений.

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

Завершая введение, заметим, что сколько-нибудь полное перечисление публикаций по данной теме является трудной задачей. Не претендующий на обзор список литературы, имеющей непосредственное отношение к теме доклада, дан в работах С24,28,2S], по которым составлен доклад. Вместе с тем следует подчеркнуть, что основу исследований по данной теме составляют труды таких ученых, как Ррнольд В.И., Бахвалов Н.С., Беллиан Р., Годунов С.К., Канторович Л.В., Келлер Г., Оден Дж., Рейнболдт В., Самарский A.A. и др.

ПРОДОЛЖЕНИЕ РЕШЕНИЯ ПО ПЯРЙМЕТРУ

Пусть требуется найти компоненты у^^х), уг<х),..., уп(х) вектор-функции у(и) скалярного аргумента х, принадлежащего конечному интервалу <а,Ь>, из решения краевой задачи для системы обыкновенных дифференциальных уравнений, содержащей скалярный параметр (2,

х€<а, Ь>

ауг/с1х=^ <х,у^ , уг,...,уп, 0)

с}уп/£1к=Гп<к,у<, у?, . .., у„, В)

с краевыми условиями

д^ <у4 (а), . . ., уп<а), у {<Ь), ... , у„ <Ь>, Е)=0 а1<у1<а),...,уп(а),у1<Ь>,...,уп<Ь>,0)=0

Э„<У4 <а>. ■ • • ,у„<а>, у{ <Ь).....уп(Ь),0>=0

Здесь и Э^ , 1 = 1, Н, ... , Г1 - достаточно гладкие функции по совокупности аргументов в некоторой области их определения. Функции т- * - ?представляют п независимых краевых условий. Полагая и 3 ^ компонентами вектор—функций f и д соответственно, ¡запишем краевую задачу в векторном виде :

хе (а,Ь), dy/dx=f (х, у, 0) , (1)

д (у <а) ,у<Ь),0)=0 . <£)

В частности, краевые условия <Е> могут задаваться независимо при х=а и х=Ь :

1 (у (а), 0) =0, г (у <Ь), 0) =0, <3>

где 1 и 1- — вектор-Функции размерности г\0 и Г|-пв, О ( п 0 < п , имевшие компоненты Ц , 1 . . . , 1 и ' ' • ■ • ' —По соответственно. Если правая часть (1) является со - периодической функцией по х, то для отыскания СО - периодического решения краевой задачи ставятся условия :

у (О) =у < 10 ) , Ш > 0, (А)

т.е. а=0, Ь= и> . Возможны и другие варианты задания краевых условий. В дальнейшем, для удобства записи <2) и связанных с

- э -

этим выражений, будут использоваться обозначения :

d =у(а), jj =у (Ь). (5)

Обычно, в определении правых частей <1) и краевых условий (S) участвует множёство параметров. Среди них выделяется параметр Q, если представляет интерес поведение решения y<x,Q> краевой задачи <1), <2) в зависимости от значений параметра О. Эффективным средством численного построения этой зависимости служит метод Ньютона-Канторовича-РаФсона (метод квазилинеаризации) в сочетании с процедурой продолжения решения по параметру.

Пусть требуется чиспенно построить зависимость y(x,Q) на от— — Col резке CQ, Q3 по Q и пусть Y (х) - вектор-ч>ункция, которая может

служить начальным приближением решения краевой задачи (1), (2> при Q = Q, Тогда, согласно методу квазилинеаризации, строится последовательность вектор-Функций <х>, к=1,2,..., сходящаяся к y(x,Q) при к—>°о, которые опредепяются из соответствующей последовательности линейных краевых задач :

х£(а,Ь>, ц(х> =Ytlí+^<x) , к=0, 1,£, ...

du/dx = йГК1(х)*м + FtKl<x), (6)

SM*u<a> + Т1Й*и(Ь> + Ч tUl = О,

Г(Л Г и"]

й (х) = fy (х, U(x), Q), U <х> = Y <х),

F^4x> = f(x,U(x),Q) - fl^(x)*U<x),

8М- g .<LHa), U(b>, D>, T^ = g <U(a), U<b>, Q), rt f

g (U (a) , U (b), Q> - (a) - T^*U(b),

при Q = D, fy, g^, g^ - соответствующие матрицы Якоби, обозначения cL и ji даны в (5). Предполагается, что краевая задача (1), • <£) имеет решение при Q = Q, дифференцируемое по параметру Q , a Y (х) - вектор—Функция, принадлежащая области притяжения y(x,Q) в методе квазилинеаризации, так что при любом к > О краевая задача (6) однозначно разрешима. Иными словами, однородная краевая задача

х £ (а, Ь), к=0, 1, 2, . .. dz/dx = R^(x)*z, S^wzla) + (Ь) = О

где

имеет только тривиальное решение. Следовательно, существует ограниченная матричная функция Грина однородной краевой задачи (в), в том числе и краевой задачи

к £ (а,b), dz/dx = й(х)*г, S*z<a> + T*z(b) = О, <7>

с матрицами ß(x),S и Т , определенными на решении y<x,Q) : Q = D,

fi(x) « fy (x,U<x),Q), U (x) = у (x, Q) ,

S = д^ШСа^ШЫ.О), T = (a), U (b>, Q) .

Для дальнейшего важно отметить, что отсюда следует однозначная разрешимость линейной краевой задачи относительно производных yg<x,Q> по Q решения (1), <2> при Q = Q :

х S (а, Ь), v(x) = y0<x,Q> , U (х) = y<x,Q> ,

* (8) dv/dx = fl<x)»v + Р(х> , S*v(а) + T»v<b> + Ц/ = о

Здесь Р<х) и - векторы производных по Q :

Р(х) = f6<x,U<x>,G) , Ц/= g^<U<a>,U<b>,<3>

Как известно, условие однозначной разрешимости краевой задачи (в) имеет вид :

dettS*4>ta) + Т*Ф(Ъ)Э £ О,

где 4><х) — фундаментальная матрица решений уравнения

х С (a,b), dz/dx 0<x>*z.

При этом существуют матричные функции Грина G0(x) и G(x,t) :

х, t £ (а,Ь>, В0(х> = - Ф(х)*С5*Ф(а) + ^«(b):"',

I 0^(x,t>, если а <= х <= t G<x,t> = I

I Gj(x,t), если t ( х (= b

Ga (x, t) = - Ф(х)*С5*Ф(а) +■ Т»Ф(Ь)3 *Т*Ф<Ь)*Ф (t>

Ga(x,t) = <Hx>*-CE - С 3«ФСа> + Т*Ф(Ь>] *Т*Ф<Ь)}*Ф <t>,

где Е - единичная матрица, и решение краевой задачи (S) может

быть представлено в виде :

Ь

v(х) = jo(x, t)P(t)dt + Ge<x)y (9)

а

Из ограниченности норм G(x,t) и 0в(х), а также норм правых частей Р(х) и Ц7 , вытекает ограниченность нормы вектор — функции y^(x,Q>, обозначенной в (в) через v(x>. <

Знание производных позволяет эффективно задать начальное приближение решения (1), (2) при Q = Q + hQ, еспи шаг по Q достаточно мал, например, в виде :

Y^(x) = у (х, Q) +■ y^(x,G)*hQ

(10)

Ск1

Затем вновь строится последовательность Y (х), к=1.£,..., согласно методу квазилинеаризации сходямаяся к y(x,Q + hQ), находится вектор-функция y^(x,Q +• hQ) из решения линейной краевой задачи, и так далее. Тем самым численно строится зависимость решения (1),(2) от параметра Q на заданном отрезке С □,□ 3 по Q.

Практически, при реапизации метода продолжения решения по параметру добиваются того, чтобы максимальное значение нормы вектор-функции невязок Z<x> оказалось меньше заданного числа EPS s

tKH] М

max i IZ(x) I I < EPS, Z(x> = Yk (x) - Y J(x) , (11)

3C

rwi

и тогда Y ^(x) принимается за решение краевой задачи (1),(2). Кроме того, вводится ограничение на число итераций, за которое условие (11) должно быть выполнено. Формально, в силу диФФеренци-руемости решения по параметру Q , итерационный процесс сходится в смысле (11) за число итераций к ( к* , к^ > О, если шаг no Q достаточно мал. Если к = , а максимум нормы вектора невязки по-прежнему больше EPS, то попытка построить решение краевой задачи при Q = Q + hQ с начальным приближением (Ю) объявляется несостоятельной. Очередная попытка продолжить решение состоит в том, что уменьшается шаг no Q, например, в два раза, и ичется решение краевой задачи с начальным приближением вида (10), где шаг no Q равен hQ/2, и так дапее, до тех пор, пока начальное приближение не окажется в области притяжения решения.

Отметим, что итерационный процесс должен сопровождаться контролем за тем, чтобы очередное приближение не выходило за область определения вектор-Функций f(x,y,Q), g(0i,jl,Q> ( ипи за выполнением каких—либо иных заданных ограничений), опять—гаки за счет выбора достаточно малого шага по параметру Q. Естественно, что при этом требуется указать минимально допустимый шаг : если при очередном делении "пробного" шага по параметру Q пополам "пробный" шаг стал меньше допустимого, то продолжение решения краевой задачи (1),(2) по параметру □ прекрашается.

Вполне возможно, что наблмдаемов уменьшение шагов по пара-

метру □ связано с близостью О к собственному числу краевой задачи (В) О = а^. . При этом краевые задачи (6) и (В) становятся плохо обусловленными, з частности, с приближением (3 к собственному числу неограниченно растет норма решения краевой задачи (В). Поэтому задание очередного приемлемого начального приближения У^(х> сопровождается все большим измельчением шага по а . В результате имеет место медленное приближение к точне ветвления П = по па— рамеру а решения (!>,(£>. Тем не менее, иногда удается избежать затруднения вычислительного характера в такого рода ситуациях за счет регулярного обращения к процедуре параметризации, которая будет рассмотрена в дальнейшем.

Если условие (11) оказалось выполненным при к < кЛ , то очередной "пробный" шаг по О увеличивается по сравнению с предыдущим реализованным шагом продолжения решения в заданное число ЕРБО раз, ЕРЗО >1. И так до тех пор, пока шаг не станет больше заданного максимально допустимого шага по О, который в этом случае берется в качестве "пробного" шага.

Возможно, что в стартовой позиции (0=0) число итераций превысит заданное ограничение, либо очередное приближение не будет удовлетворять условиям, связанным с областями определения f и д. Тогда возникает проблема задания нового начального приближения краевой задачи, которая носит самостоятельный характер.

Завершим описание технологии варианта продолжения решения замечанием о том, что начальное приближение в виде (10) берется только на первом шаге. Для последующих шагов Эффективным является гадание с привлечением предыдущего решения и его производных по параметру 0. Пусть найдены решения у(х,(3 ) и ) краевых задач (1),(£) и (8> при 0 - О^, 1 = 1,2,..., 0^=0, и пусть у(х,(3^ ) и следующие решения указанных краевых задач при продолжении решения по параметру 0 . Задав первый "пробный" шаг Ь(3 =ЕР30*<(3^4| - в качестве начального приближения У^(х) решения краевой задачи (1),<£> при О = + Ы3 возьмем следующую вектор-функцию :

t = 1 + МЗ/<0. 0.), к > 1, 1+1 1

уСо^(и> = <ь-1)2 *(1+2*и*у<х,с!£> + *(з-г^)*у(х,а^(>+ аг>

+ (а. о,)*ь*("Ь—1>*с<ъ—1>*уА(х,о*) + <х,а;„ >э г 7 & ' г ' ьм

с погрешностью порядка <аг-Н ~ ^ * " ® дальнейшем значение 0 адаптируется к условию (11), где к ( к^. .

Процессу продолжения решения по параметру 0 можно дать следующую геометрическую интерпретацию. В <Г1+1) - мерном евклидовом пространстве (х,у, , ..., уп ) графиком решения краевой задачи (1), (£) при заданном значении параметра 0 является гладкая простран-

ственная кривая :

Уг = уг<к,0>' 1 = 1'2....... 0 £ с 3

Непрерывно меняя О, мы получим гладкую поверхность Бц , целиком состоящую из пространственных кривых 1-4,. При этом каждому значению параметра О соответствует единственная пространственная кривая, если собственные числа краевой задачи (8) не принадлежат отрезку [ Й,В ] . Очевидно, площадь ЭШ) полосы поверхности Бя между двумя пространственными кривыми, задаваемыми решениями у<х,й) и у(х,0>, СЗе С 5,5 также может рассматриваться в качестве параметра краевой задачи (1),(2>. С использованием евкпидовой нормы вектора эпемент площади сШ определится по Формуле !

г г г , г </а

с!Э = ¿(3* 5с(1+1^П )*ИуйП - <ЛУа> 3 ах , (13)

а

Здесь f = f <х, у (х, О), О), у^ = у^(х, у <х, 0), П>, (Г, у^> - скалярное произведение. Отсюда, с учетом О), следует, что Б <13) - непрерывно дифференцируемая Функция, причем

Ь

аз/ао > 51'Уф''^х > о .

а

Более подробно вопрос о диФФеренцируемости решения краевой задачи (1),<£) будет обсуждаться в связи с параметризацией.

ХАРАКТЕРНЫЙ ПРИМЕР НЕЛИНЕЙНОЙ КРАЕВОЙ ЗАДАЧИ

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

х€<°»1>1 Уг<0) =0, Уц (1> =0

(14)

2

ау, /ах = у2 , ау^/ах = - 0/<1-у( )

описывающую положение гибкого натянутого электрода ппеночного электростатического репе в зависимости от разности потенциапов, приложенных к неподвижному и гибкому электродам. Здесь 0 - параметр, пропорциональный квадрату разности потенциалов, □ >= О . Значения прогиба отнесены к ширине зазора репе, а расстояние от центра симметрии прогиба отнесено к полудлине электрода. В краевых условиях учтена симметрия прогиба. Предполагается, что зазор много меньше длины электрода. Из Формулы, задающей точное решение (14), следует, что решение существует не при всех значениях

Q > О, а лишь для GK.35. При этом каждому значению Q, О (GH.35, отвечают два решения, описывающих устойчивое и неустойчивое ( с большим прогибом ) равновесия гибкого электрода. Обозначим через ¡1/1 значение у^ (О) Функции у^(х>, монотонно убывающей от у^ (О) до

0 с ростом х от О до 1. Зависимость Q от |Ч , характеризующая указанную особенность рассматриваемой краевой задачи, имеет вид :

Q([*) = .5*<l-frf)*E\/(vï + „-^l+M 3 <15>

V -I

График Функции GKjVj ) приведен на рис.1 . Максимальное значение Q( pU, равное 0.35 , достигается при JV» = = . 30B3 и определяет точку ветвления краевой задачи (14), т.е. И*. = .35.

Физическая интерпретация возрастающей ветви графика функции Q<f<1 ), О < (И < , состоит в следующем. С ростом разности потенциалов ( параметр Q > увеличивается максимальный прогиб гибкого электрода ( параметр jv) > в состоянии, когда реле не замыкает контакты. Максимальное значение параметра Q определяет так называемое "напряжение срабатывания", а значение = — максимально возможную высоту контактного столбика, соответствующую значению

1 — ffl, помещаемого на неподвижном электроде в центре. Нисходящая ветвь графика Функции Q(f>1>, fi < 1, определяет в зависимости от высоты контактного столбика так называемое "напряжение отпускания", при котором гибкий электрод скачкообразно занимает устойчивое положение равновесия, задаваемое возрастающей ветвью графика Функции Q(f4>. Стрелками на рис.1 отмечена петля гисте-реза, характерная для работы пленочного электростатического репе.

Проследим на этом примере за различными ситуациями, связанными с численным построением решения краевой задачи (14) методом продолжения решения по параметру. Заметим, что при □ = Q = 0 краевая задача имеет точное решение у(х,0) = О. При этом функции

v^(x> = .5*<1-х4), v2(x) = -х

являются компонентами вектор-Функции v(x> = yg(x,Q> при Q = О, то есть для достаточно малых Q вектор-Функция начального приближения решения (14) принимает вид :

Y^(x) = С .5*Q»<l-x2), -Q*x ]Т ,

где Т — знак транспонирования. Следовательно, имеются все предпосылки для численного построения решения краевой задачи (14) в зависимости от Q методом продолжения решения по параметру В . Однако, как видно иа рис.1, продолжение решения по параметру Q возможно лишь до некоторого значения □ = Q^ < Q^, поскольку при Q — d(Vl/dQ = vj <0>, равно как и норма v(x, Q), неограниченно растут. С другой стороны, график Функции 0((Ч) свидетельствует о том, что

* * *

РИС.1. Математическая модель работы пленочного электростатического реле.

Графики первой компоненты у fix) двух решений краевой задачи (14) при О < Q < Здесь f<l4 - максимальное значе-

ние у^(х> первого решения, (Иг. _ максимальное значение у^(х> второго решения, О < [И* < , достигаемые при х=0 . Функ-

ция описывает равновесную Форму гибкого эпектрода под

воздействием электростатического притяжения.

Гра*ик Функции П(("<> (15), (И = у^(0>. Стрелками отмечена петля гистерезиса, характерная дпя работы Физического ппеночного электростатического реле с относительной высотой контакта 1- (И^. При а = а имеет место ветвление решений краевой задачи (14) I краевая задача не имеет решений, если □ ) , и два решения, если О < □ ( О,.

• * *

при выборе р! в качестве параметра математической модели пленочного электростатического реле, метод продолжения решения по параметру поавопяет нам численно построить О < {И> на интервапе (0,1) в силу однозначности функции □ ( ) .

Для того, чтобы воспользоваться (VI в качестве параметра для продолжения решения представим краевую задачу (14) в виде интегрального уравнения относительно Функции и(х> = у^(х):

Ь

и (к) = О* £в<х,1;>РСи(*:> (15)

а

2

где Р(и) = 1/(1-и) , Б^^) - Функция Грина линейной краевой за-

р р

х е <0, 1) , А т. /бх = О , с1г/с1х(0> = О , 2(1) = О, то есть

I 1 -Ъ, если СК=х<=^ 6(х,1;> = I

I 1-х, если Ь (х<= 1 .

Поскольку, нак следует из (15) при х = О ,

Ь

(М = В* ^б(0,Ъ>РСи<Ъ) ЗйЪ, а

то параметр (3 легко исключается из (15). В результате искомое представление краевой задачи (14) с в качестве параметра име-

ет вид : О < (V» < 1 ,

Ь Ь

и(х> = ¡Ч*уВ(х, )РСи(1: ПсИ;/ £Б(0,^РСи(1:> ЭсН;, < 16)

а а

При этом зависимость О от ^ находится по Формуле :

Ь

0 = / ^ВСС^ЮРСиМгПсН:', (17)

а

где и(Ъ) - решение (16).

Тольно что предложенное преобразование имеет место и в более о5«вм случав, В частности, оно было испопьаовано для построения зависимости С1((<\> математической модели пленочного эпвктро-стлтичаского реле с учетом жесткости гибкого электрода, сформулированной как нелинейная краевая задача ■ следующем виде : п = 4,

х е (о, 1) , ау/йх - д»у + а*я<у) ,

(18)

Уг.(0) - ум(о> - о , у4(1> - уг(1> - о ,

где й - постоянная матрица, Я(у) - вектор-функция векторного ар-

гумента у:

1 0 1 0 0 1 1 0

й = 1 0 0 1 0 1 R(y)= 1 0

1 0 0 0 1 1 1 0

1 0 l/tf2 0 0 1 1 F (у^ >/Зг

Параметр пропорционален жесткости гибкого электрода, О < # ( <00. Обозначив у/(х> через и(х), вновь преобразуем краевую задачу к интегральному уравнению (15), где теперь Б(х,1;> - функция Грина линейной краевой задачи 5

х 6 (О, 1)

х = 0: du/dx = d3u/dx3

,2 U и 0 V *d u/dx4 = d u/dx'

г

<ia' >

х = 1: u = du/dx = 0

имеющая вид

G (x,t > =

1 - t -<i/5h<l/J))*Cch(l/J) - ch(t/Jf) + + <ch < < 1—t) / g) -1) *ch (x/tf)3 , 0 <= x<=t , 1-х - (i/sh(l/j())*Cch<l/a'> - ch (x/Jf) + + (ch ( (1-х) /&>-l)*ch (t/g) 3 , t < x<=l .

Отметим, что характер зависимости СЭ(^>, = у^ (О), краевой задачи (18) остается тем же, что и на рис.1, при различных значениях ТС, О ( % ( ой

Численный способ решения уравнения типа (16) будет обсужден в дальнейшем. При этом не предполагается использование явных выражений функции Грина. Поэтому Формально он может быть использован при численном исследовании нелинейной краевой задачи вида (18), когда для линейной части оператора (возможно с переменными коэффициентами) существует матричная функция Грина. В ряде случаев, например, для уравнения вида

Pr,(x,d/dx)u + Q*F (х, u) = О,

где Pp(x,d/dx) - линейный дифференциальный оператор п—го порядка с непрерывными коэффициентами, F(x,u) — достаточно гладкая функция своих аргументов, такой подход может оказаться Эффективным и в отношении вычислительных затрат.

ПАРАМЕТРИЗАЦИЯ

Рассмотренный пример интересен тем, что в нем удалось обойти ветвление решения по параметру О (ветвление типа "поворот") за счет использования компоненты у^(к) при х = О в качестве параметра {И для продолжения решения вместо параметра О. Регулярную процедуру выбора параметра для продолжения решения на один шаг и связанную с этим переформулировку краевой задачи, мы будем называть параметризацией. Сочетание метода продолжения решения краевой задачи по параметру О с параметризацией расширяет возможность численного исследования зависимости решения краевой задачи от параметра О. Параметризация устраняет вычислительную проблему ветвления решений типа "поворот" по параметру О, возможность которого следует предусмотреть при организации численного эксперимента, если не доказано обратное. В данном примере, благодаря параметризации, численно строится гладкая поверхность Зц , являющаяся графиком решения (14) в пространстве (х^^у^) с элементом плвщади , выраженном через производные решения по параметру : О < (^<1, Ь (

ЙБ = $С<1+1 1 |г)*1 1у)„11г - ^f,y(/¡)гl'Zdx.

а

(Подразумевается евклидова норма вектора.> Выбранный при параметризации параметр для продолжения решения на один шаг ( им может оказаться и параметр (3 ) мы будем назыввать "текущим" параметром.

Б дальнейшем будет предполагаться, что решение краевой задачи (1),(2> вентор-Функция у(к,0> определяет в (п+1) -мерном евклидовом пространстве (к,у) гладкую поверхность Бя . Сочетание метода продолжения решения с параметризацией позволяет численно построить Бц в той области изменения параметра О, в которой указанное предположение выполняется. При этом допускается, что Эр принадлежат несколько пространственных кривых задаваемых ре-

шениями (1), (2) при одном и том же значении параметра 0, то есть имеют место ветвления решений типа "поворот" . Вполне возможно, что краевая задача (1),(2) определяет несколько гладких поверхностей Бя . И тогда речь идет о численном построении одной из них, содержащей пространственную кривую у(х,СЗ> стартовой позиции процесса продолжения решения.

Априорное определение параметра , как это было в рассмотренном примере, далеко не всегда возможно. Возвращаясь к краевой задаче (1), (2), будем рассматривать в качестве возможных "текущих" параметров для продолжения решения как параметр 0, так и компоненты вектор-функции у(х) в некоторой точке х « с интервапа (о,Ь). Получив решение краевой задачи (1),(£>, если О — "текущий" параметр, или решение "параметризованной" краевой задачи, в которой О подлежит определению по заданному значению "текущего" параметра, следует вновь опредепить"текущий" параметр краевой задачи

(1), (£), по которому решение будет продолжено на очередной шаг.

Остановимся на процедуре выбора "текучего" параметра. Простые рассуждения показывают, что если решение краевой задачи (1),

Л/ -V _ —

(2) дифференцируемо по О при СЗ = СЗ, Об СО, Н] , то это решение

Л/

дифференцируемо и по параметру ¡VI = ук(с), с е(а,Ы, при (VI = = = ук<с, О), если к-ая компонента у!<(х) вектор-функции v(н) = = уц(и,0), к = с , отлична от нуля. Иными словами, у(х,СЗ) можно интерпретировать как решение "параметризованной" краевой задачи:

х е(а,Ь), йу/а* = I1 <х,у,а<у-1>>,

(13)

д(о/ , £ ,□( =0, Ук<с> - р , с&(а,Ь)

с параметром (V* при р-1 = £5 , определяющей вектор-функцию у(х, |И) и О <!*)). При зтои у(х, рл) и дифференцируемы по (VI , т.е.

линейная краевая задача

х £ (а, Ь), и (к) = У^(х, (Ч ), 9 = Ор ( р1>,

dw/dx = Я(х>*м + 0*Р(х>, 3*м(а> +• Т*к(Ь> +■ 8»^ = О, (20)

|л1к(с> = 1, с ® (а, Ь)

г-

где ул = ул , однозначно разрешима. Здесь и О р — производные

по 5<1 . Остальные обозначения те же, что и в (3).

Действительно, решением (20) является вектор—Функция

и(х) = у(х)/ук(с> = 0(И( рЬ*у(х), >^<0 + О, (20'

и это решение единственно. Если предположить, что существует два решения (£0), которые отметим инденсами 1 и 2, то ик разность является решением краевой задачи :

хе(а,Ь), Ц (к) - иС41(х>-мСг:1(х), Ч = 0С,]( р р )

dW/dx = Я (X) + q»P<x), 5*И(а) + Т*И(Ь)-+ ц* <С = О, (£1)

Ик(с> =0, се (а, Ь),

где М|((с) - к-ая компонента Ы(х) при х = с . Пусть д = О. Тогда Ы(х) <= О, в силу того, что краевая задача (7) имеет только тривиальное решение и, следовательно, вектор-Функция ии(х) (20') - единственное решение (£0). Пусть д ^ О. Полагая в (21) И(х> = » Я*С0(х), сраау убеждаемся, что Ь)(х) = у(х) в силу однозначной разрещимости нраевой задачи <8) . Но (Он (с) = 0 согласно краевому условию (21), что противоречит предположению ф О. Таким

)

образом, в (20') определено единственное решение (20).

Очевидно, верно и обратное утверждение. Если решение "параметризованной" краевой задачи (19) дифференцируемо по параметру [V) при ул = и ^ О, то это решение дифференцируемо и по параметру 0, т.е. краевая задача (в) однозначно разрешима при 0 = 0 = 0(^1), а ее решение пропорционально кч(х) :

у(х> = м(х)/Шр, ( >, С3р( (Ч ) ф О.

Отсюда следует, что решение краевой задачи дифференцируемо и по параметру Д = у^(с|>, 4 6 (а, Ь), если 4 О >

Таким образом, решение у(х,0) краевой задачи (1), (£>, полученное методом продолжения решения по параметру С! (и, следовательно, дифференцируемое по 0 ), может быть продолжено по параметру (VI = у^(с> как решение у(х, Г"), 0<р«1) "параметризованной" краевой задачи (19), если Уц<с> ^ О . Выберем для этого достаточно малый шаг Нр по ¡и. Метод квазилинеаризации, в котором начальное приближение решения у(х,(« ), 0((Ч ), ул = ул + Нр , имеет, например, вид :

УГ°\х> = у(и,р > + Ур(х, Р> >*Нр< , = 0)

позволяет найти решение у<х,ул >, 0({*>, дифференцируемое по параметру . И так далее, если всякий раз из решения "параметризованной" краевой задачи будет следовать, что О ^ ( (4 ) ^ О.

йдаптация шага по "текущему " параметру организуется точно так же, как и по параметру О. Предварительно, после того,как най-дсжо решение у(х,0г-^|> краевой задачи .(1), (£> методом продолжения решения "параметризованной" краевой задачи (19) по "текущему" параметру Д и определен новый "текущий" параметр (И , производные решений (1), (£), отвечающие £3 = и <3 = > пересчитываются

как производные по параметру

Как известно, обусловленность краевой задачи (в), или (£0), определяется оценкой нормы матричной Функции Грина. С учетом (14) решение краевой задачи (£ф имеет следующее интегральное представление :

Ь

м(х) - $В(х, + В0(Х)«И ,

а

где матричные Функции Грина и В9(м) пропорциональны Б(х^>

и Бе(х) соответственно !

В(х^) = 6(х^)/ук(с), ^„(х) - Б0(х> /V ¡¡(с), vк<c> ф О .

Поэтому переход от "текущего" параметра 0 с производными v(x) решения краевой задачи (1), (2) по параметру О при 0 = 0 к "теку-

* «■ *

РИС.2 Возможные типы ветвления решений краевой задачи (1), (2>, исследуемые при продолжении решения в сочетании с параметризацией.

/

а

(3«

Зависимость О от "текущего" параметра в окрестно-

сти собственного числа краевой задачи (8> : а) "поворот", б)"перегиб" с вертикальной касательной.

* # *

- ее -

щеку" параметру с производными м<х) решения "параметризованной" краевой задачи <1Э> по параметру (VI при ¡и - = ук(с,0) целесообразен, если 1ук(с>1>1. При этом краевая задача <20) оказывается лучше обусловленной, чем <в).

Принимая во внимание зто обстоятельство, можно предложить простую процедуру выбора "текущего" параметра краевой задачи (1), <£). Пусть из решения краевой задачи <в) известна вектор-Функция у<х) и

IУк<с) I = глах гоах IV • (х) I = гпах I IV <х> I I .

зе а *

Если ^к(с).1>1, то (vi = ук<с> становится "текущим" параметром. В противном случае решение продолжается по параметру С! . Если найдено решение краевой задачи (20) , то определяется максимальное значение нормы м(х) ,

Ы{,(с1)1 «= шх! 1н(») II , с ос

и 1и-(с1>1 сравнивается с 1(3^,(^)1. Пусть |н|<с1)1>1 и 1м^<с1>1> >10^(^)1. Тогда в качестве "текущего" параметра выбирается' = " у^<<3>, а 2(х) «■ к(х) /»^(б) ~ вектор-Функция производных решения "параметризованной" краевой задачи <19) по параметру Я при % = Д т у^<с1, {<1 ), ал< й > = СЭ^ ¡4)/и^ <с!>. В противном случае "текущий" параметром либо остается , либо становится О.

Рассмотрим случай, когда 0 = - собственное число краевой задачи (в), принадлежащее отрезку С 13, б ] , на котором численно исследуются решения краевой задачи <1>,(2). В силу предположения о гладкости поверхности Эц решение у <к, СЗ^) может быть найдено как решение у(х,ул>, И(р1) "параметризованной" краевой задачи (19) при «• то есть 0<|»1*3 • При этом О^у**) =

= О, а краевая задача

хе<а,Ь>, с!н/с|х т й(х)*м, 8»и<а) + Т»и(Ы » 0,

<22)

мк<с) = 1, се<а,Ь>

обязана быть однозначно разрешимой. Ветвление решений краевой за-задачи <1>,<2> по параметру 0 имеет тот же тип, что и в краевой задаче <14> рассмотренного примера. Как и в (14), параметризация позволяет непрерывным образом продолжить решение в окрестности 0 по параметру ^ = ук<с). Варианты графиков функции £3<(1> при атом представлены на рис. 2.

Пусть в процессе продолжения решения краевой задачи "текущий" параметр , выбираемый согласно указанной процедуре параметризации меняется от некоторого значения до Тогда пространственные кривые 1~ч заполняют полосу поверхности Эц , площадь АЗ которой выражается Формулой < с использованием евклидовой нормы

е=хтора )

Ь

" " Ч |-)*1 1ул,1 I

ь а »4/2

р I, % а а '/л

йЗ = ] (1 + 1 I Г1 I ) #1 !У(Ч,| I - (Г, у > з Ьх ,

{л' а

Площадь полосы 3 поверхности Эя, ограниченная графиками решений у(х,0) и у(х,0) , О С СО , а], будет суммой таких Д 3, отвечающей числу смен "текущих" параметров. Легко убедиться, что решение "параметризованной " краевой задачи дифференцируемо по 5 :

Ь

2. 2. 11/2 у5 = уп/ !са+'т 1 '*''1 ~ <1Г> 3 .

причем

Ь

I |у31 I ( 11Ур1И/ »у I .

а

Таким оЗразом, 2 моиет служить универсальным параметром для про-должэнин решения краэвой задачи, вели выполняются указанные ра-нея предположения о поверхности 84.

ОБ ИСПОЛЬЗОВЙНИИ ортогональных ПРОГОНОК

ПРИ КВЙЗИЛИНЕЙРИЗЙЦИИ ПАРАМЕТРИЗОВАННОЙ КРАЕВОЙ ЗЙДЙЧИ

Из предыдущего описания параметризации следует, что продолжение реяенин краевой задачи (1), <2) можно организовать, обращаясь к методу квазилинеаризации нелинейной краевой задачи <1),<Э) или соответствующей "параметризованной" краевой задачи (19). Как известно, эффективным методом численного решения линейных краевых задач, которые в данном случае возникают на итерациях, является метод ортогональных прогонок С.К.Годунова. В связи с тем, что "параметризованная" краевая задача несколько отличается от "стандартного" вида, кратко остановимся на методе квазилинеаризации и решения линейных краевых задач применительно к (19).

Пусть уСоЗ(х) - начальное приближение решения у(х,щ> и ) краевой задачи при ^ = ¡V) . Согласно методу квазилинеаризации строится последовательность вектор-Функций У^(х) и значений , .......... из решений краевых задач I

х 6 <а, Ь), и = ч » ,

йи/йх - А)*и + я*РС^(х) + Я^(х),

(£3)

3М*и<а> + ТВД*и(Ь> + О,

и^(с) = р , с 6 (а, Ь).

при ftí = . Здесь

pM(x) = f0 (х, и (х) , 9 >, U < х ) = YC|tl(x), в = Q^ RM(x> = FW(x) - В*РГ|£3(х>, g» (L! (а), LH b) , В),

□стальные обозначения введены ранее.

Очевидно, что для ufx) имеет место представление с

ц(х> = Ф(х)*u(a) + q*V(x) + Z(x), (24)

где Ф(х) - Фундаментальная матрица ревений,

TlO

йФ/dx « <V J<x>*$, Ф(а) = Е,

а

Ь

-4 г,.!

V <х) «= Ф(х)* J Ф (t)Pl J(t)dt, V (а) =0, а

I м

2<х) = Ф(х)* J Ф <t)R <t)dt, Z (а) =0. а

Здесь Е - единичная матрица. Согласно (24) вектор-Функции и<х) и q определяются по методу прогонки, в которой компоненты вектора и(а> и q играют роль параметров прогонки. При этом решения соответствующих задач Коши относительно Ф(х), Vtx) и Z(x) следует сопровождать процедурой ортогонализации, предложенной С.К.Годуновым . Формально параметры прогонки находятся из системы линейных алгебраических уравнений :

Т,"1'^*Ф<Ь)3*и(а) + TM*V(b)3*q -

» Tt¿1*z<b), (25)

n

^T%^<c>*Uj<a) + V^<c)*q = 1 - Zjíc), j»l

где Ф^(с) - компоненты матрицы Ф(с), V^(c) и Z^(c) - - е компоненты векторов V(C) и Z(C>, u^(a)- j-я компонента вектора и (а). Разрешимость системы (25) следует из разрешимости краевой задачи.

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

ЧИСЛЕННОЕ ИССЛЕДОВАНИЕ ДИСКРЕТНОЙ МОДЕЛИ НЕЛИНЕЙНОЙ КРАЕВОЙ ЗАДАЧИ

Метод продолжения решения по параметру в'сочетании с параметризацией Формулируется достаточно просто применительно к системе трансцендентным уравнений со скалярным параметром П вида

«Ч <и,, иг,.„, и > - О, Ф2<и<, иг,..., и„, а > = о,

ф„(и,, иг,..., им, а > =о,

или

Ф(и,0> = О, (26)

г,кэ ФШ,С1> - достаточно гладкая вектор—фуннция размерности N п компонентами Ф^ (и, О), Фг Ш, О , . . . , Ш, □), каждая из которым ется функцией векторного аргумента и размерности N с компонентами , 1)г, . . . , и скалярного параметра О . При этом численный эксперимент состоит в построении вектор-функции 11(0), задаваемой решением системы (26) в предположении, что при О е- Г 0, О 3 графиком 11(0) является гладкая пространственная кривая Сц в (N+1) -мерном евклидовом пространстве Ш.|, иг,..., Иц , О). В дальнейшем будут предложены некоторые способы приближенного представления нелинейной проблемы (1),(2) в виде системы (26), которую в этом случае условимся называть дискретной моделью краевой задачи. Введем в рассмотрение сетку по х :

а = щ < хг < ... ( х^ = Ь , II - = Хдц- х^ . <Э7)

При этом компонентами вектора и будут служить приближенные сеточные значения вектор-фун*ции у(к,СЗ> . Остановимся на деталях метода продолжения рэыения системы (26) в сочетании с параметризацией.

Как известно, гладкость Cq означает, что ранг матрицы Якоби [фу 3 в окрестности Ся равен N. При этом в качестве "теку-

щего" параметра системы (26) для продолжения решения может быть использован как параметр 0, если ^ О, так и 11^, если опре-

делитель матрицы [Фу ^ с вычеркнутым к—м столбцом отличен от

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

Ф(Х,^ ) = О , (^Ф^ ■ Ф О , (28)

подчеркивая, что ^<1 — выбранный "текудий параметр", а X - вектор размерности N . среди компонент которого может быть параметр 0.

Пусть X ' - начальное приближение решения Х(р > системы

\

(£8) при ¡VI = ^ . Тогда согласно методу Ньютона-Канторовича-РаФ-сона ( в дальнейшем метод Ньютона ) строится последовательность, векторов X®, = £,..., отвечающая последовательности систем линейных алгебраических уравнений

№ = х^ +

(Н9)

относительно поправок Z ^ к приближениям X ^ при р/1 = f . Эта

последовательность сходится к Х(р), если X^ принадлежит области

притяжения X ( > . Далее, находится вектор производных Хр,(р) как решение системы

ФХ<Х( f<«>, (vn*Xp= - « (Х( p¡ ), pi ), tíctO^ÍX ( ¡vi ), {vi ) ф О (30)

при ри = Jrt , и выбирается новый "текучий" параметр Я Для задания начального приближения V Col системы

«(V, ¿ ) = О (30' )

при (3 = А + НД , Д ), например, в виде :

vW = viXj + v¿( А )*НЛ

где НЛ - достаточно малый "пробный" шаг по Я ■ Индекс компоненты Хк( JV, ) определяется процедурой параметризации.

Лпя выбора "текущего" параметра системы (£Б> можно предложить тот же способ, что и при продолжении решения краевой задачи (1>,(£). Пусть Х^ , Хг,..., X- компоненты вектора .X, a dX</d)q, dX^/djVf, . .., dX^/djvi — компоненты вектора Xjy| при (И = jí . Определим максимальную по модулю компоненту вектора X (Щ :

ldX(¡/d¡m I = tnax (I dX{/d]*il,___, ldX#/dj*il> .

Если ld)(R/d{»i I > 1 и ldXK/dj^ I > IdO/dp I , то ^ = UR - новый "текущий" параметр системы (30') , где (повторяя сказанное отно-ситепьно системы (£8) ) V - вектор размерности N , составпенный, как и X , из компонент вектора U и параметра Q за исключением "текущего" параметра ^ . В соответствии с опредепением вектора V

Л'

компонентами вектора V при Л = Л являются

úX-x/úfi = (dX^'/d^ )/(dXn/dffl ) , i = l,£,...,N , i^k

Л*

и l/(dXn/dj«| ) при fW| = JV1 Напомним, что среди компонент вектора V может быть и параметр Q . Согласно Формулам Крамера из (30)

следует, что

1с)Я /с1(У1| = I с^'ЬФу (V, )/с!еЪФд<Х,^) I ,

причем ) 1 в сипу определения "текущего" параметра.

Поэтому

'я ) I ) I с1е1;Фх (X, р ) I > О ,

т.е. матрица Фу при указанном правиле выбора "текущего" параметра обязана быть невырожденной. Если МХ* /{¡¡Л I < МСЗ/с);*! I при = (VI то "текущим" параметром выбирается либо (3 ( если IсЮ/.сЗ^ I > 1 ) , либо (»1 остается "текущим" параметром.

Если Фу — разреженная матрица, то предложенный способ параметризации позволяет учитывать ее структуру при решении систем линейных алгебраических уравнений <2Э> и <30), что существен*" для систем с большой размерность» N . Действительно, после определения "текущего" параметра ? = и к , к=1,2, , или Я = С, для построения Фу достаточно вычеркнуть к-ый столбец матрицы ЕФу, Ф^З , сохранив при этом разреженность матрицы Фу .

Заметим, что в процессе продолжения решения системы <26) параметр 0 может принимать значение О = Су , при котором = О. Однако, если Cq — гладкая пространственная кривая, то рещение системы в окрестности (З^ может быть продолжено по "текущему" параметру ^ , отличному от О в виде решения системы (2В). Зависимость О от в этом случае имеет тот ив вид, что и на рис. 2, Таким образом, параметризация позволяет построить часть гладкой пространственной кривой между двумя особыми точками, если таковые существуют. Длина дуги Ся , отсчитываемая от точки (О, ШО)) до точки <0,и(0)> , принадлежащих Ся , является универсальным параметром системы <26> , каждому значению которого соответствует единственное решение системы (£6), в то время как для заданного

_ ек

-значения О Е С О, □ 3 может существовать несколько решении. При численном построении Cq процедуры адаптации шага по "текущему" параметру ич задания начального приближения решения "параметризованной" системы <28) те же, что и при численном построении поверхности Sq , задаваемой решением краевой задачи 41), <2).

Обратимся теперь к краевой задаче <1>, <2) и рассмотрим вариант приближенного представления краевой задачи в виде системы (26). При этом будем считать, что на каждом из отрезков [х^х^З решение краевой задачи приближенно представлено кубической параболой , определяемой значениями -пункции и ее производной на концам отрээка :

«^Ч^Ы3' 5Э<><г> = Уг»У<х-г>, Зз<х'+4) = У,+<*У<Х^> г г г г г.-"

(х ¿> = Г -- ^(х^, у ,0), ¿53/ЙХ<0«^)= f = >0>-

В середине отрезка кубическая парабола принимает значение у ,

1*1/2 7 (v/ \ 1 + <

у = . 5* (у + у ) + . 125*^*(1= — г > . (31)

Воспользоваваиеь Формулой Симлсона , получим

г 1 + 1/2 141

^ (х,у(х),С!>ах й (Ь>/6)*(Г + 4»? + f

к;

)

в

f^/■г= f<кr.5»h.,ylW/2 ,0) (32)

результате интегрирования системы (1) от х- до х^ приходим к следукщей системе уравнений размерности N=(«1+1 >*п относительно приближенных сеточных значений у"*, у2 у*7**'вектор-функции у(х):

д(у\утЧи) =0

у4 - у* + (Ы/Б)*^-1 + +• fг ) = О , (33)

У - У

® + <ьг/в>*(гг + 4*^2 +• f 3 ) = о ,

уг*« + + г"»* ) = О ,

Очевидно, дискретная модель (31)—(33) позволяет найти приближенное решение с погрешностью порядка Нглах , Нглах = гаах на решениях краевой задачи (1),(£). Матрица Ф^ системы (33) имеет ненулевые блоки размеров п х п только по главной и ниже расположенной блочных диагоналях, а также блок в правом верхнем углу в силу краевого условия (2) . Учет этой структуры позволяет достаточно экономично продолжать решение системы (33) с применением параметризации.

Ту же дискретную модель можно рассматривать как следствие коллонационных условий в средних точках каждого из отрезков Сх^ , , 1 = 1, 2, ■ . -, т,

ЙБ^ /йх = f при х = .5*(х^ + х^ц),

или . . ..

»•»* « 1 гИ 1+4/2

1.5*(у - у1 )/Ь; + .25* (Г + f ) = f

Добавив сюда краевые условия, получим систему (33).

Заметим, что использование коллокационных условий дает возможность сформулировать дискретную модель, аналогичную (31)—(33), нелинейной краевой задачи для системы дифференциальных уравнений, не разрешенных относительно производных, вида :

х е- (а, Ь) , Р<х, у, р, 0) = 0 , д (о( , £ ,0) = 0 , (34)

Здесь р = с1у/с|!< , Р - достаточно гладкая еектор-»уннция по совокупности аргументов в некоторой окрестности решения <34) при условии, что йеЪРр ^ О ■ Другие обозначения сохранили свой смысл. По—прежнему предполагается, что графиком решения (34) является гладкдя поверхность вс^,. Дискретную модель краевой задачи (34) относительно приближенных сеточных значений у1 и рг, у1 = у(х£>, р1 = р(х^), можно представить в виде :

д(у4 ,ут*',0> = О,

Р^.у'.Р'.с» = О

1 = 1,2,...,гл, Р(х-1+Ь1/2,у%*1/г,р1*'^г ,0) = О (35)

¿И г*1 г'+-«'У ' Р = 0

где

» 7-И Г 1*4

у = . 5* (у + у' ) + .125*Ь^*(р' - р 1 )

{-м/2 г-и I г 1+?,

р = 1.5*(у - у + . 25* (р + р' >.

(36)

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

Естественный обобщением <1), (2) является краевая задач?, в

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

л

параметров ч размерности п :

х Е (а,Ь) , dy/dx = f<x,y,q,Q),

(37)

g(o< , £ ,q,Q> = О ,

где теперь вектор-функция g имеет размерность п + п . Очевидно, дискретная модель (37) имеет вид (31)-(33) с V-образной блочной структурой матрицы Якоби за счет производных no q .

Приведем Формулировку дискретной модели краевой задачи (1), (2) с аппроксимацией порядка Нтах', численный анализ иоторой не потребует существенного увеличении вычислительных затрат. С этой целью вместо 5^(х) будем использовать параболу 5-ой степени Sg(x), принимающей на концах отрезка Сх{ , х{»(3 заданные значения вместе с первой и второй производными. Вместо Формулы Симпсона обратимся к интерполяционной Формуле Котеса, точной на многочленах 5-ой степени. При этом дискретная модель запишется в виде следующей системы трансцендентных уравнений :

9«У4,У^,0> = 0 ,

1=1,;

Здесь

» I »+<Л)

у - у + (^/90) *(7*f + Зг*f + (38)

1 + </2 ^Ъ /Ц 7 + 1

+ 12*Г + 3£*f + 7*f ) = о

}=0, 1,2,3,4 Г**'"1 = «■

причем

у14,'Л = (45Э*у'* + 53*у*+1 ) /512 + (Ь£/1024) * (189*Г* - ЗЭ*Г1+1> +

+ (Ь£/2048)*<27»^ ' + -3*ГЧ *f ),

у'М/г = > 5* (у' + у1'+< ) + (5*Ь{/3£)*(Г* - Г1'*« > +

г г 4*4

+ (Ь^/ЗД)»^ у *f + «Г 1 ) ,

г+ЗЛ % г-М * 1-М

у = (53*уг + 459*у > /512 + (^/1024)*(39*Г - 189*1^ ) +

г г 1*1

+ (Ь£/£04а>»<9*^ *f + 27*^ *f ) .

Можно полагать, что на решениях <3в) матрица СФц 3 является

возмущенной матрицей СФ0 , Фф 3 системы (33) с нормой возмущений порядка Иток* и использовать последнюю для вычисления компонент патрицы Ф при решении системы (38). Определенный таким образом модифицированный метод Ньютона позволяет Эффективно строить репе-ния дискретной модели (38) практически с теми же затратами, что и дискретной модели (33).

При численном исследовании дискретной модели важно предусмотреть возможность адаптации первоначально заданной сетки в связи с проблемой аппроксимации решений краевой задачи <1>,<£> в процессе продолжения решения. Воспользуемся правилом распределения узлов сетки при интерполяции функции и (х> € С1*(а, Ь) эрмитовым кубическим сплайном, при котором абсолютная величина погрешности приближения не превышает заданного значения £ . Таковыми являются узлы сетки (27), подчиненные условию :

£ = (Ь£/384)»1й и/ах"(8) I, ВбСх^.х^З.

Пусть теперь и(к) - компонента вектор-Функции у(к), представляющая численное решение краес.-<й задачи (1), (2), полученное на первоначально заданной сетке по х :

а = ( < ... < = Ь , Н^ = - -Ц . (39)

Предположим, что в соответствии с порядком аппроксимации дифференциальной краевой задачи дискретной моделью, абсолютная величии ч ■

на погрешности, связанная с ц(х), пропорциональна Нглах »Ш и/с1хи1, Нглах = гпах Н^. Если построенное численное решение оказалось достаточно близким к точному, то, используя приближенное выражение производной, можно попытаться повысить его точность путем пересчета краевой задачи (1),(2) на вновь Сформированной сетке (27) с варьированными узлами сетки (39) за исключением и - Ис-

ходя из равенств :

£. = »1с|Ч и/с1х" (8) I , Ве[к,,к;н] , <40)

выполняющихся для каждого из отрезков 3 , где хг,х^ ,..

..хт подлежат определению, а С > О - некоторая постоянная, составим сумму

V" 'ц ч </Ч V« ь ми

С * и/ахч <В> I = г(1»е. « С »¿Ш и/йх (х> I *с1х

¡ = 1 а Отсюда следуют равенства, аналогичные (40) !

ц „ V« Ь. Ч ц V«

5 1с1 и/йх4 (X) I *с1х = (1/гп)*$ IЬ и/ах (х) I *с1х (41)

х г а

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

На практике четвертая производная определяется по Б^х). Предполагая сетку (39) достаточно густой, принимается, что распределение погрешности решения краевой задачи достаточно хорошо учитывается кусочно—постоянной Функцией, определяемой на каждом из отрезков 3 четвертой производной, вычисленной, напри-

мер, в средней точке отрезка. При этом графиком функции

1ч ч ^/Ч <5(г) = $1,с1 Ц/ах (х)| *ЙХ , а ( г ( Ь ,

а

будет ломаная линия, а задача построения сетки (27) сведется к

решении достаточно простого уравнения относительного i-oro уэла s

/

S?(х^> ■ <i-l>» су <Ь> /м .

Рассмотренное правило адаптации сетни дополняется еще одним попезным приемом. Если расстояние Хц^ц — х к , к=1 ,£,..., гп, между узлами сетки (27) оказалось больше заданного dof dc< Ь-а , то вводятся дополнительные узлы. При этом меняется значение гл, а. вместе с ним и размерность системы, представляющей дискретную модель краевой задачи. После Формирования сетки <27) решение краевой задачи на сетке <39) используется для задания начального приближения на сетке (27) с целью пересчета решения.

ДИСКРЕТНЫЕ МОДЕЛИ, ОСНОВАННЫЕ Hfl СПЛАЙН-ИНТЕРПОЛЯЦИИ

Рассмотрим частный случай краевой задачи (1),(2), имеющей Формулировку типа (18) s

хе (а, Ь) , dy/dx = В(х)*у + R(x,y,0) ,

(42)

Ly(а) = О , Му(Ь) = О

Здесь В(х) - <пкп) ~ матричная Функция, Rix,у,В) - вектор-Функция размерности n; L - (п х п.») - матрица ранга гц, М - (пьх и) - матрица ранга rit, п . Предпопагается, что однородная краевая задача

х в (а, Ь) , dz/dx = B(x)»z

<42')

L»z(а) = О , R*z(b) = О

имеет только тривиальное решение. Для краевых задач небольшой размерности (г> = 2,3) с достаточно гладкими решениями, позволяющими не обращаться к процедуре адаптации (27) сетки в процессе продоп-жения решения по параметру, эффективными оказались дискретные модели, основанные на приближенном представлении вектор-Функции R(x,y(x),Q) в виде интерполяционного сплайна по х на сетке (27). При этой коэффициенты сплайна будут Функциями приближенных сеточных значений у(х) :

га+1

R(x,y(x),Q) я S (х) = 4K(x)*R(x<,yK, О) ,

к=1

гяе*?к(х>, к=1, 2, . .., т+1, — Фундаментальные сплайн— функции. После замены в (42) вектор-Функции R(x,y(x),G!) на S(x) прибпиженное решение краевой задачи (42) записывается в виде :

гл+1

у (х ) - Ч^к (х>*И(х^,у*,С!) , (43)

к=1

где матричная Функция Ч'к(х) находится из решения линейной краевой задачи : к=1, £,..., гп+1 ,

и 6 (а, Ь) , с! Ч'к / йх = В(х)»Ч/К+ <ЯК(х)*Е ,

(44)

^^«(а) = О , М**^к(Ь) = О

В силу предположения Относительно (42), матричная функция Ч"к(х) определяется при этом однозначно. Лиснретная модель (42) следует из (43), если потребовать выполнения равенства в узлах сетки(27):

1 = 1, 2, . .., га+1,

гя+1

у'х = ф^(х^)*Н(х(4,уК,а) (45)

к=1

Вычисление матриц системы (45) производится в два

этапа. Вначале по заданной сетке определяются фундаментальные сплайн-функции. Как правило, это кубические интерполяционные сплайны, позволяющие на решениях (42) пёлучить аппроксимацию порядка . Значения параметров сплайна *Рк(х) в этом случае находятся из решения линейной алгебраической системы с трехдиаго-нальной матрицей, имеющей диагональное преобладание, то есть вычислительные затраты на этом этапе, еспи используется метод прогонок, весьма незначительны. Столь же эффективно может быть организовано численное решение серии линейных краевых задач (44) на втором этапе, например, методом ортогональных прогонок. ( Естест— венно, что достаточно хорошая обусловленность линейной краевой задачи (44) предполагается). Еспи учесть, что далее предстоит многократное обращение к решению системы (45), связанное с продолжением решения по параметру,.то усилия, затраченные на Формирование дискретной модели предлагаемым способом, не выглядят чрезмерно большими. Особенно это относится к случаю, когда в виде (42) представлена краевая задача с линейным дифференциальным оператором высокого порядка и правой частью, зависящей от небольшого числа компонент у, как это имеломесто в краевой задаче (18).

Те же коэффициенты системы (45) опредепяются через матричную функцию Грина Г(х^) линейной однородной краевой задачи (42'):

Ь

= ^("(х^)*«?,.^)^, (46)

а

При этом дискретная модель (45) сразу следует из интегрального представления краевой задачи (42) после замены вектор-Функции

Я<х,у(х>,0) на сплайн. Иногда Формула (46) позволяет найти точные выражения коэффициентов Ч^к'х^) через параметры сплайна ^(х). Задача упрощается, если использовать интерполяционные сплайны класса С(а,Ь). В некоторых случаях (для получения численных оценок решения) полезно использование кусочно-постоянной вектор—функции 5(х), заменяющей Я(х,у(х),□>. И тогда

«ы

^к(х^) = ^ Г (х^, 1; )с№ «К

В качестве примера рассмотрим дискретную модель краевой задачи (14) на сетке (£7). Поскольку Функция Грина 6(х,1;), введенная в (15), имеет простой вид, то коэффициенты дискретной модели а1К 1 опРеяеляеные по Формуле :

1

О

имеют точные выражения через параметры интерполяционного кубического сплайна 1?ц(х), которое мы приводить не будем. Укажем, что аналогичная ситуация имеет место, например, при Формулировке дискретной модели краевой задачи, часто встречающейся в приложениях:

х £ (0,1) , + (}/х)*с1и/с1х + 0»Р(и) = О ,

х = 0 : йи/йх = О ; (47)

х = 1 : ач/йх +• 50*и = О ,

где ф — 1 , гс > О , Р(и> — заданная функция. Обращаясь к интегральному уравнению (15), получим систему трансцендентных уравнений вида (45) :

икаи(Х|<), к = 1, £,...,т+1,

«1+1 (48)

и» = П* } а1|С *Р(иК) , 1=1,£, . . . ,т+1 . к=1

Очевидно, система <48) будет дискретной моделью краевой задачи <47>, если при вычислении коэффициентов учесть вид Функции

Грина. В случае краевой задачи (14) и„||= О в силу краевого условия. Далее к системе (48) применяется тот же общий подход, позволяющий изучить зависимость решения от параметра 0.

Пусть "текущим" параметром выбрана компонента и^ . Тогда "параметризованную" дискретную модель системы (48) можно представить в виде :

i = l, 2,---- m, , frt = uj ,

m+1 rn+1

ui = .ДдИ. *F<UK) , (49)

k=l k+-l m+1

a = fl /V"a|< *F(uv) . k+1

Именно в таком виде изучалась дискретная модель краевой задачи (18). При этом параметром fl служила компонента и< , а коэффициенты определялись по Формуле :

а1к= SG<Xi't>dt .

агк

где G(x,t) - Функция Грина (18'). В результате из решения системы (49) методом простой итерации находились двустороннние оценк>- "ьз— пряжения срабатывания" пленочного электростатического реле по значениям оценок Q = Q*.

ПРОГРАММНАЯ РЕАЛИЗАЦИЯ ПРОГРАММНЫЕ КОМПЛЕКСЫ "SYSTEM","BPR-Q", "INTENS"

Предложенные диснретные модели краевых задач и вариант метода продолжения решения по параметру послужили основой для разработок комплексов программ, позвопямщих чиспенно исследовать системы нелинейных уравнений (26) (комплекс "SYSTEM") и нелинейные краевые задачи <1>,(2) (комппекс "BPR-Q"). Благодаря достаточно обшей Формулировке проблемы, эти комплексы программ легко адаптирух!тся к различным приложениям, что подтверждается многолетней практикой их использования. Отметим также комплекс программ "NEWSТ", разработанный Лукьяновой Р.Г., ориентированный на чиспенное исследование нелинейных краевых задач для систем обык--ювенных дифференциальных уравнений 2-го порядка, где для построения дискретной модели применялась сплайн-интерполяция правых частей уравнений, обсуждавшаяся выше. Приведем краткие ха-зактеристики "SYSTEM" н "BPR-Q".

1."SYSTEM". Комплекс программ "SYSTEM" ориентирован на чис-шнный анализ решений системы нелинейных уравнений общего вида и фоизвольной размерности, содержащей параметр,с учетом возможной множественности ревений. Данная проблема возникает при определе-(ии стационарных реаений автономных систем уравнений и их устой— швости в зависимости от параметра при математическом моделиро-1энии стационарных режимов в процессах макрокеталиэа, горения и -.д. На основе алгоритмов, разработанных в Институте математики ¡О РАН предлагается вариант метода продолжения решения по параме-•ру (метод Ньютона в сочетании с параметризацией) и анализом ус-

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

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

Начиная с 1384 г., версии "SYSTEM" использовались при решении различных прикладных задач, в частности, в Институте математики СО РАН , Институте катализа СО РАН , НИФХИ им. Карпова РАН (г.Москва) и др. В зависимости от математических моделей, число уравнений менялось от порядка 10 до порядка 1000 (для дискретных моделей краевых задач) . При этом временные затраты определялись не только размерностью системы, но и организацией процесса продолжения решения, учетом структуры системы и т.д. Размерность рассматриваемых систем уравнений ограничена оперативной памятью ЭВМ.

В настоящее время комплекс программ реализован на IBM РС/ЙТ, исполнитель Гайнова И.А. ( операционная система MS-DOS, язык программирования FORTRAN, version 5.0). Результаты счета сохраняются на диске. На экран дисплея выводится краткая информация о решаемой задаче и цветные графики выбранных пользователем компонент решения. Подобраны демонстрационные примеры, иллюстрирующие различные возможности комплекса.

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

2."BPR-Q". Комплекс программ "BPR-Q" ориентирован на численное исследование нелинейных краевых задач для систем обыкновенных дифференциальных уравнений первого порядка общего вида на конечнон

отрезке по скалярному аргументу в зависимости от параметра задачи. На основе алгоритмов, разработанных в Институте математики СО PRH, предлагается вариант метода продолжения решения по параметру (квазилинеаризация в сочетании с процедурой параметризации) нелинейной краевой задачи в предположении, что графиком решения в пространстве компонент вектора неизвестных и сналярыого аргумента является гладкая поверхность. При этом одному и тому же значению параметра может соответствовать несколько решений краевой задачи. Комплекс программ позволяет численно построить гладкую поверхность, выделив решения при заданном значении параметра и в точках ветвления решений типа "поворот", в заданной области изменения параметра. К рассматриваемой пробпеме относятся краевые задачи, моделирующие стационарные режимы работы каталитических реактороз, процессы на зерне катализатора, нелинейные колебания, а также многие другие математические модели, описывающие явления, в которых возможен гистерезис.

В комплексе программ "BPR-Q" исследуется система нелинейных уравнений с параметром, представляющая дискретную модель краевой задачи, решение которой с высокой точностью аппроксимирует решение краевой задачи на сетке по скалярному аргументу с произвольным разбиением (задаваемая точность аппроксимации 2,b,S) . В основе определения дискретной модели — интегральное представление краевой задачи на каждом из отрезков разбиения. В комплекс включена программа адаптации сетки с добавлением узлов, сопровождающей процесс 1 продолжения решения, что позволяет строить решения с большими градиентами. Используемый в комплекса вариант параметризации (выбор параметра для продолжения решения среди компонент вектора неизвестных и параметра системы) сохраняет разреженность матрицы Якоби, что делает метод экономичным. Кроме того, экономичность метода обусловлена выбором варианта продолжения решения, допускающего достаточно большие шаги по параметру. Процедура адаптации шага позволяет избежать авостов, а также предоставляет возможность попасть в достаточно малую окрестности бифуркационного решения, отличного от "поворота".

Версия "BPR-Q" используются для решения прикладных задач, начиная с 1985 г., в частностиj в Институте математики СО PRH, Институте катализа СО РРН, Институте неорганической химии СО РДН и других организациях. Практика использования "BPR-Q" свидетельствует об эффективности комплекса как инструмента численного исследования нелинейных проблем. Применение комплекса имеет естественное ограничение : предлагаемый алгоритм эффективен для построения достаточно гладких решений без сильной осцилляции. Ограничение на размерность рассматриваемых систем дифференциальных уравнений связано с об'емом оперативной памяти ЭВМ. В настяяее время комплекс программ "BPR-Q" реализован на IBM PC/RT, исполнитель Гайнова И.О., с теми же возможностями, что и у версии "SYSTEM".

В последующей версии комплекса предполагается : разработать режим диалога и соответствующий вывод графической информации, модернизировать алгоритм численного дифференцирования. Готовится вариант версии "BPR-Q", основанный на дифференциальной прогонке решения линейных краевых задач на итерациях.

3."INTENS". Как пример эффективного использования идеи параметризации укажем на комппекс программ "INTENS" , разработанный в связи с исследованием параметрической чувствительности решений задачи Коши дпя автономной системы уравнений размерности п на конечном отрезке СО, t03, to>0, по скалярному аргументу t :

t€(0,to>, dy/dt = f(y,Q),

(50)

у = IJ(Q) при t=0.

Здесь f(y,Q) и y<Q) - достаточно гладкие вектор-функции размерности п , □ - скалярный параметр. Пусть y(t,Q) - решение (SO) при некотором В £ CB,Q3, и пусть компонента yn(t,Q> решения достигает мак'симапьного значения f* при t =Л, Л£(0, t0)., т.е. ri-я компонента f(t,Q) обращается в О при t = Аi

уп =¡«1, fn(y, Q) = О при t =А. (51)

Представляет интерес изучение зависимости р) и Л от параметра Q. С этой целью на отрезке СО,t03 строится решение задачи Коши дпя уравнений Пуанкаре

te(0,-t0>. u (t) = ya(t,Q),

du/dt = fM (t,y(t,Q),Q)*u + fQ(t,y(t,Q),Q),

u =<?g(B) при t = О,

и определяется u(Ä > как Функция параметра □. Возможна ситуация, когда в окрестности некоторого значения □ = Q», Q« 6CQ, S], малым изменениям параметра Q будут отвечать бопьвие изменения параметра f, то есть lu(h )1»1, что свидетельствует о появпении численной неустойчивости решений задачи Коши (50). Эта пробпема возникает, например, при интегрировании сингулярно возмущенных систем, описывающих явления типа "теплового взрыва". Выбрав [V» в качестве "текущего" параметра, мы тем самым переходим от задачи Коши на отрезке СО,А ] по t к краевой задаче на отрезке СО,13 по х, х = t/Л , присоединив к (50) условия (51) :

XÉ (0,1), dy/dx = a*f<y,B),

y = <J>(Q) при x - 1, (52)

У = И i fri<y»a> = 0 ПРИ x = 1

где fK — заданный параметр, a Q(f«l ) и Д ( pi ) подлежат опредепению. При этом решение краевой задачи (53), очевидно, будет обладать свойством : IdQ/dfi 1«1. В процессе продопжения решения (53) параметры Q и Л также могут быть "текущими" параметрами.

Предлагаемый комплекс программ "INTENS" является адаптацией "BPR-Q" путем учета специфики краевой задачи (52). В частности, начальное приближение решения при стартовом значении параметра Q определяется из решения задачи Копи (50). В дапьнейшеи построение зависимостей f4(Q) и Д((3) осуществляется по методу продолжение гг.— шения (52) с параметризацией. Кроме того, задача Коши (50) испопь— зуется для построения решения на отрезке по t.

ЗАКЛЮЧЕНИЕ

Таким образом, в докладе представлены 1) разработка варианта метода продолжения решения в сочетании с параметризацией, 2) его программная реализация в виде комплексов программ, позвопяю— щих проводить численное исследование нелинейных краевых задач и систем трансцендентных уравнений в общей Формулировке, и 3) численное исследование конкретных матетематических моделей из различных областей приложений, отраженные в списке публикаций. К ним относятся математические модели пленочной электромеханики, макро-катглиза, отдельные задачи из таких разделов, как нелинейные колебания, теорип горения. Ряд примеров взяты из публикаций ,и результаты их численного исследования представпят интерес в качестве тестов для проверки эффективности применения программных комплексов.

Основные результаты.

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

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

3.Предложены дискретные модели нелинейной краевой задачи, определяемые на сетке с использованием Формул Котеса, а также на основе сплайн-интерполяции.

4. Предложен вариант адаптпции сетки в процессе продолжения решения дискретной подели краевой задачи.

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

8. Разработаны комплексы программ "BPR-Q", "SYSTEM" и "INTENS", прошедшие многолетнюю апробацию на численном исследовании математических моделей из различных областей приложений.

7.В области пленочной электромеханики в рамках математических моделей были найдены основные характеристики работоспособности ряда приборов. Описание устройств и результаты расчетов приведены в монографии С2&3.

8.В области математического моделирования стационарных процессов на зерне катализатора и в химических реакторах были проведены исследования для Института катализа СО РАН математической модели многокомпонентной диффузии при наличии химической реакции в зерне катализатора в процессе синтеза жидкого топлива и химко-технологи-ческой схемы с рециклами, впервые построена бифуркационные зависимости для стационарных решений двухфазовой диффузионной модели кипящего слоя с определением критических значений параметров.( Примеры построения автором бифуркационных кривых, иллюстрирующих множественность решений различных нелинейных проблем, имеются в [243

и (2БЗ.). По математической модели, предложенной в Институте неорганической химии СО РАН, был исследован процесс роста кристаллической Фазы карбида кремния в непроточной газотранспортной ячейке.

Исследования проводились в соответствии с утвержденным планом научно-исследовательских работ Института математики СО РАН. Начиная с 1983 года, тематика исследований была внлючена в план совместных с Институтом катализа СО РЙН работ в рамках программ ГКНТ СССР "Разработать эффективные численные методы расчета периодических решений" и "Методы расчета, оптимизации и управления нестационарными процессами" (под научным руководством С.К.Годунова со стороны Института математики). В частности, в соавторстве были представлены следующие отчеты по этим темам ( номера гос. регистрации 018300055962 и 01830069286 соответственно )s

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

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

этап« Расчет стационарного Фронта реакции.-1985, 60 с. этап: Расчет математической модели ХТС, задача СтеФана.-1985,

165 с.

этап: Численные исследования множественности стационарных решений нелинейных краевых задач методом Ньютона с параметризацией. -1985, 275 с.

Апробация результатов.

Помимо публикаций результаты работ докладывались на ряде конференций и семинарах!

Всесоюзная конференция " Химреактор-6", Дзержинсн,- 1977.

Всесоюзный семинар по исследованию Фронтов химических реакций", Новосибирск, 19S5.

Всесоюзная конференция "йктуальные проблемы вычислительной и прикладной математики". Новосибирск, 1987.

Всесоюзная конференция "Математическое моделирование : нелинейные проблемы и вычислительная математика". Звенигород, 1988.

Рабочий семинар "Численные методы теории бифуркации". Пущино,

1388.

Региональная конференция по математическому моделировач^' i химико-технопогическин процессах. Яремче, 1990.

3-я Всесоюзная шкопа—семинар по математической кинетике, химической и магнитной гидродинамике. Красноярск, 1390.

Школа—семинар по комплексам программ математической физики. Ростов-на-Дону, 1990.

Международная конференция "Нестационарные процессы в катализе " . Новосибирск, 1990.

3 я Всесоюзная конференция "Динамика процессов и аппаратов химической технологии ". Воронеж, 1990.

VII Всесоюзная конференция "Математические методы в химии". Казань, 1990.

Республиканский семинар "Динамина процессов и аппаратов непрерывной технологии : моделирование, идентификация, управление и проектирование ". Яремче, 1991.

Всесоюзная конференция "Химреактор-11". йлушта, 1992.

Международная конференция по электронным материалам. Новосибирск, 1992.

Шкопа семинар по комплексам программ математической физики. Новосибирск, 1992.

В течение последних пяти лет автор выступал с докладами по рассматриваемой теме на семинарах в ИМ СО РЙН, ПК СО РЯН, ИТПМ СО РАН, ВЦ СО РЙН, ИКиГ СО РЙН, ИНХ СО РЙН, ИВТ СО РЙН, ф-тет Вычислительной математики и кибернетики МГУ, НИФХИ им.Карпова(Москва), Институт хим.-физики РЙН(Черноголовка).

Работы по программным комплексам "SYSTEM", "3PR-0", и "INTENS" в настоящее время ведутся в ранках программы ГКНТ по приоритетным направлениям развития химичесной науки и технологии,(" Развитие и применение методов математического моделирования и вычислительного эксперимента для химико-технологических систем и производств", головная организация НИФХИ им.Карпова, Москва).

ПУБЛИКАЦИИ

I. Дятлов В. Л., Фадеев С. VI., Шведова К. В. Некоторые результаты расчета статических характеристик пленочных электростатических реле. В книге Вычислительные системы. Материалы ко 2-ой Всесоюзной конференции. Секция 4.-Новосибирск, 1969. -С. 72-75.

2. Лукьянова Р.Г., Фадеев С.И., Шведова К.В. Расчет статических параметров механической модели пленочного электростатического реле // Вычислительные системы. -Новосибирск, 1970. Вып. 40. -С. 3-35.

3. Фадеев С.И. Численный метод решения одного интегрального уравнения типа Гаммерштейна в связи с расчетом характеристики пленочного электростатического реле// Вычислительные системы.—Новосибирск, 1972. Вып. 50. -С. 3-29.

4. Дятлов В. Л., Кирилюк А. Г. ,Коняшкин В. В., Потапов Б. С. , Фадеев С. И. Исследование пленочных консольных структур// Техника индексаций.-Киев: Наукова думка, 1976.-С. 39-46.

5. Доронин В.П., Лукьянова Р.Г., Фадеев С. И. Численное построение всех решений краевой задачи для системы нелинейным дифференциальных уравнений // Методы сплайн-функций. —Новосибирск, 1977. Вып. 72 : Вычислительные системы. С. 99-114.

6. Доронин В. П. , Лукьянова Р. Г., Слинько М. Г. , Фадеев С. И. Шеплев B.C. Неединственность стационарных режимов при протекании экзотермической реакции при псевдоожиженном слое катализатора.// Тезисы докладов 6-ой Всесоюзной конференции по моделированию химических и нефтехимических процессов и реакторов,"Химреактор-6", Дзержинск,-1977,- т.2, С. 289-298.

7. Фадеев С.И. О численном решении линейной краевой задачи для обыкновенных дифференциальных уравнений методом диФФеренци- • альной прогонни// Методы сплайн-Фуекций.-Новосибирск,1978. Вып. 75t Вычислительные системы. С. 80-95.

8. Фадеев С. И. Об электростатическом притяжении пленочной металлизированной ленты// Моделирование в пленочной электромеханике. -Новосибирск, 1981 . Выл. 84: Вычислительные системы.-С 64-73.

9. Лукьянова Р. Г., Фадеев С. И. , Шеплев B.C. К расчету экзотермической реакции первого порядка на зерне катализатора // Методы сплайн-Функций.-Новосибирск, 1981. Вып.87s Вычислительные системы.-С. 99-113.

10. Фадеев С.И. Расчет характеристик многослойной пленочной емкостной структуры// Моделирование в пленочной электромеханике. -Новосибирск, 1982. Вып. 95: Вычислительные системы.-С 66-78.

II. Дятлов В.Л., Фадеев С.И. Перспективы развития пленочной электромеханики // Проблемы обработки информации. -Новосибирск, 1983. Вып. 100: Вычислительные системы. -С. 129-139.

12. Фадеев С.И. О решении систем трансцендентных уравнений с параметром методом Ньютона. Новосибирск.1984.-25 с.-(Препринт/ / СО ЯН СССР. Институт математики, 72 )

13. Фадеев С.И. О решении систем трансцендентных уравнений с параметром методом Ньютона // Методы сплайн-функций.—Новосибирск, 19В5. Вып.108: Вычислительные системы.-С. 78-93.

1A. Reshetnicov S.I., Fadeev S.I., Sheplev V.S., Methods for 'finding steady-state soltions of mathematical models for chernico-thechnological schemes// React. Kiriet. Catal. Let t. —198Б. v. ЗО, 2,— P. 275-281.

15. Фадеев С.И. О применении параметризации при решении нелинейных нраевых задач// Аппроксимация сплайнами.-Новосибирск, 1987. Вып.121s Вычислительные системы. -С. 102-12А.

1Б. Лукьянова Р.Г., Фадеев С.И. Применение параметризации при чиспенном решении нелинейных краевых задач// йктуапьные проблему вычислительной и прикладной математики. Тезисы докладов Вс^ияа-ной конференции.-Новосибирск, 1987. —С. 182-183

17. Фадеев С. И. Программа численного решения нелинейных краевых задач для систем обыкновенных дифференциальных уравнений с параметром. Ч. 1-3.-Новосибирск, 19S8. -1А2 с.-( Препринт/ СО РАН СССР. Институт математики, 31-33 ).

10. Яятлов В.Л., Фадеев С.И. Задача о поршне в переменном электростатическом поле. Расчет электрокинетического давпения//Мо-делирование в пленочной электромеханике. -Новосибирск, 1989. Вып. 131: Вычислительные системы.—С. 71-87.

19. Фадеев С. И. 0 методе дифференциальных встречных прог онок. // Моделирование в пленочной электромеханике. -Новосибирск, 1989. -Вып. 131s Вычислительные системы. -С. 88-119.

20. Иванов Е.Д., Фадеев С.И. Исследование бифуркаций стационарных решений двухфазной диффузионной модели кипящего слоя// Математическое моделирование каталитических реакторов. —Новосибирск: Наука. Сиб. Отделение, 1989.- С. 106-120.

21. Решетников С.И., Фадеев С.И. Численный метод построения стационарных режимов в зерне катализатора// Математическое моделирование каталитических реакторов. -Новосибирск: Наука. Сиб.Отделение, 1989. - С. 215-232.

22. Фадеев С.И. Численное исследование модели теплового взрыва. -Новосибирск, 1989. -38 с.-( Препринт СО ЯН СССР. Институт математики, 22 ).

23. Фадеев С.И. Численное исследование модели теплового взрыва// Приближение сплайнами.-Новосибирск, 1990. -Вып. 137: Вычислительные системы. - С. 122—1А7.

24. Фадеев С. И. Программа численного решения нелинейных краевых задач для систем обыкновенных дифференциальных уравнений с параметром // Вычислительные методы линейной алгебры. Под редакцией Годунова С.К. - Новосибирск! Наука, Сиб.Отделение. -<Труды Института математики.- 1990. -Т. 17. -С. 104-198 ).

25 Фадеев С. И. Программы численного исследования нелинейных проблем методом продолжения решения по параметру// Нестационарные процессы в катализе. Тезисы докладов. Международная конференция, -Новосибирск.-1990. С. 223-224.

26. Дятлов В. Л. , Коняшкин В. В. , Потапов B.C., Фадеев С. И. Пленочная электромеханика.-Новосибирск: Наука, Сиб. Отделение,1991. -24В с.<4.3. Численный анализ нелинейных краевых задач для систем обыкновенных дифференциальных уравнений. С. 140-238).

27. Иванов Е.й.,Фадеев С.И. Применение параметризации при решении системы нелинейных уравнений с параметром//Хиническая технология.-Киев! Наумова думка, 1991.- С. 95-101.

28. Фадеев С.И. Организация численного эксперимента при иссле-доь^ ии нелинейных краевых вадач // Пленочная электромеханика. -Новосибирск, 1991. -Вып. 143i Вычислительные системы. -С. 78-10Б.

29. Иванов Е.ft., Фадеев С.И. Численное исследование нелинейных моделей химической технологии методом продопжения решения по параметру// 7-я Всесоюзная конференция "Математические методы в химии", Тезисы докладов.-Казань.-1991.-С. 12-14.

30. Березин К), й., Решетников С. И., Фадеев С. И. Математическая модель ХТС как нелинейная краевая задача// 7-я Всесоюзная конференция "Математические методы в химии", Тезисы докладов.-Казань.—1991. -С. 141-143.

31. Иванов Е.Я., Фадеев С. И. Исследование моделей химической технологии методом продолжения решения по параметру. //Доклады Всесоюзной конференции "Химреактор-11",-Харьков.-1992. С. 251-255.