CC BY
62
Вычислительные технологии
Том 9, № 1, 2004
О ЧИСЛЕННОМ РЕШЕНИИ СИСТЕМ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ С НЕЛОКАЛЬНЫМИ УСЛОВИЯМИ
К. Р. АйДА-ЗАДЕ Национальная академия наук Азербайджана, Баку e-mail: Kamil_Aydazade@rambler.ru
An approach to numerical solution of systems of linear ordinary differential equations with multipoint unshared conditions is suggested. Formulas and an algorithm of the solution are given. Efficiency of the suggested approach is illustrated numerically for linear and non-linear problems with non-local conditions.
Введение
В работе предложен и исследован численный метод решения систем линейных обыкновенных дифференциальных уравнений с переменными коэффициентами при многоточечных неразделенных условиях, являющихся разновидностью так называемых нелокальных краевых условий. Такие задачи возникают из задач оптимального управления с промежуточными условиями, для решения которых используются необходимые условия оптимальности в форме принципа максимума [1, 2].
Предлагаемый метод сдвига краевых условий в определенном смысле развивает известный метод переноса, исследованный в [3-5]. Так, в частном случае многоточечных неразделенных условий для двухточечных неразделенных краевых условий К. Мошинским предложено [6] сведение их к разделенным краевым условиям за счет удвоения порядка системы уравнений, после чего можно применить классические схемы метода переноса [5]. В случае многоточечных условий можно развить этот подход, но тогда порядок системы существенно возрастет. В [7] исследованы методы переноса для различных классов дискретных краевых задач, но использование этих результатов для непрерывных задач не всегда возможно. Более подробный анализ предлагаемого подхода будет приведен в третьем разделе.
Для данной работы важно предполагаемое условие неавтономности рассматриваемой динамической системы. В случае автономных систем для решения задачи эффективно использовать фундаментальную матрицу решений дифференциальной системы, в результате чего решение задачи приводится к решению системы линейных алгебраических уравнений [2].
Отметим общую проблему линейных систем дифференциальных уравнений с нелокальными краевыми условиями. Проблема связана со сложностью получения конструктивных необходимых, а тем более достаточных условий существования, единственности решения
© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2004.
этих краевых задач. Этим исследованиям посвящено много работ, начиная с работ Тамар-кина, Валле-Пуссена и др. [8, 10], библиография этих работ приведена в [9].
1. Постановка задачи
Рассмотрим следующую линейную неавтономную систему:
¿(t) = A(t)^(t) + B(t), t G [io,ifc], (1.1)
где ¿(t) G En — искомая вектор-функция; матрица A(t) ф const, t G [t0, tk] размерности (n x n), вектор-функция B(t) G En заданы. Пусть n начально-краевых условий системы в общем случае имеют смешанный характер, а именно n0 условий являются многоточечными неразделенными:
k n
£ ) = 7, V =1, 2,..., no, (1.2)
j=0 i=1
ni и n2 — разделенные условия, заданные соответственно на левом и правом концах интервала:
n
J]7vVi(to) = n0, v =1, 2,...,ni; (1.3)
i=1
n
X^t^tk) = 7, v = ni + 1,ni + 2,...,ni + n2 = na. (1.4)
i=1
Здесь n = n0 + ni + n2 = n0 + na; упорядоченные в порядке возрастания моменты времени tj G [t0, tk] и коэффициенты aj,^ 7, v = 1, 2,..., n0 , j = 0,1,..., k, 7^, 70, v = 1, 2,..., ni, 7^, 7k, v = ni + 1, ni + 2,..., na, i = 1, 2,..., n, заданы.
Предположим, что условия (1.2)-(1.4) линейно независимы, а в целом все данные, участвующие в краевой задаче (1.1)—(1.4), таковы, что она имеет единственное решение. Условия, при которых это имеет место, для различных видов краевых задач исследовались во многих работах (например, [6-8]).
Ясно, что условия (1.3), (1.4) — частные случаи условий (1.2), но, учитывая, что условия (1.2) для численного решения краевой задачи представляют наибольшую сложность, будем различать условия вида (1.2) от (1.3), (1.4). Что касается разделенных краевых условий вида (1.3), (1.4), то относительно них хорошо развит аппарат их переноса из одного конца в другой [3, 5]. В связи с этим будем считать, что в зависимости от соотношения количества условий при t = t0 и t = tk уже произведен перенос меньшего числа условий в конец, в котором задано большее число условий. Поэтому вместо условий (1.3), (1.4) зададим na = ni + n2 условий в одном из концов, например в правом, записав их в векторном виде
7*x(tk) = nv, v = 1, 2,..., na, (1.5)
или в матричном виде
7*z(tk) = 7, (1.6)
где вектор Yv = (7vi,..., 7vn)*, v = 1, 2,...,na; * — знак транспонирования; 7 — матрица размерности (na x n), столбцами которой являются векторы 7v, v =1, 2, ...,na.
Неразделенные условия (1.2) также запишем в векторном виде
к
) = Д,, V =1, 2,...,по, (1.7)
з=о
и матричном виде
к
) = Р, (1.8)
3 = 1
где а^ = (а3и1, а3и2---, а3ип )*, V = 1, 2 ,...,п0, а3 — матрица размерности (п0 х п), строками которой служат векторы а3и, V =1, 2,..., п0.
Здесь приводятся векторная и матричная записи условий (1.5)—(1.8), так как обе эти формы будут использованы в следующих разделах.
Решение поставленной задачи заключается в определении вектор-функции х(£), удовлетворяющей линейной системе дифференциальных уравнений (1.1), п3 — начальным условиям разделенного типа (1.5) или (1.6) и п0 — многоточечным неразделенным условиям типа (1.7) или (1.8).
2. Численный метод решения и его обоснование
Предлагаемый ниже подход численного решения поставленной задачи с многоточечными неразделенными условиями (1.1), (1.7) или (1.8) основан на идее классического переноса краевых условий. Он заключается в поэтапном уменьшении используемых точек в нелокальных условиях за счет сдвига участвующих в них всех данных в каком-либо направлении: влево или вправо, в зависимости от специфики задания разделенных условий. В результате осуществления поэтапного последовательного сдвига относительно всех неразделенных многоточечных условий получаются локально заданные краевые условия типа (1.3) и (1.4).
Рассмотрим одно из условий (1.7), в частности V-е:
к
) = Р. (2.1)
3=0
Определение 1. V-е условие (2.1) будем называть сдвинутым вправо относительно х(£)-решения системы (1.1) вектор-функциями а3,(¿), 3 = 0,1,...,к, и функцией (¿), такими, что
а1 (¿о) = ¿3, 3 = 0,1,...,к, &(¿о) = Р, (2.2)
если для любого £ из самого левого подынтервала, в данном случае £ £ [¿0 ,¿1], имеет место
к
а0;(г)х(г) + ^ а1*(г)х(г3) = & (¿). (2.3)
3 = 1
Функции а3,(¿), 3 = 0,1,..., к,в(¿), удовлетворяющие соотношениям (2.2), (2.3) при £ £ [¿0,£1], определены в общем случае не единственным образом, как это имеет место для других методов переноса краевых условий [5].
Из (2.2), (2.3) ясно, что соотношение (2.3) при £ = ¿о совпадает с условиями (2.1). Отсюда следует, что при £ = ¿1 справедливо соотношение
к
а0*(*1) + «!*(*!) + £ аЖМ^) = ^ (¿1). (2.4)
7=2
Введя следующие переобозначения:
О, = (¿1) + (¿1), «V = (¿1), 3 = 2, 3,..., к, Д/ = в^(¿1), из (2.4) получим новое условие
к
)= в*, (2.5)
7=1
эквивалентное (2.1). Легко заметить, что число значений искомой траектории ж(^) в различных точках (моментах времени), участвующих в условии (2.5), на единицу меньше по сравнению с условием (2.1), а именно отсутствует значение в самой левой точке ж(^0). Следовательно, для дальнейшего сдвига нелокальных условий вправо будут использоваться на одну функцию меньше, а именно «(¿),«2(¿), •••,ак(¿).
Повторив подобную процедуру сдвига вправо с условием (2.5) на интервале [¿1,^2] с помощью некоторых функций «(¿), 3 = 1, 2,..., к, и в(¿) можно получить новое условие, эквивалентное (2.5), но не содержащее значения искомой траектории ж(^) при Ь = ¿1. Продолжая к раз сдвиг вправо V-го нелокального условия для конечного интервала [¿к-1, ¿к], получим условие
ак-1^)^) + «к(¿)х(*к) = в^(¿), * € [¿к-1 М (2.6)
которое при ^ = ¿к после переобозначения примет вид
«к ^¿к) = в* • (2.7)
Условие (2.7) эквивалентно V-му условию (2.1), но уже является локально заданным на правом конце интервала, и его можно приписать к (1.5) или (1.6) в качестве (п3 + 1)-го условия. Таким образом, поэтапный сдвиг v-го нелокального условия завершился. При этом соответственно уменьшилось число неразделенных многоточечных краевых условий вида (1.7) или (1.8) с п0 до п0 — 1.
Из определения 1 и описания самой процедуры сдвига несложно понять, что поэтапный сдвиг краевых условий можно осуществлять как для каждого условия в отдельности, так и по всем условиям (1.7) или (1.8) одновременно или по группам условий, используя в этом случае вместо вектор-функций (¿), в^(¿), 3 = 0,1, ...,к, матричные функции а7(¿), 3 = 0,1,..., к, и вектор-функции в(¿).
Замечание 1. При определенной специфике начально-краевых условий (1.2)-(1.4) может быть эффективней сдвиг нелокальных условий влево. Например, в случае п1 > п2 или большой заполненности матриц при меньших значениях индекса 3 (т.е. большего участия в многоточечных неразделенных условиях значений траектории на начальной стадии) следует осуществлять сдвиг условий (1.2) влево.
Остается решить вопрос выбора функций а7 (¿), 3 = 0,1,..., к, в(^), осуществляющих сдвиг v-го условия (2.1), которые, как указывалось выше, не единственны. В частности, в качестве одного из возможных вариантов выбора таких функций можно использовать те, что приведены в следующей теореме.
Теорема 1. Пусть функции аV(г), V = 0,1,...,к, в(г) определяются из решения следующих (п + 2) нелинейных задач Коши на интервале \Ьо,1\]:
а0°(г) = Бо(Ь)а0 - А*(1)а1, а1 (¿о) = 6?°; (2.8)
мм (г) = Бо(г)м(г), м(го) = 1; (2.9)
Д, (г) = Б0(г)/3„ (г) + а°° (г) в (г), д, (го) = ; (2.10)
а{ (г) = м (г)а{, 3 = 1, 2,...,к; (2.11)
где Б0(г) = (0°*Аа° — а°)/(а0°*а0° + 01)и выполнено условие \\а£\\Е« + Щ = 0. Тогда эти функции осуществляют сдвиг V-го нелокального условия (2.1) на интервале [г0,1\], т.е. для них выполняется соотношение (2.3). Причем имеет место
К(г)\\Б„ + 02(г) = к\\Е« + О = соп^, г £ %М, (2.12)
откуда следует устойчивость задач Коши (2.8), (2.10).
Доказательство. Пусть для х(г), являющейся решением задачи (1), (2), имеет место соотношение
к
а1*(г)х(г) + ^о?и*х(з) = Д(г), г £ [гоМ (2.13)
3=1
где а°(г), Д,(г) — пока произвольные дифференцируемые функции, удовлетворяющие лишь условию
а°°(го) = а°°, Д,(го)= Щ. (2.14)
Продифференцируем (2.13) и учтем (1.1) (аргумент г у функций для краткости записи иногда будем опускать):
• о* I о* ♦ _ о
V х \ V х — в V,
[о?; + а°°*А]х + [—^ + а°°*В] = 0. В силу того что х(г) ф 0 и пользуясь произвольностью функций а°(г) и Р(г), потребуе
равенства нулю выражений в квадратных скобках:
= —ао°*А, Д, = а* В. (2.15)
Умножив обе части (2.13) на произвольную функцию М(г), такую, что
м (го) = 1, (2.16)
запишем
к
м (г)ао°*(г)х(г) + ^ м (г)оР*х(Ьз) = м (г)р (г). (2.17)
3=1
Введем обозначения:
д0°(г) = м(г)аКг), 3(г) = м(г)о1, з = 1,2,...,к,
д2(г) = м (г)ри (г), (2.18)
после чего (2.17) примет вид
к
g?*(t)x(t) + X gf(t)x(ij ) = g2(t), (2.19) j=i
причем из (2.14), (2.16) следует
(to) = äj, j = 0,1,..., k, g2(to) = ßv (2.20) Потребуем от M (t) выполнения условия
g?*(t)g?(t) + g2(t) = const, t g [to,ti]. (2.21)
Продифференцировав (2.18) и (2.21) и учитывая (2.15), получим
g?* g? + g?*g? + 2^2^)2 = 0; (2.22)
g?* = Ma?* + Ma??* = Mg?* - Ma??*A = Mg?* - g?*A; (2.23)
<72 = Mß, + M/i = Mg2 + Ma ?*B = —g2 + giB. (2.24)
MM g2 + Ma ?*B = M My2 v M
Подставив (2.23), (2.24) в (2.22), получим
М „ о*о о *л„о, М о *о „ о*л*„о,9/М 2, П
„1 — „1 + „1 — „1 Л „1 I М„2 + „2„1В 1 =
Отсюда, сгруппировав и использовав обозначение
(¿) = („0*а„0 — „2„о в) / („0*„0 + „2),
получим
Щ = (¿), Мф = 5о(¿)М(¿). (2.25)
Подставляя (2.25) в (2.18), (2.24) запишем
¿о = £о„0 — а*„0, „2 = + „1В.
Переобозначив функции „7(¿) через а7 (¿), 3 = 0,1,...к, „2(^) через в(¿), получим доказательство теоремы. □ Замечание 2. Важно отметить следующее. Как указывалось выше, матричные функции, осуществляющие сдвиг краевых условий, определяются не единственным образом. Например, функции, определяемые задачами Коши (2.14), (2.15), формально отвечают условию (2.13), но одна из линейных задач Коши:
¿(¿) = А(г)ж(г) + в(¿), ж(го) = хо,
а°^) = —А*^)«0(¿), «О^О) = «О,
или обе одновременно в зависимости от собственных значений матрицы А неустойчивы. Выполнение же условия (2.12) относительно функций, осуществляющих сдвиг краевых условий, очень важно при проведении практических расчетов, так как оно исключает наличие быстрорастущих решений задач Коши.
Таким образом, сдвиг одного неразделенного краевого условия на весь интервал [¿о, ¿к] слева направо и приведение его к локально заданному краевому условию сводятся к решению задач Коши относительно (п + 2) нелинейных дифференциальных уравнений (2.8)-(2.10). Решение осуществляется поэтапно на каждом последовательном подынтервале, что в целом эквивалентно однократному решению этой системы на всем интервале [¿о, ¿к], причем объем вычислений не зависит от расположения и числа точек ¿7-, 3 = 0,1,..., к. Ясно, что приведение п0 неразделенных многоточечных условий (1.2) к локально заданным краевым условиям вида (1.3) или (1.4) потребует решения задачи Коши в общей сложности относительно (п + 2)п0 дифференциальных уравнений на всем интервале [¿о, ¿к]. Для получения окончательного решения краевой задачи после сдвига всех условий вправо решается система п алгебраических уравнений и определяется значение вектора ^(¿к). Далее используется какой-либо численный метод решения задачи Коши относительно системы (1.1), в данном случае с начальными условиями на правом конце.
При проведении расчетов можно столкнуться со следующими трудностями, которые присущи не только предлагаемому методу сдвига, а вообще всем методам переноса. Условие (2.14) обеспечивает то, что среди функций, осуществляющих перенос, нет сильно-растущих, но часто возникает ситуация, что перенесенные условия в один конец, скажем, ¿ = ¿к, оказываются "почти" линейно-зависимыми, т.е. соответствующая система линейных алгебраических уравнений относительно значения ^(¿к) плохо обусловлена. Для предотвращения этой ситуации, как известно, лучше использовать методы матричного переноса с сохранением нормы матрицы переносимых краевых условий. Предложено поступить таким же образом для получения функций, осуществляющих сдвиг многоточечных неразделенных краевых условий.
Теорема 2. Пусть матричные функции а7(¿), 3 = 0,1,...,к, размерности (п0 х п) и п0-мерная вектор функция в(¿) определяются решением следующих задач Коши на интервале [¿0, ¿1]:
а0^) = ^(¿К^) — а^А^), а0 (¿о) = а0; (2.26)
М(¿) = ^(¿)М(¿), М(¿0) = Е; (2.27)
в (¿) = ^(¿)в (¿) + а0^)В (¿), в (¿о) = в0; (2.28)
а7(¿) = М(¿)а7, 3 = 1, 2,..., к,
где 50(¿) = (а0Аа0* — а0Вв*)(аоао* +вв*)-1, га^а0 = п0, матричные функции М(¿), б"0^) и единичная матрица Е имеют размерность (п0 х п0). Тогда а7(¿),3 = 0,1,к, и в(¿) осуществляют матричный сдвиг краевых условий (1.8) вправо:
к
(¿^(¿) + £ а7(¿^(¿7) = в(¿), ¿ € [¿о, ¿1] , (2.29)
7=1
причем имеет место
К^ЦЕпохп + ||в(¿)||Епо = ||а0||е„0хп + = СОПй1, ¿ € [¿о,¿1]. (2.30)
Если же матричную функцию S0(t) взять в виде
S0 (t) = a0(t)A(t)a*0 (a0 (t)a0* (t))-1, (2.31)
то вместо (2.30) будет иметь место
l|a0(t)\\2 = а0 = const, t £ [t0,t\\.
Из этой теоремы следует, что если матрицы а3, участвующие в нелокальных условиях (1.8), хорошо обусловлены, то это свойство сохраняется при сдвиге краевых условий.
Доказательство теоремы, которое не приводится из-за громоздкости, использует ход доказательства теоремы 1 и матричные преобразования, примененные в [1, 4] c учетом специфики многоточечных условий.
Из (2.26)-(2.28) ясно, что объем вычислений при использовании матричного сдвига определяется решением задачи Коши для системы дифференциальных уравнений порядка (n0n + n2 + n0) в целом на всем отрезке [t0, tk\, что превышает объем вычислений при векторном сдвиге по формулам (2.8)-(2.11) на (n2 — 1) дифференциальных уравнений. Ясно, что при большем числе неразделенных многоточечных условий разница между объемами вычислительной работы при векторной и матричной схемах сдвига становится существенной. Но при этом преимущество матричного сдвига заключается в сохранении линейной независимости в целом всех краевых условий.
В случае нелинейных систем с нелинейными неразделенными краевыми условиями можно использовать метод последовательной линеаризации с применением в дальнейшем предложенного подхода к линейным нелокальным задачам. Ниже приведены результаты численного решения линейной и нелинейной краевых задач.
3. Сравнение метода
Как указывалось во введении, для решения задач с неразделенными двухточечными условиями
оох(го) + ок х(гк ) = в (3.1)
многие авторы рекомендуют использовать метод Мошинского [6], который заключается в следующем. Вводится 2п-мерной вектор:
'( 2М ) = ( $ + , — г) ) '
относительно которого исходная задача с неразделенными условиями становится задачей с разделенными условиями на интервале г £ [го, (го + г к )/2]:
¿1(1) = А(1)Х1(1) + В (г);
(3.3)
¿2(г) = —А(го + гк — ^(г) — в (го + гк — г);
аог1(1о) + ок ¿2(1о) = в; (3.4)
Ло + Ло + гкл , л
¿1( 2 ) — ¿2 (^^)=0. (3.5)
Далее метод переноса локальных краевых условий (3.4) в точку t = (t0 + tk)/2, использующий идею, примененную выше (см. например [5]), потребует решения задачи Коши относительно системы дифференциальных уравнений (2nn0+п0)-го порядка на полуинтервале [t0; (tk +10)/2] и еще решения задачи Коши относительно системы порядка (n0n + n0) для переноса условий из серединной точки t = (t0 + tk)/2 в один из концов интервала [t0,tk], в котором заданы оставшиеся локальные условия типа (1.5), (1.6). Таким образом, объем вычислений при переносе краевых условий (3.1) методом Мошинского можно оценить порядком системы уравнений, относительно которых решается задача Коши:
3
^Ыош = [(2nn0 + т) + (n0n + n0)]/2 = 2 nn0 + n0 (3.6)
(здесь произведено деление на 2, так как задачи решаются на половине интервала).
Как указывалось в разд. 2, объем вычислений предложенным выше подходом векторного сдвига краевых условий не зависит от числа точек, участвующих в нелокальных условиях, и равен
Квреекдтл = (n + 2)n0 . (3.7)
Ясно, что ^ВреДЛ для любых размерностей системы дифференциальных уравнений и числа неразделенных условий меньше, чем ^ош, кроме частного случая, когда n = 2 и n0 = 1, при которых имеет место Умош = КфедЛ.
Если же сравнить объем вычислений методом матричного сдвига
СГ/л = n0n + ni + n0,
то даже для него при n0 < n/2 (что чаще всего и встречается в практических задачах) имеет место
Кредл < ^Мош-
Что касается многоточечных условий, как нам известно, метод Мошинского никто не использовал для таких задач, так как он становится очень громоздким и неэффективным при числе точек более чем два. А самое важное, объем вычислений методом Мошинско-го зависит от числа точек и определяется системой дифференциальных уравнений, как несложно показать, порядка « kn0n (k — число точек, участвующих в краевых условиях), в то время как объем вычислений предложенным методом определяется и в этом случае формулой (3.7).
4. Численные эксперименты
Приведем результаты численных экспериментов для линейной и нелинейной систем.
Пример 1. Рассмотрим систему дифференциальных уравнений четвертого порядка, заданную на интервале [0; 3]:
xi = 2x\ — x2 + x3 — 6t2 + I6t + 3 cos t — sin t — 12,
x 2 = x3 + 2x4 — 3t2 + cos t + 2 sin t — 2, x 3 = xi — 3x2 + x4 + t2 + I2t + 5sin t — cos t — 11, x4 = 2x2 — x3 + t2 — 8t — 3 sin t — 3 cos t + 7
со следующими условиями:
xi(0) - 2x3(0) - xi(2) + 2x3(2) = 10.7551 ,
2x2(1) + 3x4(1) + x2(3) - 2x4(3) = 13.1827, 2x1(3) - 3x4(3) = -2.3245, x2(3) - x3(3) = -4.8389.
Функции
x1(t) = 2t2 - 3t + sin t + 3, x2(t) = t2 + 2t + sint - 1, x3(t) = 3t2 - 4t - 2cost + 2, x4(t) = 3t - sint + cost + 1
— точные решения поставленной задачи. Первые два условия являются многоточечными неразделенными, другие два — разделенные и заданы на правом конце. Очевидно, что правильнее использовать поэтапный сдвиг вправо первых двух условий. В табл. 1, 2 приведены результаты последовательных двух этапов сдвига вправо первого нелокального условия, в табл. 3 — результаты сдвига вправо второго условия. По результатам сдвига нелокальных условий после решения алгебраической системы относительно x1(3), x2 (3), x3(3), x4(3) исходной системы дифференциальных уравнений в обратном порядке от t = 3 до t = 0 получено искомое решение краевой задачи (табл. 4). Расчеты проводились с использованием метода Рунге — Кутты четвертого порядка при числе разбиений интервала, равном 120.
Таблица 1
Первый этап сдвига вправо первого условия
N < ®12 ®13 а14 аи а13 в1
0 1.0000 0.0000 -2.0000 0.0000 -1.0000 2.0000 10.7551
10 0.7822 -1.1358 -1.5404 0.6663 -0.7761 1.5522 11.0885
20 0.4402 -1.5409 -0.5844 1.1533 -0.4697 0.9393 7.5394
30 0.1861 -1.4956 0.1429 1.3522 -0.2880 0.5760 4.4522
40 0.0008 -1.3501 0.6851 1.4587 -0.2016 0.4031 3.0696
50 -0.1749 -1.1449 1.1996 1.5313 -0.1627 0.3255 3.6592
60 -0.3867 -0.7857 1.7820 1.5247 -0.1497 0.2994 7.1135
70 -0.6577 -0.0528 2.4026 1.2422 -0.1510 0.3019 14.8670
80 -0.8679 1.2068 2.5524 0.3060 -0.1443 0.2885 24.7102
Таблица 2
Второй этап сдвига вправо первого условия
N а11 а22 а13 а14 в1
80 -1.0122 1.2068 2.8409 0.3060 24.7102
90 -0.8235 2.3169 1.7586 -1.0482 25.7089
100 -0.4482 2.4434 0.4069 -1.7797 15.0440
110 -0.1462 2.2317 -0.5587 -2.0739 1.1171
120 0.1096 1.9578 -1.3420 -2.2376 -16.3000
Таблица 3
Сдвиг вправо второго условия
N а21 а22 а23 а24 а22 < в2
40 0.0000 2.0000 0.0000 3.0000 1.0000 -2.0000 13.1827
50 -0.0468 1.2050 0.4685 3.4128 1.4930 -2.9860 13.8519
60 -0.2318 0.2467 1.3838 3.5893 1.8502 -3.7003 13.7380
70 -0.5563 -0.1541 2.5230 3.4098 2.0111 -4.0222 17.8679
80 -1.0116 0.4177 3.7386 2.8461 2.1685 -4.3370 31.2198
90 -1.4574 2.1592 4.4000 1.2784 2.2273 -4.4547 51.3945
100 -1.3506 3.9019 3.0815 -1.2022 1.7478 -3.4956 53.9578
110 -0.7864 4.1489 0.8739 -2.7394 1.0798 -2.1596 31.4939
120 -0.2841 3.7739 -0.7773 -3.3483 0.6893 -1.3786 2.9908
Таблица 4
Решения краевой задачи
N Полученное решение Точное решение
ж1 Х2 Х3 х4 ж1 Х2 Х3 х4
0 2.9998 -1.0005 -0.0000 2.0005 3.0000 -1.0000 0.0000 2.0000
10 2.6222 -0.1904 -0.7500 2.4718 2.6224 -0.1901 -0.7503 2.4715
20 2.4793 0.7293 -1.0047 2.8982 2.4794 0.7294 -1.0052 2.8982
30 2.5566 1.7442 -0.7755 3.3000 2.5566 1.7441 -0.7759 3.3001
40 2.8414 2.8416 -0.0803 3.6987 2.8415 2.8415 -0.0806 3.6988
50 3.3240 4.0116 1.0571 4.1162 3.3240 4.0115 1.0569 4.1163
60 3.9975 5.2476 2.6086 4.5731 3.9975 5.2475 2.6085 4.5732
70 4.8590 6.5466 4.5440 5.0877 4.8590 6.5465 4.5440 5.0878
80 5.9093 7.9093 6.8323 5.6745 5.9093 7.9093 6.8323 5.6746
90 7.1531 9.3406 9.4438 6.3437 7.1531 9.3406 9.4438 6.3438
100 8.5984 10.8485 12.3522 7.1004 8.5985 10.8485 12.3523 7.1004
110 10.2566 12.4441 15.5361 7.9440 10.2567 12.4442 15.5361 7.9440
120 12.1410 14.1410 18.9799 8.8688 12.1411 14.1411 18.9800 8.8689
Пример 2. Рассмотрим нелинейную систему дифференциальных уравнений, определенную на интервале [0; 1.2]:
± 1 = ж3 - *3 + 4* - 1, Ж2 = Ж4 - *2 + 2* + 1,
. = ж1 15*6 - 30*4 + 32*2 - 1 Ж3 = - х2 + х2 + 5(*4 - 2*2 + 2) ,
. = Х2 10*5 - 20*3 + *2 + 20* - 3 Ж4 = - ж2 + ж2 + 5(*4 - 2*2 + 2)
со следующими линейными условиями:
Ж1(0) + 2x2(0) + 2x2(0.8) - Ж3(1.2) = 14.448, Ж1 (0.4) - 2x3(0.4) + £4(1.2) = -2.368, Ж2(0.4) + £4(0.4) + £1(1.2) = -1.8, ж2(1.2) + ж3(1.2) = 1.168.
Точным решением задачи являются функции:
хг{г) = 2г2 - 1, Х2(г) = г2 - з, х3(г) = г3 + 1, х4(г) = г2 - 1.
Первые три условия являются многоточечными неразделенными, четвертое — задано на правом конце. Для решения нелинейной задачи использовался метод линеаризации в окрестности текущей функции х0 = х0(г), в результате чего система уравнений, как легко проверить, принимает вид
х1(г) = х3(г) - г3 + 4г - 1,
х2(г) = х4(г) - г2 + 2г + 1, х1 - х022 2х0°х2 - 2х°° 1Ы6 - зог4 + 32г2 - 1
хз() = (х02 + х02)2х1 + х2 + х02)2х2 - (х12 + х02) + 5(г4 - 2г2 + 2) '
. 2х0х0 -х02 + х02 2х0 10г5 - 20г3 + г2 + 20г - з
х4() = (х12 + х02)2х1 + (х12 + х22)2х2 - (х02 + х02) + 5(г4 - 2г2 + 2) • За начальное приближение взяты функции
х°1(г) = 2г + 1, х0(г) = г- з, х33(г) = г + 1, х4(г) = г - 1,
которые не удовлетворяют краевым условиям.
В табл. 5-8 приведены результаты решения линеризованной задачи с многоточечными неразделенными условиями на первой итерации. Полученное решение было снова использовано для линеаризации. В табл. 9-13 приведены результаты применения предложенного в работе метода, полученные на второй итерации линеаризации. Максимальное отклонение по траектории на третьей итерации при числе разбиений, равном 150, и использовании метода Рунге — Кутты четвертого порядка не превышало величину 10-6.
Первая итерация
Таблица 5
Первый этап сдвига вправо первого условия
N а 11 а12 а1з а14 а12 а1з 01
0 1.0000 2.0000 0.0000 0.0000 2.0000 -1.0000 -14.4480
20 0.9849 1.9761 -0.1579 -0.3161 1.9752 -0.9876 -14.0215
40 0.9429 1.9093 -0.3041 -0.6107 1.9075 -0.9538 -13.1625
60 0.8821 1.8111 -0.4308 -0.8696 1.8106 -0.9053 -12.0624
80 0.8123 1.6943 -0.5355 -1.0875 1.6990 -0.8495 -10.9144
100 0.7414 1.5702 -0.6190 -1.2656 1.5835 -0.7918 -9.8766
Таблица 6
Второй этап сдвига вправо первого условия
N а211 а22 а13 а14 а133 01
100 0.7414 3.1537 —0.6190 —1.2656 -0.7918 -9.8766
120 0.6825 2.9587 -0.6938 —1.6671 -0.7461 -8.7554
140 0.6215 2.7464 -0.7486 -1.9977 -0.6967 -8.0012
150 0.5920 2.6393 -0.7697 -2.1381 -0.6720 -7.8030
Таблица 7
Сдвиг вправо второго условия
N а1! а22 а2з а24 а34 в2
50 1.0000 0.0000 -2.0000 0.0000 1.0000 —2.3680
70 0.9333 -0.0297 -2.0318 0.0023 0.9410 —2.1899
90 0.8785 -0.0574 -2.0555 0.0090 0.8868 -2.1451
110 0.8338 -0.0815 -2.0731 0.0194 0.8371 -2.3117
130 0.7971 -0.1009 -2.0865 0.0326 0.7913 -2.7634
150 0.7664 -0.1156 -2.0969 0.0478 0.7491 -3.5550
Таблица 8
Сдвиг вправо третьего условия
N а31 а32 а3з а34 «31 в3
50 0.0000 1.0000 0.0000 1.0000 1.0000 -1.8000
70 0.0151 1.0810 -0.0012 0.9118 1.0850 -1.4763
90 0.0296 1.1671 -0.0051 0.7982 1.1720 -1.0515
110 0.0420 1.2525 -0.0114 0.6553 1.2563 -0.5386
130 0.0516 1.3286 -0.0198 0.4815 1.3303 0.0340
150 0.0578 1.3847 -0.0296 0.2798 1.3847 0.6138
Вторая итерация
Таблица 9
Первый этап сдвига вправо первого условия
N а 11 а02 а!3 а14 а12 а?3 в1
0 1.0000 2.0000 0.0000 0.0000 2.0000 -1.0000 -14.4480
20 0.9868 1.9752 -0.1578 -0.3157 1.9725 -0.9863 -14.0080
40 0.9498 1.9069 -0.3036 -0.6079 1.8960 -0.9480 -13.1093
60 0.8946 1.8098 -0.4289 -0.8605 1.7847 -0.8924 -11.9553
80 0.8265 1.7012 -0.5304 -1.0684 1.6556 -0.8278 -10.7602
100 0.7464 1.5951 -0.6090 -1.2361 1.5236 -0.7618 -9.7052
Таблица 10 Второй этап сдвига вправо первого условия
N а211 а22 а13 а14 а133 в1
100 0.7464 3.1187 —0.6090 —1.2361 -0.7618 -9.7052
120 0.6597 2.9423 -0.6774 —1.6221 -0.7111 -8.6483
140 0.5516 2.7481 -0.7242 -1.9497 -0.6621 -8.0080
150 0.5013 2.6437 -0.7410 -2.0958 -0.6397 -7.8699
Таблица 11
Сдвиг вправо второго условия
N а121 а22 а23 а24 а24 в2
50 1.0000 0.0000 —2.0000 0.0000 1.0000 —2.3680
70 0.9110 0.0160 —2.0420 -0.0013 0.9466 —2.3659
90 0.8205 0.0240 -2.0800 -0.0045 0.9015 -2.4867
110 0.7275 0.0118 -2.1144 -0.0075 0.8639 -2.8004
130 0.6469 -0.0348 -2.1401 -0.0059 0.8305 -3.3629
150 0.6078 -0.1006 -2.1495 0.0051 0.7963 -4.2114
Таблица 12
Сдвиг вправо третьего условия
N ah a32 а3з «34 < Аз
50 0.0000 1.0000 0.0000 1.0000 1.0000 -1.8000
70 -0.0082 1.0746 0.0007 0.9193 1.0927 -1.4596
90 -0.0130 1.1532 0.0026 0.8185 1.1942 -1.0077
110 -0.0099 1.2332 0.0049 0.6922 1.3004 -0.4532
130 0.0033 1.3097 0.0059 0.5335 1.3981 0.1758
150 0.0179 1.3725 0.0044 0.3406 1.4682 0.8147
Таблица 13
Полученные после второй итерации и точное решения задачи
N Полученное решение Точное решение
Х1 Х2 Хз Х4 Xi Х2 Хз Х4
0 -0.9999 -3.0001 0.9997 -0.9997 -1.0000 —3.0000 1.0000 -1.0000
15 -0.9711 -2.9857 1.0015 -0.9853 -0.9712 —2.9856 1.0017 -0.9856
30 -0.8848 -2.9424 1.0136 -0.9422 -0.8848 -2.9424 1.0138 -0.9424
45 -0.7408 -2.8704 1.0464 -0.8702 -0.7408 -2.8704 1.0467 -0.8704
60 -0.5392 -2.7696 1.1104 -0.7693 -0.5392 -2.7696 1.1106 -0.7696
75 -0.2801 -2.6399 1.2158 -0.6397 -0.2800 -2.6400 1.2160 -0.6400
90 0.0367 -2.4815 1.3730 -0.4813 0.0368 -2.4816 1.3733 -0.4816
105 0.4111 -2.2943 1.5924 -0.2941 0.4112 -2.2944 1.5927 -0.2944
120 0.8430 -2.0782 1.8844 -0.0782 0.8432 -2.0784 1.8847 -0.0784
135 1.3326 -1.8334 2.2594 0.1665 1.3328 -1.8336 2.2597 0.1664
150 1.8798 -1.5598 2.7278 0.4400 1.8800 -1.5600 2.7280 0.4400
Заключение
Исследовано численное решение систем обыкновенных дифференциальных уравнений с многоточечными неразделенными краевыми условиями. Предложенный метод решения основан на поэтапном переносе, названном в работе сдвигом нелокальных условий, приводит к решению вспомогательных задач Коши, как и в классических методах переноса. Основное достоинство метода заключается в независимости объема вычислений от количества "точек", участвующих в краевых условиях.
Проведенные численные эксперименты подтвердили устойчивость метода и возможность получения им результатов с высокой точностью и его применимость к нелинейным системам с использованием методов линеаризации.
Автор приносит благодарность аспиранту В.М. Абдуллаеву за разработанные им алгоритмы, программы и проведенные многочисленные эксперименты.
Список литературы
[1] МОИСЕЕВ Н.Н. Численные методы в теории оптимальных систем. М.: Наука, 1971. 424 с.
[2] БраЙСОН А., Хо Ю-Ши. Прикладная теория оптимальных систем. М.: Мир, 1972. 544 с.
[3] Метод "прогонки" для решения разностных уравнений / С.К. Годунов, В.С. Рябенький. Введение в теорию разностных схем. М.: Физматгиз, 1962. С. 283-309.
[4] AБPAMOB A.A. Вариант метода прогонки // Журн. вычисл. математики и мат. физики. 1961. Т. 1, № 2. С. 349-351.
[5] AБPAMOB A.A., Бураго Н.Г. и др. Пакет прикладных программ для решения линейных двухточечных краевых задач // Сообщения по программному обеспечению ЭВМ. М.: ВЦ АН СССР, 1982. 63 с.
[6] MüSZYNSKI K. A method of solving the boundary value problem for a system of linear ordinary differential equations // Algorytmy. 1964. Vol. 11, N 3. P. 25-43.
[7] Самарский A.A., Николаев Е.С. Методы решения сеточных уравнений. М.: Наука, 1978. 590 с.
[8] VALLEE-POUSSIN Ch.J. Sur l'equation differentielle lineare du second order defermination d'une integrale par deux valeurx assignees. Extension aux eqution d'orde n // J. Math. Pura et Appl. 1929. N 9. P. 125-144.
[9] Ешуков Л.Н., Веков A.A., Степанов A.n. Проблемы и библиография теории краевых задач для обыкновенных дифференциальных уравнений // Тр. Рязан. ра-диотехн. ин-та. 1972. Вып. 42. С. 164-192.
[10] Кикурадзе И.Т. Некоторые сингулярные краевые задачи для обыкновенных дифференциальных уравнений. Тбилиси: Изд-во ТГУ, 1975. 352 с.
[11] БеллмАН Р., КАЛАБА Р. Квазилинеаризация и нелинейные краевые задачи. М.: Мир, 1968. 138 с.
Поступила в редакцию 8 июля 2003 г.
