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

0.1. Постановка задачи.

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

Математически проблема исследования решений нелинейной краевой задачи со скалярным параметром Q: dy fa = f{x, y, Q), x?[a, b], g (y{a), y (b), Q) = o, где у, f (x, y, Q) и g (y (a), y (b), Q) — непрерывно-дифференцируемые вектор-функции по совокупности аргументов в некоторой области Q. ставится следующим образом. Пусть известно, что при Q = Qq краевая задача (0.1) имеет решение у (х, Qo). Тогда областью Q является окрестность у (х, Qo). По непрерывности решение краевой задачи (0.1) существует и в окрестности Q. Таким образом, ставится вопрос об изучении зависимости решений y (x, Q) нелинейной краевой задачи (0.1) при значениях Q из области Q. В дальнейшем мы будем предполагать, что графиком решения у (х, Q) в (7V+ 1)-мерном евклидовом пространстве (х, у) является гладкая поверхность Eg. Отметим, что при этом не исключаются случаи, когда некоторым значениям параметра Q отвечает несколько пространственных кривых, принадлежащих Eg, которые являются графиками решений краевой задачи (0.1). В силу гладкости Eg это означает, что краевая задача (0.1) может иметь ветвление решений только типа «поворот».

Как правило, невозможно проверить предположение гладкости поверхности Eg. Кроме того, далеко не всегда известно, что решение краевой задачи (0.1) существует при всех значениях Q из окрестности U. Утверждается лишь, что в силу непрерывности существует некоторая окрестность точки Q = Qo, принадлежащая О, где вектор-функция y (x, Q) может быть построена. В то же время, значение Q не может, очевидно, выходить за область определения вектор-функций f (x, y, Q) и g{y{a), y (b), Q) по Q. Поэтому численное исследование свойств у (х, Q) носит характер численного эксперимента.

Задача определения решения нелинейной краевой задачи (0.1) при Q = Qq, вообще говоря, является неформализуемой проблемой: каждый раз конкретная нелинейная краевая задача требует индивидуального подхода для поиска решения. Существуют различные способы попыток построения решения при Q = Q0. Отметим два способа.

1. Пусть требуется найти решение нелинейной краевой задачи: = хе[а, ь], ^ д (у (а), у{Ь)) = О методом Ньютона. Одним из способов решения этой проблемы является погружение (0.2) в однопараметрическое семейство краевых задач путём ввода искусственного параметра (метод гомотопии) [15]: = F (x, y, Q) = 0, хе [a, b], ^ ^.

GX (y (a), y (b), Q) = 0.

При этом всегда вектор-функции F и G можно задать таким образом, что.

F (s, у, 1) = /(*, у), G (y (a), y (b), 1) = д (у (а), у (Ъ)).

Вектор-функции F и G называются операторами гомотопии. Если при Q — 0 задание начального приближения в методе Ньютона для краевой задачи (0.3) не вызывает затруднений, то можно попытаться найти начальное приближение для решения краевой задачи (0.2) продолжая решение (0.2) от Q = 0 до.

Q =.

2. Пусть параметр Q в краевой задаче (0.1) является одним из параметров математической модели и требуется численно построить зависимость решения задачи (0.1) от параметра Q на отрезке [Qi, Q2]. Однако, далеко не всегда известно начальное приближение решения при Q — Q. При этом возможны две ситуации. a) Известно, как можно задать начальное приближение решения краевой задачи (0.1) при Q = Qo, Qo [Qi, Q2- Тогда выполняется попытка продолжить решение краевой задачи (0.1) на отрезок [Qi, Q2] от Q = Qq. b) Такой возможности нет. В этом случае предпринимается попытка погрузить решение краевой задачи (0.1) при Q — Q в однопараметрическое семейство введением искусственного параметра, о котором говорилось выше.

Отметим, что, как правило, априори заранее неизвестно, чем завершится продолжение решения по параметру по методу гомотопии. Например, решение краевой задачи (0.3) может не существовать при Q = 1.

Продолжение по параметру — это способ численного изучения зависимости решения краевой задачи от параметра по методу Ньютона, использующий производные по параметру для задания очередного хорошего начального приближения решения. Пусть при Q = Q1 известно решение у (х, Qi) краевой задачи (0.1). В этом случае имеется возможность определить производную v (x, Qi) = yQ (x, Q 1) решения y (x, Q{) по параметру Q при Q = Qi из следующей краевой задачи: = A (x)v + R (x), х? [a, b], ^ ^.

Sv (a) + Tv{b) + <�р = 0, где.

Aix) = R (x) = fQ (x, y (x, Qi), Qi),.

S = 9a (y (a), y (b), Qi), T = gp (y (a), y (b), Qi), p = 9Q{y{a), y (b), Q i), a = y (a), P = y{b).

Выберем достаточно малый шаг AQ по параметру Q и примем за начальное приближение y^ix, Q + AQ), например, вектор-функцию:

Ql + AQ) = у (х, ОД + AQ • v (x, ОД, (0.5) чтобы определить решение краевой задачи (0.1) при Q — QIAQ и так далее. При этом обычно вводится ограничение на число итераций, за которое определяется решение нелинейной краевой задачи (0.1). Если число итераций достигло заданного числа щ, то начальное приближение можно улучшить за счёт уменьшения шага по параметру, например, в два раза. Тогда начальное приближение выбирается в виде.

УЩ (a, Qi + = У (*, Qi) +^-v Or, Яг).

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

Как правило, нелинейная краевая задача характеризуется множественностью решений, т. е. одному и тому же значению параметра Q может соответствовать несколько решений краевой задачи (0.1). В качестве примера рассмотрим краевую задачу d2y.

2/(0) — 0, у{1) = 0.

0.6) с точным решением: у (х) = 2 In ch b Q.

2 б2 ch&a-' ^ ch26″ Введём уо = y (0) = 2 In chb. Тогда (см. рис. A):

QM = l Ь + 2 In (1 + Vl — e-w) «.

Как видно из рис. А, существуют значения параметра Q, при которых краевая задача (0.6) имеет два решения, есть значения Q — при которых нет ни одного решения, и есть особая точка Q = Q* — точка «поворота», которой соответствует одно решение.

Уо 4.

2 .

Q* «0.88 Q.

Q* 2.

Рис. ^4.

Перепишем краевую задачу (0.6) в следующем виде: dyi dx dy2 dx 2/2, же [0,1], = -Qe*.

0.7).

2/2(0) = 0, 2/i (1) = 0.

При Q = 0 краевая задача (0.6) или (0.7) имеет только тривиальное решение. Как было отмечено, в этом случае при Q = 0 можно определить производную решения по параметру Q из краевой задачи, аналогичной (0.4): dvi г ,.

— r = v2, хе 0,1, ах = еи (*) деЫ*)п = -1, (°-8) dx.

2(0) = 0, vi (l) = 0. Из (0.8) получаем единственное решение х2 1 = «у + 2' =.

Следовательно, решение краевой задачи (0.6) можно продолжить по параметру Q. В процессе продолжения решения по параметру Q от 0 до Q* с приближением значения Q к неограниченно растут по абсолютной величине элементы производной uq{x1 Q): yq (x, Q) ->¦ oo при Qt Q*.

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

Возможности метода существенно расширяются за счёт использования параметризации — процедуры поиска параметров, которые, наряду с Q, могут быть использованы для продолжения решения. Если в окрестности особой точки Q* выбрать в качестве параметра продолжения, например, уо, то кривая Q (yo), а также и решение у (х, Q (yo)) строятся без затруднений. Другими словами, удобнее рассматривать следующую краевую задачу с параметром V = 2/0: d2y dx2 dQ QeV = 0, х£[0,1], = 0, dx y'{0) — 0, у (1) = 0, y{0) = ft, которая имеет то же решение, что и краевая задача (0.6), но, в отличие от неё, однозначно разрешима при любом значении параметра /i.

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

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

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

0.2. Метод множественной стрельбы для решения нелинейных краевых задач.

Как известно, основными методами решения нелинейных краевых задач являются различные варианты метода Ньютона (квазилинеаризации), метод стрельбы, метод множественной стрельбы.

Опишем метод Ньютона (квазилинеаризации) [5] для решения нелинейной краевой задачи, в которой формально параметр не указан и, следовательно, не ставится задача об изучении решения краевой задачи от параметра: д{у{а), у{Ъ)) = 0.

Здесь предполагается, что вектор-функции f (x, y) и д{у{а), у (Ь)) — достаточно гладкие по совокупности своих аргументов. Предположим, что известно достаточно хорошее начальное приближение решения у{х) краевой задачи (0.9). Построим последовательность вектор-функций у^(х), к = 1,2,., определяемую из решений линейных краевых задач: = + х € [а, 6], (оло).

5И3/[*+Ч (а) + ТЦу^{Ъ) =.

• А^М = «ж,/^)), FW (x) = /(ic.yWW).

S[k] = 9a (y[k](a)Jk](b)), TW = ^(/](а), У[*](Ь)), ^ =^yW (a) +TVW — 9(у[к]Ш[к]Ш a = s/(a), /5 = #). В дальнейшем будет показано, что при достаточно хорошем начальном приближении у^ (ж) и хорошей обусловленности линейных краевых задач (0.10) последовательность у^(х) сходится к решению у (х) нелинейной краевой задачи (0.9).

Методы коллокации [19] заключаются в разложении неизвестного решения у (х) в ряд по некоторой системе известных базисных функций hj{x): 0 m.

Ук{х) = 11ац^{х), k=l, 2,., N. (0.11) з=i.

Разобъём отрезок [a, b] на m частей, а = хi < X2 < ¦ ¦ ¦ < хт+ = 6.

Потребуем, чтобы вектор-функция (0.11) удовлетворяла дифференциальному уравнению (0.9) в точках к = 2,3,., т. Кроме того, потребуем, чтобы выполнялись краевые условия. В результате получаем систему из N (m — 1) + N = Nm нелинейных уравнений для Nm неизвестных постоянных а^. Из полученной системы методом Ньютона определяются коэффициенты ay-разложения (0.11) и, следовательно, определяется решение краевой задачи (0.9). В зависимости от выбора базисных функций hj (x), j = '1, 2,., m, и расположения точек к = 2,3,., т, получаются различные варианты метода коллокации.

В [27] описан метод ортогональной коллокации. Здесь в качестве базисных функций hj{x) выбраны ортогональные полиномы Якоби степени j, определённые на отрезке [а, Ь]. Разобъём отрезок [а, Ь] на т частей и на каждом из отрезков [xi, xi+i] построим полиномы Якоби Pj (x) степени j, j — 0,1,2,., п. Далее, решение краевой задачи (0.9) будем искать в виде п а-е[ач, ач+1], k = l, 2,., N. (0.12).

3=О.

Потребуем, чтобы вектор-функция (0.12) удовлетворяла дифференциальному уравнению (0.9) во всех точках отрезка [2^,^+1], которые являются нулями полинома Якоби Р*, степени п. определённого на отрезке [х^, Xi+i. и краевым условиям. Эти условия состоят из Nnm + N уравнений. Кроме того, необходимо выполнение условий непрерывности вектор-функции (0.12) во всех точках отрезка [а, 6], т. е. должны выполняться условия y[(xi+l) = у{+1(хш), г = 1,2,., т'~ 1, к = 1, 2,., N.

Таких условий N (m — 1). Следовательно, получаем систему из Nnm + TV + N (m — 1) — N (n + 1) m уравнений для N (n + 1) m неизвестных постоянных.

4 у.

В [22] рассмотрен метод сплайн-коллокации. Представим решение краевой задачи (0.9) в виде системы интегральных уравнений X у (х) = y{xi) + Jf (t, y (t))dt, xi.

Xj.

Сеточные значения вектор-функции y (x) удовлетворяют равенствам:

Xi+l y{xi+i) = y{xi) + J f (t, y (t))dt, i = 1,2,. ., m, ^ ^.

9(y (xi), y (xm+i)) = 0.

Заменим интеграл в (0.13) приближённо на квадратурную формулу, например, по формуле Симпсона. С погрешностью порядка /г5 имеем:

Xi+l 7 fit, y (t))dt ъ ^ [f (Xi, у') + 4 f (xt, f) + f (xi+1, yi+1)] ,.

CCj, где Xi = - hi = Xi+1 — Xi, у1 = у (xi). Для определения у1 воспользуемся кубической параболой Эрмита, аппроксимирующей у (х) на отрезке [xi, xi+i с погрешностью /г, 4: х = хi + hiT, 0 < т < 1. у (х) «8{х) = (1 — т) у* + туш + т (i — т) [(1 — т) а1 + тЩ, где а* = у'- yi+1 + hif (Xil у% V = yi+1 — у{ - hif (xi+1, yi+1). Полагая здесь г = 0.5, получим: yi = y (xi) = 0.5(yl + yi+1) + 0.125h.

1(хиуг) — f (xi+1,yi+1).

0.14).

После замены интеграла в (0.13) на квадратурную формулу Симпсона приходим к следующей системе уравнений:

У'.

Ум + Т f (xi, У1) + 4 f{xi, f) + f{x i+i, yi+l) = о, г = 1,2,., ш (О'15).

9ivym+l) = 0.

Размерность этой системы равна N (m + 1) и она определяет компоненты векторов ук, к — 1, 2,. ., т + 1.

Формулы (0.14), (0.15) могут быть получены из условий коллокации в средней точке. Действительно, в силу выбора параметров, а Ьг, из условий коллокации в средней точке следует формула (0.15). Условия коллокации в средней точке имеют вид.

Xi + Xi+i (Xi + Xi+1 {Xi+Xi+i\ или s'(xi) = /(хг, уг). Отсюда получаем следующие равенства: 1.

— У1 + yi+1 + (1 — 2г) ((1 — т) а1 + тЪг) + r (1 — т)(-аг + 6г) s'(xi) 1.

Л".

А:

2/ + + 4 (2yi+1 — - Гц (f (xi+ь 2/*+1) + /(Ж, 2/0) 1 hi l (yl+1 -y') -(f (*i+byi+1) + 1(Хг, У1) f{Xi, yl). = f (Xi, t).

Последнее равенство совпадает с формулой (0.15): hi yi+1 ~ У1 6 0.

Заметим, что из формул (0.14), (0.15) следует, что данный метод сплайн-коллокации можно использовать для краевых задач с разрывными правыми частями.

В результате метод сплайн-коллокации приближённо преобразует краевую задачу к системе нелинейных уравнений с разреженной структурой матрицы Якоби. При этом эффективность метода Ньютона зависит от учёта этой структуры.

В отличие от метода коллокации метод стрельбы преобразует нелинейную краевую задачу к задаче Коши для системы уравнений (0.9). При этом вектор начальных данных определяется как решение системы нелинейных уравнений, определяемой краевыми условиями. Рассмотрим задачу Коши для уравнения (0.9): g = /(*,"), *е.М,.

У = Р, х = а решение которой обозначим через у (х, р). Здесь р — вектор начальных данных, компоненты которого требуется «пристрелять» таким образом, чтобы выполнялись краевые условия (0.9), т. е.

Ф (р) = д (р, у (Ь, р)) = 0. (0.17).

В этом случае, очевидно, вектор-функция у (х, р) будет давать решение краевой задачи (0.9), которое, как мы предполагаем, существует. Следовательно, проблема сводится к решению системы нелинейных уравнений (0.17) относительно компонент вектора р. Применение здесь метода Ньютона — один из вариантов корректировки «стрельбы» [16, 19, 30].

Пусть р° — начальное приближение в методе Ньютона для решения системы (0.17). Тогда для уточняющей р° поправки Ар0 имеем систему линейных алгебраических уравнений с матрицей Фр (р°):

Ф p (p°)V = -Ф (р°), где.

Фр (Л = 9а (р°, у (Ь, Р0)) + 9р (р°, у (ЬУ))уР (Ь, р0). Затем за начальное приближение берётся вектор р1 = к нему ищется поправка Ар1 из системы линейных алгебраических уравнений.

Ф^р1) V = -Ф (р1), и так далее. Поправка Арк к приближению рк находится из системы линейных алгебраических уравнений.

Фр (рк)Арк = -Ф (р*), к = 1,2,., где рк) = 9а (ру (Ь, рк)) + др (ру (Ъ, рк))Ур (Ь, рк1 (0.18) так что рк+1 = рк + Арк. При этом предполагается, что решения задач Коши.

J = /(*,"), *еМ,.

У =р к = 0,1,2,., V ' > х = а определены на [а, Ъ) и, кроме того, р° — достаточно близкое к решению системы (0.17) начальное приближение.

Как следует из (0.18), для организации итераций по методу Ньютона требуется вычисление при х = Ь матрицы производных вектор-функции у (х, рк) по компонентам вектора начальных данных. Поэтому задача Коши (0.19) решается совместно с уравнениями для вариаций, определяющими матричную функцию V (x) = ур (х, рк): = fy (x, y (x, pk))V, xe[a, b], V /. х = а.

После того как компоненты вектора р будут найдены, решение краевой задачи (0.9) может быть восстановлено на всём отрезке [а, Ь] с помощью решения задачи Коши (0.16).

Отметим, что применение метода стрельбы может быть связано с вычислительными проблемами, на которых мы в дальнейшем остановимся. Эти проблемы во многом решаются с помощью метода множественной стрельбы [19, 30]. Приведём краткое описание метода. Разобьём отрезок [а, Ь] на т частей: а = хг < х2 < ¦ ¦ ¦ < хт < xm+i = Ь, и на каждом из отрезков [xi, xi+i] рассмотрим задачу Коши.

T~ = f (x, y), xe[xi, xi+1], dx м (0.21) у =pM, г = 1, 2,., ra, 1 ;

00 где рМ — вектор начальных данных. Одновременно строится решение задачи Коши У{х, рЩ для матричного уравнения: у- = fy{x, y (x, p[l]))V, xe[xi, xi+1], у -" = j (°-22).

00 ——• 00 'i где у (х, рЩ — решение задачи Коши (0.21). Очевидно, вектор-функции у (х, рЩ, i — 1,2,., т, будут давать в совокупности решение краевой задачи (0.9) при выполнении условий, требующих непрерывности решения в точках разбиения отрезка [а, 6], и краевых условий: ф[1] = ^(рМ р^Ч) = о, фИ=у (®2,рМ)-рИ = 0, (0 23) ф[т+1] = у[хт+ър[т]) — р[т+1] = Q.

В результате мы имеем систему нелинейных уравнений относительно компонент векторов., p[m+1l. Обозначим через р и Ф векторы, составленные из векторов рМ, pl2. и фМ} фИ,., ф[т+Ч соответственно, что позволяет представить систему (0.23) в векторном виде, аналогичном (0.17):

Ф (р) = 0. (0.24).

Итак, в отличие от метода стрельбы, в методе множественной стрельбы «пристреливаются» компоненты га+1 векторов, составляющих вектор р. При этом «стрельба» реализуется с помощью метода Ньютона применительно к системе (0.24) точно так же, как в случае системы (0.17). При этом на итерациях используется матрица Якоби, определяемая на решениях задачи Коши (0.22), которая в данном случае имеет вид.

9а (р[1р[т+1]) др (р[1р[т+1]).

V{x2,pW) -1.

Л, ,, рМ) -I р).

— I.

V (xm+hpM) -I.

После того как компоненты вектора р будут найдены, решение задачи (0.9) может быть восстановлено на всём отрезке [а, Ь] с помощью решений задач Коши (0.21). Можно утверждать, что метод множественной стрельбы позволяет корректно решать хорошо поставленные краевые задачи.

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

0.3. Методы решения линейных краевых задач.

В методе Ньютона (квазилинеаризации) для нелинейной краевой задачи (0.9) на итерациях решаются линейные краевые задачи. По этой причине приведём краткое описание основных методов решения линейных краевых задач.

Рассмотрим линейную краевую задачу: dx.

T = A (x)y + f{x), х G [a, b], (0.25).

0.26).

Ly (a) = ?, Ry (b) = r, гдеy (x), f (x) — вектор-функции размерности N-?rt — векторы размерности к и N — к соответственно, a A, L, R — матрицы размерности N х N, к х N, (N — k) xN соответственно. В дальнейшем предполагается, что ранг матрицы L равен к, а ранг матрицы R равен N — к.

Простейшим по форме методом решения краевой задачи (0.25), (0.26) является метод стрельбы [4]. Найдём частное решение уо (х) неоднородного уравнения = A (x)y0 + f (x), удовлетворяющее левому краевому условию. Затем найдём полную систему линейно независимых решений однородного уравнения t принимающих на левом конце отрезка [а, Ь] нулевое значение. Обозначим эти векторы через yi (x), У2(%), ¦ ¦ Уы-к{%) — Тогда, по принципу суперпозиции, любой вектор

N-k.

У{х) = У0(®-) + Е СгУг{х) «=1 удовлетворяет уравнению (0.25) при любых постоянных Q, i = l,., iV — к, ж левому краевому условию, но, вообще говоря, не удовлетворяет правому краевому условию. Чтобы найти искомое решение краевой задачи (0.25), (0.26) необходимо выделить из полученного многообразия решение, удовлетворяющее правому краевому условию. Таким образом неизвестные постоянные q, i = 1, 2,., N — к, определяются из системы уравнений.

Ry (1) = R.

N-k.

2/0(1) + Е <�чу%(Х) i=1 г.

В предположении однозначной разрешимости задачи (0.25), (0.26) определитель этой системы отличен от нуля. В [19, 30] рассмотрены другие варианты метода стрельбы.

В [7] был предложен эффективный метод решения краевой задачи (0.25), (0.26), получивший название ортогональной прогонки С. К. Годунова. Разобъём отрезок [a, b] на достаточно короткие участки, а = xq < xi <. < хп = Ь, длина которых должна быть тем меньше, чем больше разница между величинами различных собственных значений матрицы А. Выполним интегрирование уравнений для уо, ух,., уы-к на отрезке [xq, xx. Полученные в точке х = х векторы проортогонализируем. На следующем отрезке будем интегрировать уравнения, для которых начальными данными являются векторы, полученные в результате предыдущей ортогонализации. Значения решений на правом конце второго отрезка опять проортогонализируем и т. д., пока не дойдём до точки х = Ь. Проводимые ортогонализации мешают «сплющиванию» решений. Метод позволяет вычислить значения у (х) во всех точках х^ i = 0,1,., п.

В [26] предлагается в каждой точке х = Х{ определять углы между векторами у, у2, уы-к и производить ортогонализацию, если углы становятся очень малыми:

Параметр, а может изменяться от 0 до |. Если, а = 0, то ортогонализацию необходимо производить на каждом шагу, если, а = то ортогонализация не требуется. Критерий (0.27) требует много дополнительных вычислений, поэтому там же предлагается более простой критерий: где М — некоторое постоянное число.

В [13] была предложена модификация метода ортогональной прогонки, в которой ортогонализация по Грамму-Шмидту заменена на процедуру разложения матрицы в произведение двух множителей, один из которых имеет ор-тонормированные столбцы, а другой — верхний треугольный. Как известно из линейной алгебры, такая процедура, носящая название <5^-разложения, обладает повышенной устойчивостью к погрешностям округлений при реализации метода на ЭВМ. В такой модификации ортогональной прогонки доказана её устойчивость к вычислительным погрешностям в предположении хорошей обусловленности решаемой краевой задачи.

В [21] рассматривается метод встречных дифференциальных прогонок. Как известно [6], для существования единственного решения задачи (0.25), (0.26), непрерывно дифференцируемого на отрезке [а, Ь], необходимо и достаточно mm arccos hj.

0.27) max.

Уг, Уг) > M. выполнения условия: det.

LY (a) RY (b) т^О, где Y (x) — фундаментальная матрица решений уравнения: dy dx.

А (х)у.

0.28).

0.29).

Рассмотрим на отрезке [а, 6] две задачи Коши относительно матричной функции иа (х) размерности к х N и матричной функции щ (х) размерности (N — к) х N: dua (x) иа{х = L, а (0.30) dx dub{x) —щ{х)А (х), щ (х) I = R. ах х = Ъ.

Определив из (0.30) матрицы иа (х) и щ (х), составим квадратную и (х) и транспонированную ит (х) матрицы: и (х) иа (х) щ (х) ит (х) = Ua (x) ul (x).

0.31).

Очевидно, матрицы и^(х) и и^{х) являются решениями задач Коши: duTa (x) dx duj{x) dx.

— Ат{х)ита, ul (xl = L.

I x = a.

— AT (x)ul, ul (x) |.

I x = b RT.

0.32).

Пусть UT (x) — фундаментальная матрица решений уравнения: dy dxАт (х)у.

0.33).

V,.

Представим решения задач Коши (0.32) в виде: ита (х) = ит{х) [ит (а)}~ чть (х) = ит{х) [UT{b)YlRT.

Отсюда находим, что ua{x) = LU~a)U (x), щ{х) = RU~l{b)U (x), или и =.

U (x).

LU~l (a) RU~b).

Далее заметим, что U{x) обращает в тождество уравнение: dU dx.

— UA (x).

Следовательно, матрица U~l{x) — является фундаментальной матрицей решений уравнения (0.29), т. е. U (x) = У1(ж) (достаточно продифференцировать тождество U~1(x)U (x) = /). Теперь матричную функцию и (х) запишем в виде: и{х).

LY (a) RY (b).

Y-x).

При выполнении условия (0.'28) имеем: detu (t) ф 0, следовательно, и{х) фундаментальная матрица решений уравнения (0.33). Будем искать решение уравнения (0.25) в виде: у{х) = и" «1(ж)г-(а:) = Y (x)v{x), где.

Y (x) = Y (x).

LY (a) RY (b) фундаментальная матрица решений уравнения (0.29). Согласно методу вариаций произвольных постоянных, вектор-функция 'и (ж) должна удовлетворять уравнению: или = «(*)/(*). ¦ (0.34).

Рассмотрим теперь две задачи Коши на отрезке [а, Ь] относительно вектор-функции va (x) размерности к и вектор-функции vb (x) размерности N — к: dl- = Ua (x)f (x), va (x dx dvb dx L x = a ub (x)f (x), vb (x.

0.35) r.

Очевидно, составной вектор v (x) vb (x) обращает в тождество уравнение (0.34). В результате решение краевой задачи (0.25), (0.26) определяется как решение системы линейных алгебраических уравнений: и{х)у{х) = v (x). (0.36).

Действительно, по построению матрицы и{х) и вектора v (x) из решений задач Коши (0.30) и (0.35), решение системы (0.36) обращает уравнение (0.25) в тождество. Кроме того, из представления (0.36) в виде: иа{х)у (х) = va (x), иъ (х)у{х) = vb (x), а также в силу начальных условий задач Коши (0.30) и (0.35), сразу следует выполнение краевых условий (0.26).

Удобно объединить иа (х) и va (x), представив прогонку в виде одной задачи Коши для однородной системы. Имеем: прямой правый ход прогонки: т d dx и v, т, а Т.

АТ (Х) 0.

TW 0 и V.

LT f.

0.37) х а.

Матрица системы содержит (N + 1)-ый нулевой столбец. Аналогично для вычисления щ{х) и уь (х) из решения задачи Коши определяется прямой левый ход прогонки: d dx г т Щ ' Ат (х) 0″ Г Т «Ч Г Т 1 Щ — RT~ т. Vb. т. vb. J т vt rT х = b.

0.38).

По завершении встречных прогонок (0.37), (0.38) остаётся найти решение у (х) краевой задачи (0.25), (0.26) из системы линейных алгебраических уравнений (0.36).

Следующий вариант метода прогонки [2] применяется при решении многоточечных краевых задач для уравнения (0.25). Пусть в точках Х±, Х2, ¦ ¦ Хщ-1″ xm отрезка [а, Ь] заданы условия: p^y{xs) = а8, det (.

1,2, m.

0.39).

Здесь сps — матрица размерности N х ksois — вектор размерности ks, s = 1,2,., m, ki + к2 + • • • + km = TV. Рассмотрим условие: det (.

Известно, что если определить на [а, Ь] функцию — решение уравнения удовлетворяющее условию то функция удовлетворяет уравнению и условию.

Фг’Ог) = <�Рг, 1г{х) = Ф f{x)y{x).

7г' = Ф?(*)/(®-).

0.40) (0.41.

0.42).

7ifa) = ai. (0.43).

Таким образом, для получения решения многоточечной задачи (0.25), (0.39) достаточно решить задачи Коши (0.40), (0.41) и (0.42), (0.43) для каждого из условий (0.39) и обратиться к системе алгебраических уравнений ФП®-)" «71 (ж) «.

Ф т{х)-у{х) = ¦у (х) = Фт (х).. 7m М .

Эти формулы часто очень неудобны для практической реализации, так как в ряде случаев при численном решении задач (0.40), (0.41) и (0.42), (0.43) получаются сильно растущие решения и, что ещё хуже, столбцы матрицы Ф становятся почти линейно зависимыми.

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

Ф- = -AT (x)^i + Ф<(Ф?'Ф0−1ФГАГ (Я-)Ф<) (0.44) удовлетворяющую условию ф i (Xi) = <�Рг- (0.45).

Прямой подстановкой убеждаемся, что функция.

AM = *?(х)ф) удовлетворяет уравнению и условию.

Я = А{х)ЩЪ? ФО^А + ф Jf (x).

0.46) (0.47).

Теперь решение многоточечной задачи (0.25), (0.39) можно определить из системы уравнений: ФТМ «Ш '.

Фт (ж) • у (х) = ¦у{х) = * An (гс) .

Легко проверить, что.

ФТФ)' = о, и, следовательно,.

Эти соотношения позволяют считать задачи (0.44), (0.45) и (0.46), (0.47) более удобными для расчётов, чем задачи (0.40), (0.41) и (0.42), (0.43).

Замечание. Целесообразно перенести краевое условие из точки х в точку Х2, объединить это условие с условием, уже имевшимся в Х2, перенести полученное таким образом условие в х^ и т. д., а не переносить каждое условие независимо.

Отметим ещё один из методов прогонки для нахождения периодических решений уравнения (0.25) [3]: у{а) = у (Ь). (0.48).

Пусть известна какая-либо m-точечная задача г = 1,2,., т. (0.49).

Здесь Si — матрица размерности ki х N, ^ — вектор размерности hi, ki + + ¦•• + km = N, такая что при заданных Si, S2,., Sm задача (0.25), (0.48) хорошо обусловлена при любых? i, £2,£т. Зафиксируем эти Si, S2, Sm и будем искать £,£2,—Лт такие, чтобы выполнялись условия (0.48). Этого достаточно для решения поставленной задачи. Прогоним каждое из условий (0.49) по всему отрезку [а, Ь]. При этом мы получим Si (x) в каждой точке, а так как £{(х) нам не были заданы, то возьмём общее решение для каждого из? г (х). Отметим, что для? i (x) в известных методах прогонки выписываются линейные уравнения. Мы предположим, что задача (0.25), (0.49) хорошо обусловлена при любых? i, ?2? An5 поэтому решение уравнений для каждого? i (x) будет численно устойчивым. Общее решение для £{(х) содержит ki произвольных постоянных, всего подлежит определению к + +. + кт = N констант ci, С2,., сдг. Пригнав каждое из условий (0.49) в точку х = а, мы получим N линейных соотношений между величинами 2/1 (а), у2(а), ., ум{а>), ci, c2, сдг. Пригнав каждое из условий (0.49) в точку х — Ъ, мы получим N линейных соотношений между величинами у (Ъ), У2 (Ь),., уы (р), с, С2,., cn-Вспомнив (0.48) видим, что мы получили 2N линейных соотношений между 2N величинами. Если предположим, что задача (0.25), (0.48) была хорошо обусловлена, то мы должны получить хорошо обусловленную алгебраическую задачу решения линейной системы порядка 2N.

Замечание. Высказанные соображения остаются, очевидно, в силе и при решении более общей, чем (0.25), (0.4−8), краевой задачи. Указанным приёмом можно решать многоточечную задачу, в которой линейно связаны значения у в различных точках. Если вместо (0.25), рассматривается однородное уравнение ух) = А{х,)у (х), где, А зависит от параметра А, который выбирается таким образом, чтобы задача (0.25), (0−48), имела нетривиальное решение, то пристрелку по, А ведём таким образом, чтобы упоминаемая выше линейная алгебраическая система порядка 2N (в этом случае она будет однородной) имела нетривиальное решение.

0.4. Краткое содержание работы.

В главе 1 рассматриваются основные определения и свойства нелинейной краевой задачи для системы из N обыкновенных дифференциальных уравнений со скалярным аргументом х вида = /(*,!/), х?[а, Ь], (а50) g (y (a), y{b)) = 0, где у, f и # — векторы с компонентами yi, fi и gi соответственно, г — 1, 2,., N. Вектор-функции f (x, y) и g{y{a), y{b)) предполагаются достаточно гладкими по совокупности их аргументов, так что краевая задача (0.50) по предположению имеет решение у = у{%), дифференцируемое как по х, так и по параметрам, входящим в определение вектор-функций f (x, y) и g{y (a), у{Ь)). При этом предполагается, что краевая задача (0.50) хорошо обусловлена.

Будем называть нелинейную краевую задачу (0.50) хорошо обусловленной (хорошо поставленной), если хорошо обусловлена соответствующая (0.50) линейная краевая задача du = A{x)u + F{x), xG[a, b], ^^.

Su{a) + Tu (b) = tp, где A, S и Т — матрицы Якоби:

А (х) = fy (x, y (x)), S = ga (y (o), y (b)), Т = gp{y (a), y (b)), у (х) — решение краевой задачи (0.50), a F{x) и ср — произвольные правые части, F (x)? С[а, ь]- Для удобства здесь введены обозначения векторных аргументов: а = у (а), /3 = у (Ъ). Заметим, что, по непрерывности, вместо у (х) в формулировке краевой задачи (0.51) можно взять любую другую функцию у (х) из достаточно малой окрестности решения у (х) краевой задачи (0.50), которое, по предположению, существует. При этом требуется, чтобы «возмущённые» матрицы:

Мх) = fy (x, y (x)), S = да (у (а), у (Ь)), Т = др (у (а), у (Ь)) были достаточно близки к, А (ж), S и Т соответственно.

Пусть Q — один из параметров математической модели и требуется изучить зависимость решения краевой задачи (0.50) от Q. В связи с этим перепишем (0.50) в виде: dy g (y (a), y (b), Q) = О, — Qe[Q0iQ]. Будем считать, что при Q = Qo краевые задачи (0.50) и (0.52) совпадают, так что при Q = Qo краевая задача (0.52), решение которой обозначим через y (x, Q), хорошо обусловлена. Это даёт возможность определить при Q = Qo производную решения по параметру Q — вектор-функцию v (x) = уо (х, Qo), которая является решением хорошо обусловленной линейной краевой задачи = A (x)v + R (x), х£ [а, Ь], ^ ^.

Sv (a) + Tv{b) = ф, где.

А (х) = fy{x, y (x, Qo), Qo), R{x) = fQ{x, y (x, Qo), Qo), S = ga{y{a, Qo), y (b, Qo), Qo), T = gp (y (a, Qo), y (b, Q0), Qo), Ф = -gQ (y{a, Qo), y (b, Qo), Qo)-Далее вводится определение параметризованной краевой задачи с параметром /, i, которая имеет следующий вид [22]: dy = f{x, y, Q), x?[a, b], fx = '(°-54) g (y (a), y (b), Q) = о,.

Здесь, в отличие от краевой задачи (0.52), требуется найти вектор-функцию y (x,/j,) с задаваемым значением к-ой компоненты при х = хс, так что в данном случае параметр Q относится к искомым величинам, т. е. Q = Q (ц). Обусловленность параметризованной краевой задачи (0.54) можно определить из следующей краевой задачи, решением которой являются производные z[x) = уй (х,[Ао) и в = о) решения краевой задачи (0.54) по параметру fi при ц = /л0: dz = A (x)z + R (x)6, х G [a, b],.

JjJU d9 o (0.55) dx '.

Sz (a) + Tz (b) = фв, zk{xc) = 1, xc G [a, b], где A (x), S, T — те же матрицы, a R (x) и ф — те же векторы, что и в краевой задаче (0.53).

Заметим, что при Q = Qo, которому соответствует значение ji = до = yk{xc, Qo), решением краевых задач (0.52) и (0.54) является одна и та же вектор-функция у{х, Qq) = у (х, цо). Можно показать, что при определённых условиях краевая задача (0.54), возможно, будет лучше обусловленной, чем (0.52). Тогда вектор-функцию у{х, Qo) следовало бы искать как вектор-функцию у (х, до)> решая краевую задачу (0.54), а не (0.52).

Для решения нелинейной краевой задачи (0.52), а также и для параметризованной краевой задачи (0.54) в работе используется метод множественной стрельбы [16, 19, 30]. Разобьём отрезок [а, Ь] на га частей: а — хi < х2 < ¦ ¦ ¦ < хт < xm+i — Ь, и на каждом из отрезков [xi, xi+1] рассмотрим задачу Коши dy dx=f{x, y, Q), xe[xi, xi+J, = рй, г = 1,2,., га, У.

0.56).

ОС — ЗС^ гдерМ — вектор начальных данных. Одновременно строится решение V (x, p]< задачи Коши для матричного уравнения:

Г = fy{x, y (x, p[%]), Q) V, xe[xi, xi+1], dx (0.57) ос ос^ где у{х, рЩ — решение задачи Коши (0.56). Очевидно, вектор-функции у (х, р^г г = 1, 2,., т, будут давать в совокупности решение краевой задачи (0.52) при выполнении условий, требующих непрерывности решения в точках разбиения отрезка [а, 6], и краевых условий: фМ = g (pW, ptm+1Q) = 0, фИ=у (ж2,р[1])-р[2]=0,.

0.58) ф[ш+1] = у (Хт+ьр[гп1) р[т+1] = 0.

В результате мы имеем систему нелинейных уравнений относительно компонент векторов., р[то+Ч. Обозначим через р и Ф векторы, составленные из векторов ¦ ¦ ¦ и Ф^Ч, фИ,., ф[т+Ч соответственно, что позволяет представить систему (0.58) в векторном виде: ф (р, q) = о.

0.59).

Пусть р° — составной вектор начальных приближений в методе Ньютона для решения системы (0.59). Тогда для уточняющей поправки Ар° имеем систему линейных алгебраических уравнений.

Затем за начальное приближение берётся вектор р1 = р°—Ар°, к нему ищется поправка Ар1 из системы линейных алгебраических уравнений.

Фр (р1,<�Э)Ар1 = -Ф{рЯ) и так далее. При этом на итерациях используется матрица Якоби ФР (р, Q), определяемая на решениях задач Коши (0.56), которая в данном случае имеет вид.

Фр (р, Q) = S.

У (х2,рЩ Т.

— I.

— I.

— I.

V{xm+bPM) -I.

0.60).

Здесь s = 9a (p[1p[m+1], Q), Т = ^(pW, p[m+1], Q).

После того как компоненты вектора р будут найдены, решение задачи (0.52) может быть восстановлено на всём отрезке [а, Ь] с помощью решений задач Коши (0.56).

Заметим, что «приемлемая жёсткость» задач Коши (0.56) достигается за счёт выбора величины отрезка интегрирования [a^^i+i] и> следовательно, за счёт числа разбиений [а, 6]. При этом обычно следят за выполнением условия.

D • {xi+1 — Xi) и 1, max fy (x, p®, Q) < D. (0.61).

XG[Xi, Xi+ ij.

Кроме того, всякий раз при интегрировании однородной системы (0.57) задаётся в качестве матрицы начальных условий единичная матрица (т. е. ортогональная матрица), что с учётом выбора величины отрезка [х^ a^+i] препятствует «сплющиванию» базисных решений. В частности, при использовании метода множественной стрельбы для решения линейных краевых задач типа (0.51) мы имеем, по существу, вариант метода ортогональных прогонок С. К. Годунова.

Применение метода множественной стрельбы к линейной краевой задаче (0.53) сводится к решению системы линейных алгебраических уравнений относительно производных у (х, рЩ = vq{x, pW).

Sql Ч + Tg^+i] =.

— q[m+1] = -v°(xm+l, pN)} где gM — вектор начальных данных для задачи Коши = A (x)v + R (x), X Е [xi: Xi+l] J v — gM, i — 1, 2,., га.

00 — (JO 2.

Заметим, что у{х, рЩ можно записать в виде где матричная функция V (x, p®) определяется из условий (0.57), a v°(x, p^) — решение задачи Коши dv°, 0 A (x)v +R (x), xE[xi, Xi+1], v dnx ' ' ' ' (0.63).

0 = 0.

J0 — oo^.

Тогда систему уравнений (0.62) можно переписать в векторном виде ф p (P, Q) pq = -*Q (P, Q) — (0.64).

Здесь Фp (p, Q) — матрица (0.60), Фq (p, Q) — составной вектор, составленный из векторов.

Pq — составной вектор, составленный из векторов.

Как мы убедимся, реализация метода множественной стрельбы для параметризованной краевой задачи (0.54) осуществляется практически так же, как и в случае краевой задачи (0.52). В результате аналогичных преобразований параметризованная краевая задача (0.54) сводится к следующей системе нелинейных уравнений ф[1] = Ф[1], Р[т+1], д) = о, фИ = уКрИ)-рМ = 0|.

0.65).

Ф[т+1] = 2/(^+1,РН)-р[т+1]=0, ф[ш+2] = p[f д = 0, которую перепишем в векторном виде ф (р,/*)=О.

Здесь Ф и р — составные векторы, составленные из векторов фМ, фИ,., ф[т+Ч Ф^ и ., pN) PK4) Q соответственно.

Для соответствующей линейной краевой задачи (0.55) получаем систему линейных алгебраических уравнений.

Фр (р> = где I — вектор правых частей системы линейных алгебраических уравнений, все компоненты которого равны нулю, кроме последнего, равного единице. Здесь матрица Якоби Фр системы нелинейных уравнений (0.65) имеет вид.

Фр (&М) = где Фр и Фд те же, что и в системе (0.64).

Фд" .

0. 1—1. 0 0.

Приведённые рассуждения позволяют сформулировать единообразный подход к численному построению решения как краевой задачи (0.52), так и краевой задачи (0.54) методом множественной стрельбы, что приводит к определению соответственно нелинейных уравнений (0.59) или (0.65). В дальнейшем будем считать, что д — одна из компонент составного вектора р, т. е. р — либо одна из компонент составного вектора р, либо // — параметр Q. Введём в рассмотрение вектор искомых величин X, компонентами которого являются компоненты вектора р за исключением /л. В соответствии с этими определениями р и X результатом применения метода множественной стрельбы будет система нелинейных уравнений:

Ф (Х,/л) = 0, (0.66) совпадающая с системой (0.59), если /л = Q, или с системой (0.65) после.

Г?1 Ы исключения в ней если ц —.

Совместное решение серии задач Коши (0.56), (0.57) и (0.63) на итерациях по Ньютону позволяет находить вектор Ф и матрицу [Фр Фд] и, следовательно, матрицу Якоби Фх¦ При этом Фх совпадает с Фр, если р, = Q, т. е. в матрице [Фр Фд] вычеркнут последний столбец, или с матрицей [Фр Фд], в которой вычеркнут столбец с номером k+(j — 1) N, отвечающий компоненте вектора р, если р = В результате проблема изучения зависимости решения краевой задачи (0.52) от параметра Q сводится к изучению поведения решения системы нелинейных уравнений (0.66), которая является результатом точного преобразования нелинейной краевой задачи. После завершения итераций определяется вектор производных Xм из решения системы линейных алгебраических уравнений.

Ф х (Х (^), р) Х, = -Фм№), р), (0.67) где Фм — вычёркиваемый столбец матрицы [Фх Фд]. Тем самым становятся известными производные и Q^.

В главе 2 описывается метод продолжения по параметру решения нелинейной краевой задачи, основанный на методе множественной стрельбы. Пусть при Q = Qq известны как решение краевой задачи (0.52), так и производная решения по параметру Q. Тогда для значения Qi, принадлежащего достаточно малой окрестности Qq, мы можем эффективно задать начальное приближение у°(х, Qi) решения краевой задачи (0.52) при Q = Qi, например, в виде y°{x, Qi) = y{x, Qo) + AQyQ (x, Qo), где AQ = Qi — Qq — шаг по параметру Q. Обычно размеры окрестности заранее не известны. Поэтому в процессе построения решения у (х, Qi) шаг по параметру может уточняться. В частности, если при ограничении на число итераций по методу Ньютона заданная точность приближений не была достигнута, то шаг делится пополам. Это один из приёмов адаптации шага. Получив решение у (х, Qi) и производную решения yQ (x, Q1), мы задаём начальное приближение у°(х, Q2), где Qi принадлежит достаточно малой окрестности Qi и так далее.

Указанный способ численного построения зависимости решения краевой задачи от параметра известен как метод продолжения по параметру или метод гомотопии [15]. Точно так же реализуется продолжение по параметру применительно к параметризованной краевой задаче (0.54), где в формировании начального приближения используются вектор-функции у (х, ц) и уц (х,(х).

Возможности метода продолжения по параметру существенно расширяются, если на каждом шаге, после того как решение вместе с производными решения по параметру получены, применять процедуру выбора параметра для продолжения решения на один шаг, которую мы будем называть параметризацией. Далее, в соответствии с результатом параметризации задаётся начальное приближение либо краевой задачи (0.52), либо параметризованной краевой задачи (0.54) с использованием производных решения по выбранному параметру. Обращение к параметризации позволяет организовать регулярный вычислительный процесс численного исследования зависимости решения краевой задачи (0.52) от параметра Q, включая случай, когда при некоторых значениях Q имеет место ветвление решений типа «поворот». Параметр, выбираемый для продолжения решения на один шаг, мы будем называть текущим параметром. Текущим параметром может оказаться как параметр параметризованной краевой задачи, так и параметр Q. В дальнейшем, и это не вызовет недоразумений, мы будем обозначать текущий параметр через [л.

Предположим, что на рассматриваемом интервале по Q краевая задача (0.52) с решением y (x, Q) может иметь особые точки только типа «поворот». Обозначим через Sq график вектор-функции у (х, Q) при некотором значении Q в (JV+ 1)-мерном евклидовом пространстве (х, у). Тогда при непрерывном изменении параметра Q мы получим в пространстве (х, у) гладкую поверхность Sq, состоящую из пространственных кривых Sq. Как было показано, применение метода множественной стрельбы редуцирует проблему численного построения поверхности Еq к системе нелинейных уравнений (0.59) относительно сеточных значений вектор-функции y (x, Q), т. е. относительно компонент составного вектораp (Q). Таким образом, в силу гладкости Eg система.

0.59) определяет гладкую пространственную кривую Гд в (М + 1)-мерном евклидовом пространстве (p, Q), где М — N (m + 1), являющуюся графиком p (Q). При этом в точках ветвления решений краевой задачи (0.52) типа «поворот» будет ветвление решений того же типа системы (0.59), и, следовательно, в окрестности точек ветвления одному и тому же значению параметра Q будет соответствовать несколько решений системы (0.59).

В свою очередь, гладкость Гд означает, что по теореме о неявных функциях ранг М х (М+ 1)-матрицы [Фр Фд] равен М в некоторой окрестности Гд. Иными словами, учитывая определение вектора X и параметра /л системы (0.66), для любой точки Гд найдётся такая компонента /1 вектора р, равная /io, что в окрестности fiQ т. е. в окрестности /ло система (0.66) будет иметь решение Х (/л) и производную.

ХМ.

Существуют различные варианты продолжения системы нелинейных уравнений типа (0.59), определяющих гладкую пространственную кривую, в которых предусмотрена возможность появления ветвления решений типа «поворот» [23]. В данном случае предлагается вариант продолжения с применением на каждом шаге параметризации для выбора текущего параметра, согласованный с определением параметризованной краевой задачи [24]. Выбранный способ параметризации характеризуется тем, что на итерациях улучшается обусловленность соответствующих линейных краевых задач. Отметим, что параметризованной краевой задаче соответствует система нелинейных уравнений (0.66), если параметр (л не совпадает с Q.

В процессе продолжения решения краевой задачи (0.52) по параметру условие (0.61) может нарушаться. При этом задачи Коши (0.56) могут стать слишком жёсткими. В связи с этим в методе продолжения используется процедура адаптации сетки.

Нелинейные краевые задачи в методе продолжения решаются методом множественной стрельбы, основаном на совместном решении задач Коши (0.56), (0.57) и (0.63). Заметим, что мы можем выбрать тот или иной метод интегрирования задачи Коши (0.56) в зависимости от особенностей конкретной задачи для достижения нужной точности, в то время как для вычисления элементов матрицы Якоби из решения задач Коши (0.57) и (0.63) можем выбрать метод с невысокой точностью, требуя только сходимость метода Ньютона.

Для интегрирования задачи Коши (0.56) в данной работе используется диагонально-неявный метод Рунге-Кутта третьего порядка [10]: у (х о + h) = у (х о) + h (0.5ui + 0.5мг),.

Щ = f (xо + jh, у (хо) + jhuh Qj, (0.68).

U2 = f (x0 + (1 — 7) h, у (х0) + h[(1 — 2j) m + ju2], Q), где 7 = ^e^- Метод (0.68) пригоден для решения жёстких задач и является А-устойчивым. Уравнения в задачах Коши' (0.57) и (0.63) могут быть жёсткими, поэтому метод интегрирования этих задач связан с «замораживанием» матрицы А (х) в средней точке на каждом шаге интегрирования и аппроксимацией вектора правых частей полиномом. Тогда в случае жёстких систем в предлагаемом методе не будет ограничений на длину шага, характерных для методов типа Рунге-Кутта [17].

Выполним интегрирование задачи Коши (0.56) на отрезке [жо, XQ + h]. Тогда в точках х = хо и х = xo + h становятся известными значения векторов у (хо), у{хо + К) и, следовательно, значения вектор-функций: (ж0, у (х о), Qj, f[x о + h, у (х0 + h), Q),.

R (x0) = fQ (xo, y (xo), Qy R (x0 + h) = fq (x0 +h, y (x0 +h), Qy Введём матричную функцию: fy (x, y, Q), где x = xq + 0.5 • h, у = y{x). Для определения у воспользуемся кубической параболой Эрмита, аппроксимирующей y (x, Q) на отрезке [жо, %о 4- h] х = xq •+ h • т, 0 < г < 1, у (х, Q) и 5(ж) = (1 — т) у (х0) + ту[хо + К) + т{ 1 — т) [(1 — т) а° + тЪ° где а0 = у (хо) -y (xQ + h) + hf (x0,y (xQ), Q.

6° = у (х0 + К) — у (х0) — h • f (xо + h, y (xQ + h), Q), Полагая здесь г = 0.5, получим.

У = у{х, Q) = 0.5(2/(ЯГ0) + у{хо + h)) +.

0.125 h f (xQ, у (х0), Q) — f (xо + h, y (x0 + h), Q).

Далее будем решать задачи Коши е «замороженной» матрицей, А = А{х), равной А (х) на каждом отрезке [ж0,ж0 + h], и правой частью R (x), которую будем аппроксимировать полиномом первой степени:

R (xо)(жо + h — х) + R{xо + h)(x — xq).

R (x) h.

R{xо + h) — R (xq) R (xq)(xq + h) — R (x0 4- h) xо ж + ,.

Таким образом, получаем следующие задачи Коши: dV dx V х.

A (x)V, х? [xk, Xk+i]i Xk.

0.69) dv° v dx о.

A (x)v° + R{x), = 0, x G [xk, xk+i],.

0.70) x — x}~ решения которых близки к решениям задач Коши (0.57), (0.63) соответственно.

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

А1.

А2 -I А3.

— I.

Ам+1 bi.

Ь2 Ъ3.

— I.

Ам.

-/ 0 1 ьм 0 X1 ¦ F1 '.

X2 F2.

X3 рз.

Xм FM.

У. G.

0.71).

Здесь в последней строке единица в столбце с номером (г — 1) N + jк = 1,2,., М+1, — матрицы размерности NxNbk, Хк, Fk, к = 1,2,., М, — векторы размерности N. В работе предлагается эффективный метод решения системы уравнений (0.71) связанный с использованием метода ортогональных преобразований и учётом блочно-диагональной структуры матрицы. Отметим, что метод ортогональных преобразований сохраняет обусловлено-ость матрицы.

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

В главе 3 рассматриваются 2 примера применения предлагаемого метода продолжения по параметру.

Первым примером является нелинейная краевая задача, моделирующая стационарные режимы работы каталитического реактора с кипящим слоем [22]: Ш, 0 < ж < 1, Q > О, ах dy2 so si + s2) yi — S2yz — Qf{yi)y4 dx dyz, л dy± (, = -S3Wi)2/4,.

2/2(0) — y3(0) = 0, 2,4(0) = 1, </2(l) = 0.

Здесь функции yi (x), уз (х) и у4(ж) описывают безразмерные распределения температуры твёрдой фазы, газовой фазы и концентрации соответственно, г / I 1 -1 он 4- .ял p{yi).

У1 + S4 s5 + ехр

1 + sefei + «4), где Q, Si, i = 0,1,. ., 6, — параметры модели. Данный пример характеризуется множественностью решений: с увеличением параметра so появляются области изменения параметра Q, где число решений тем больше, чем больше soПо этой причине, если использовать методы, не предусматривающие возможность появления особых точек типа «поворот», то численное исследование примера сопровождается вычислительными трудностями. Этот пример может рассматриваться как тест, определяющий эффективность численного исследования нелинейной краевой задачи.

Второй пример иллюстрирует возможность использования метода продолжения для изучения нелинейных колебаний, описываемых автономной системой. Как известно, в этом случае требуется сформулировать проблему в «стандартном» виде, привлекая условие трансверсальности. В работе рассматривается модель Лоренца [31], которая представлена системой уравнений: dx dy. пxz + Qx-y, dz 1 ху — bz. dt У.

Модель Лоренца была выведена из уравнений Навье-Стокса в задаче о тепловой конвекции, поэтому параметры a, Q и b имеют вполне определённый гидродинамический смысл: а — число Прандтля, Q — число Рэлея, а 6 — характеризует размеры системыфункция х отвечает одной из компонент скорости, а у, z соответствуют членам разложения температуры в ряд Фурье.

Модель Лоренца обладает характерными признаками нелинейной динамической модели: имеет богатый набор решений различных типов поведения, колебательные решения типа предельного цикла, хаотические решения, критическое поведение решения в зависимости от параметров модели [24]. В частности, при сг = 16, 6 = 4 существует более двадцати ветвей периодических решений [29]. В работе рассмотрены различные варианты условия трансверсальности, которое вместе с условиями периодичности сводит задачу определения периодических решений к «стандартной» постановке.

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

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

На защиту выносятся следующие результаты:

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

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

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

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

5. Выполнена программная реализация метода продолжения и рассмотрены характерные примеры нелинейных краевых задач с параметром.

Заключение

.

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

Применение метода множественной стрельбы к исходной и параметризованной краевым задачам редуцирует проблему к системе нелинейных уравнений с параметром модели относительно сеточных значений искомой вектор-функции. При этом сетка задаётся в методе множественной стрельбы, исходя из условия ограничения жёсткости серий задач Коши, определяющих систему нелинейных уравнений вместе с матрицей Якоби. Это позволяет продолжать решение системы по параметру применяя метод Ньютона. Для численного интегрирования серий задач Коши используются эффективные методы, учитывающие их жёсткость. Кроме того, решения систем линейных алгебраических уравнений на итерациях находятся методом ортогональных преобразований, учитывающим структуру матрицы Якоби.

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

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

Результаты докладывались на семинаре кафедры дифференциальных уравнений ММФ НГУ под руководством А. М. Блохина. на семинаре «Методы сплайн-функций» под руководством В. JI. Мирошниченко в ИМ СО РАН (г. Новосибирск). на семинаре по вычислительной математике под руководством В. П. Ильина в НВМиМГ СО РАН (г. Новосибирск). на IV сибирском конгрессе по прикладной и индустриальной математике (ИНПРЙМ-2000). на сибирской конференции «Методы сплайн-функций» (г. Новосибирск). на XXXIX международной научной студенческой конференции «Студент и научно-технический прогресс» (г. Новосибирск). на международной конференции «Современные проблемы прикладной математики и механики: теория, эксперимент и практика» (г. Новосибирск). на XV международной конференции по химическим реакторам (г. Хельсинки, Финляндия).

Список работ по теме диссертации.

1. Когай В. В., Фадеев С. й. Применение продолжения по параметру на основе метода множественной стрельбы для численного исследования нелинейных краевых задач // Сиб. журн. индустр. математики. Новосибирск: Издательство Ин-та математики СО РАН — 2001. -Т. 4. — N. 1. — С. 83−101.

2. Когай В. В., Фадеев С. И. Применение продолжения по параметру на основе метода множественной стрельбы для численного исследования нелинейных краевых задач // Сибирская конференция «Методы сплайн-функций». Новосибирск. 2001. С. 47−48.

3. Когай В. В. Применение продолжения по параметру на основе метода множественной стрельбы для численного исследования нелинейных краевых задач // Материалы XXXIX международной научной студенческой конференции «Студент и научно-технический прогресс». Математика. Новосибирск. 2001. С. 109−112.

4. Фадеев С. И., Когай В. В. Метод множественной стрельбы и продолжение по параметру при численном исследовании нелинейных краевых задач //IV сибирский конгресс по прикладной и индустриальной математике (ЙНПРИМ 2000). Тезисы докладов. Часть II. Новосибирск. 2000. С. 101.

5. Фадеев С. И., Когай В. В. Численное исследование нелинейных задач продолжением по параметру // Вычислительные технологии. Т. 6.

Часть 2. CD-ROM: Международная конференция «Современные проблемы прикладной математики и механики: теория, эксперимент и практика» Новосибирск. Институт Вычислительных Технологий СО РАН. 2001.

6. Fadeev S Л., Kogayi V. V., Korolev V.K. Numerical investigation of nonlinear problems by continuation parameter method. Software automated package BPR-Q for mathematical modeling of the catalytic processes // XV International conference on Chemical Reactors. Helsinki. Finland. 2001. 67−7Qp.