CC BY
6УДК 519.62:517.93:521.1 А. В. Гурьянов
НЕКОТОРЫЕ АСПЕКТЫ ЧИСЛЕННОГО РЕШЕНИЯ ДВУХТОЧЕЧНОЙ КРАЕВОЙ ЗАДАЧИ МЕТОДОМ ОРТОГОНАЛЬНОГО ПЕРЕНОСА ГРАНИЧНЫХ УСЛОВИЙ
Кафедра прикладной математики и компьютерных наук, Ивановский государственный университет, ул. Ермака, д. 39, 153025 Иваново, Россия. E-mail : a.v.gur.2008@mail.ru
Двухточечная краевая задача традиционно применяется в механике сплошных сред при описании тепломассопереноса, колебательных процессов, электрических свойств и других явлений, зависящих от граничных условий. Подобные математические модели имеют большое прикладное значение при изучении физических свойств и структуры мезогенных жидких кристаллов. Рассматриваются некоторые подходы к реализации одного из численных методов решения двухточечной краевой задачи путем решения двух задач Коши. Обосновывается выбор методов решения этих задач и связанных задач линейной алгебры. При решении задач механики сплошных сред описываемый метод имеет некоторые преимущества по сравнению с разностными схемами.
Ключевые слова: двухточечная краевая задача, задача Коши, система обыкновенных дифференциальных уравнений, система линейных алгебраических уравнений, численные методы
A. V. Guryanov
SOME ASPECTS OF THE NUMERICAL SOLUTION OF TWO-POINT BOUNDARY VALUE PROBLEM BY THE METHOD OF ORTHOGONAL TRANSFERRING
OF THE BOUNDARY CONDITIONS
Department of Applied Mathematics and Computer Science, Ivanovo State University Ermak str., 39, 153025 Ivanovo, Russia. E-mail : a.v.gur.2008@mail.ru
The two-point boundary value problem is traditionally used in continuum mechanics for describing heat and mass transfer, oscillatory processes, electrical properties and other phenomena that depend on the boundary conditions. These mathematical models are of great practical importance in the study of the physical properties and structure of the mesogenic liquid crystals. Some approaches to the implementation of one of numerical methods for solving the two-point boundary value problem by the two Cauchy problems are discussed. Methods for solving these problems and the problems of linear algebra are chosen. This method has some advantages over the difference schemes in continuum mechanics problems.
Key words: two-point boundary value problem, Cauchy problem, system of ordinary differential equations, system of linear algebraic equations, numerical methods
© Гурьянов А. В., 2014
Введение
Двухточечная краевая задача традиционно применяется в механике сплошных сред при описании тепломассопереноса, колебательных процессов и других явлений, зависящих от граничных условий. Подобные математические модели имеют большое прикладное значение при изучении свойств и структуры жидких кристаллов. Достаточно широкий класс таких задач может быть описан линейной системой обыкновенных дифференциальных уравнений (ОДУ) вида
y'(x) - A(x) y(x) = fx), a < x < b ; (1)
здесь a2R, b2R, x2R, A : [a, b] ! Cnxn, f: [a, b] ! Cn -заданные непрерывные функции, y : [a, b] ! Cn -искомая функция. Система дополняется разделенными линейными граничными условиями
фа y(a) = Ya, Фь y(b) = Yb, (2)
где Va2Cmxn, Фь2 C(n~m)xn, Ya2 Cm, Yb2 Cn-m - заданные матрицы, rank ф0 = m, rank фь = n - m.
Результаты и обсуждения
Существует целая группа методов решения задачи (1), основанных на так называемом «встречном переносе» граничных условий во внутреннюю точку отрезка [a, b]. Рассмотрим один из них. Для этого распространим граничные условия (2) на весь отрезок [a, b]:
ф, (x) y(x) = Yi(x), 9r(x) y(x) = Yr(x), (3)
где ф,: [a, b] ! Cmxn, фг: [a, b] ! C(n-m)xn, Yi: [a, b] ! Cm, Yr : [a, b] ! Cn-m, rank фг = m, rank фг = n - m, такие, что
Ф, (a) = фа, Yi(a) = Ya, фг (b) = фь, Yr(b) = Yb. (4)
Определим функции ф^), Yi(x), фг(x), Yr(x). Для этого продифференцируем (3) по x:
ф,'^) y(x) + фi(x) y'(x) = Yi' (x),
фг'^) y(x) + фг^) y'(x) = Yr'(x), (5)
x 2 [a, b]. Подставим y'(x) из (1):
ф,'^) y(x) + ф,^) (A(x) y(x) + fx)) = Yi'(x), фг'(x) y(x) + фг(x) (A(x) y(x) + fx)) = Yr'(x), (6)
x 2 [a, b];
или
(фl'(x) + ф,^) A(x)) y(x) = y,' (x) - ф,^) fx), (ф/x + фг^)A(x)) y(x) = y/(x) - фг^) fx), (7) x 2 [a, b].
Исходя из (3) получаем, что для всех x 2 [a, b] имеют место соотношения:
ф1 '(x) + ф,^) A(x) = M(x) ф,^), (8)
фг'(x) + фXx) A(x) = N(x) фXx),
Y,'(x) - фl(x) fx) = M(x) Y,(x), (9)
Yг'(x) - фг^) fx) = N(x) Y^x),
где M: [a, b] ! Cmxm, N : [a, b] ! C{n~m)x<n'm) - некоторые функции.
Умножим (8) справа на ф,*^) и ф,*(x) соответственно. Здесь * - эрмитово сопряжение. Одновременно потребуем, чтобы для всех x Е [a, b] имели место соотношения:
ф,'(x) фl*(x) = 0, фг'^) фг*(x) = 0, (10)
что влечет
ф,(x) ф,*(x) = const,
фXx) фг*(x) = const, (11)
x 2 [a, b].
Тогда (8) будет иметь вид:
ф,^) A(x) ф,*^) = M(x) ф,^) фl*(x),
фг^) A(x) ф*^) = N(x) ф^) ф*^), (12)
x 2 [a, b].
Отсюда найдем функции M и N:
M(x) = ф) A(x) ф*(х) (фг(х) ф*(х))-1,
N(x) = (pr(x) A(x) fr*(x) (fr(x) (pr*(x))-1, (13)
x 2 [a, b].
Подставив в (8) и (9) получаем две задачи Коши для переноса граничных условий из концов отрезка [a, b] в произвольную точку x 2 [a, b]: ^'(x) + ф) A(x) (I- Vl*(x) (ф) Vl*(x))-1Vl(x))= 0, yi'(x) - tyt(x) A(x) yt*(x) (yt(x) фl*(x))~1 yl(x) =
= 9t(x) fx), (14)
фг (a) = фа, yl(a) = ya, x 2 [a, b],
фr'(x) + фr(x) A(x) (I - ф*^) (фr(x) фr*(x))-1 фr(x))= 0, Yr' (x) - фг^) A(x) Фг*^) (фr(x) Фг*(x))-1 Yr(x) =
= Ф^) Ax), (15)
Фг (b) = Фь, Yr(b) = Yb, x 2 [a, b],
где I - единичная матрица соответствующей размерности. Очевидно, что решение y(x) краевой задачи (1) в этой точке можно получить исходя из системы линейных алгебраических уравнений (СЛАУ):
(16)
Без ограничения общности можно потребовать выполнения условий:
фа фа* = I, фЬ фЬ* = I. (17)
Вопросы устойчивости данного метода подробно рассматриваются в [1]. В частности, оказывается, что задачи Коши (14)-(15) численно устойчивы, если исходная краевая задача устойчива относительно малых возмущений функций ее определяющих, включая удовлетворительную обусловленность СЛАУ (16). Кроме того, выполнение условия (10) в процессе интегрирования задач Коши (14)-(15), означает, что при перемещении в бесконечно близкую точку х + ёх базис линейного многообразия в ней будет получаться ортогональным проектированием базиса в исходной точке х. Это соответствует наиболее медленному изменению базиса из всех возможных и одновременно означает, что решение задач Коши (14)-(15) также будет изменяться наиболее медленно среди всех извест-
Yl(x)
y(x) =
фг(x) Yr(x)
ных методов переноса граничных условий. Что дает возможность увеличить шаг в методе численного интегрирования. [1, с. 405] Там же автор не рекомендует в ходе численного интегрирования (14) -(15) в силу свойства (17) отказываться от вычисления (фг(х) фг*(х)У1 и (фг(х) фг*(х))-1.
В целом данный метод можно считать одним из самых эффективных методов встречного переноса граничных условий, поэтому определенный интерес представляет обсуждение некоторых аспектов его практической реализации, включая выбор и возможную адаптацию численных методов решения задач Коши (14)-(15) и СЛАУ (16).
Следует отметить, что системы (14)-(15) состоят из двух различных уравнений первого порядка: матричного - для фг(х) и фг(х) и векторного -для уг(х) и уг(х). Поскольку большинство программных реализаций методов численного решения задачи Коши существует только в векторной форме, необходимо выполнить ряд модификаций классического алгоритма. В частности, крайне привлекательным с точки зрения программной реализации выглядит размещение матрицы ф?(х) и вектора у^(х) в непрерывной области памяти, которую затем можно интерпретировать как вектор
(х) Е С(т+1)п, таким образом «механически» приводя систему (14) к векторному виду. Тогда можно интегрировать оба уравнения согласованно, вычисляя правые части при одних и тех же значениях независимой переменной. Очевидно, что сказанное выше справедливо и для системы (15). Такой подход приводит к существенному снижению трудоемкости вычислений из-за наличия у уравнений системы общих матричных подвыражений ф(х) А(х) и ф*(х) (ф(х) ф*(х))-1. Однако необходимо принять дополнительные меры для обеспечения контроля погрешности вычислений и автоматического выбора шага интегрирования.
Это связано с тем, что при реализации численных методов решения задачи Коши для систем ОДУ в алгоритмах контроля погрешности и выбора шага применяются оценки на основе нормы вектора что будет не совсем корректным применительно к компонентам вектора ^(х), содержащим элементы матрицы ф(х). Более предпочтительным в данном случае следует считать вычисление евклидовой нормы ||у(х)|| е = (у(х) у*(х))1/2
Поскольку вычисление квадратного корня - достаточно затратная процедура, необходимо преобразовать все выражения, связанные с оценкой нормы таким образом, чтобы в них входил ее квадрат, что несколько уменьшит трудоемкость оценок.
При выборе конкретных методов численного решения задач Коши (14)-(15) можно ограничиться явными методами, поскольку эти задачи хорошо обусловлены настолько, насколько хорошо обусловлена исходная краевая задача (1)-(2) [1, 2]. Наиболее популярными методами из этой группы являются методы Рунге - Кутты. Очевидно, что из-за высокой трудоемкости вычисления правой части систем (14)-(15) крайне желательно минимизировать количество таких вычислений одновременно с максимально возможным увеличением шага интегрирования. Одним из способов достижения этой цели может служить повышение порядка метода и применение оценки погрешности по Фель-бергу для автоматического прогноза величины следующего шага интегрирования. Некоторые подходы к построению схем методов численного решения задачи Коши на основе разложения в ряд Тейлора и алгоритмы автоматического выбора порядка метода и управления шагом интегрирования можно найти в [3]. В данном случае для решения систем (14)-(15) целесообразно применить методы Рунге - Кутты - Фельберга седьмого-восьмого порядков и Дормана - Принса восьмого порядка, поскольку при прочих равных условиях они позволяют проводить интегрирование системы с гораздо большей величиной шага, чем популярный метод Рунге - Кутты - Фельберга четвертого-пятого порядков.
Поскольку решение задачи Коши в методах переноса граничных условий является лишь одним из этапов решения основной задачи, следует ужесточить требования, предъявляемые к контролю погрешности применяемых методов. В частности, целесообразно исходить из предположения, что локальные погрешности интегрирования суммируются, что приводит к следующему варианту выражения для прогноза шага интегрирования:
í _
К+ = • К • ч/^, (18)
где К , Нк+1 - текущее и прогнозируемое значения шага, ек - максимально допустимое значение погрешности на текущем шаге, ёк - оценка погрешности на текущем шаге, ^ - порядок метода,
^ - «константа безопасности». Желательно ограничить диапазон однократного изменения шага, например, уменьшать его не более чем в 256 раз, а увеличивать - не более чем в 16 раз. Одно из возможных значений «константы безопасности» -15/16. В случае применения метода четвертого порядка для некоторого снижения трудоемкости классического контроля погрешности по методу Рунге можно рекомендовать при необходимости уменьшать шаг в два раза, а увеличивать - в полтора. Кроме того, следует предусмотреть минимально возможное значение шага исходя из разрядности мантиссы и возможных ограничений на время выполнения программы. Это значение можно определять автоматически, например, взяв некоторую долю от полусуммы модулей границ интервала интегрирования. Ситуацию, когда прогнозируемая величина шага становится меньше минимального значения, следует считать аварийной, поскольку требуемая точность решения вряд ли будет достигнута. Однако интегрирование может быть продолжено с шагом, равным этому минимальному значению.
Заключение
Очевидно, что задачи (14)-(15) должны решаться на одной и той же сетке, в узлах которой необходимо получить решение задачи (1)-(2). Поэтому в качестве начального значения шага можно принять расстояние между первыми узлами сетки и сохранять достигнутую на конце каждого интервала величину шага в качестве начального значения для интегрирования на следующем интервале. Для гарантированного получения решений систем (14) и (15) в одних и тех же узлах сетки независимо ни от переменной величины шага, ни от интегрирования в противоположных направлениях, необходимо предусмотреть проверку на «захват» ближайшего узла сетки на текущем шаге. В этом случае шаг должен корректироваться с целью точного «попадания» в узел сетки. Эту процедуру можно выполнять не только в сторону уменьшения шага, но и в сторону небольшого превышения, например, на 1/16 величины шага.
Существенного снижения трудоемкости вычислений можно достичь, применив для вычисления обратной матрицы (ф(х) ф*(х))-1 метод квадратного корня и метод ЬИ-разложения для решения СЛАУ (16).
Как в том, так и в другом случае для контроля погрешности необходимо выполнять итерационное уточнение решения. Кроме того, выполнение свойств (11) или (17) можно применить в качестве дополнительного средства контроля точности численного решения как задач Коши (14)-(15), так и исходной задачи (1)-(2).
С точки зрения программной реализации необходимо использовать систему программирования, поддерживающую комплексный тип данных удвоенной точности. Данные в ходе решения можно хранить в группе согласованных стеков: для узлов сетки xi, вектора ^l(xi) и решения y(xj) системы (1)-(2).
В целом реализация метода ортогонального переноса граничных условий представляет определенный интерес, поскольку, с точки зрения требований к аппаратным ресурсам, в ходе решения задачи Коши расходуется только время работы процессора, а объем занимаемой оперативной памяти практически не зависит от количества узлов сетки. Это позволяет получать решение задачи (1)-(2) даже на вычислительных установках с ограниченными ресурсами.
Список литературы / References
1. Абрамов А. А. О численной устойчивости одного метода переноса граничных условий // Журн. вычислит. математики и математической физики. 2006. Т. 46, № 3. С. 401-406. [Abramov A. A. O chislennoy
ustoychivosti odnogo metoda perenosa granichnykh us-loviy (Numerical stability of a method for transferring boundary conditions) // Zhurnal vychislitel'noy mate-matiki i matematicheskoy fiziki. 2006. Vol. 46, № 3. P. 401-406. (in Russian)].
2. Абрамов А . А. О переносе граничных условий для систем линейных обыкновенных дифференциальных уравнений (вариант метода прогонки) // Журн. вычислит. математики и математической физики. 1961. Т. 1, № 3. С. 401-406. [Abramov A. A. O pere-nose granichnykh usloviy dlya sistem lineynykh obyknovennykh differentsial'nykh uravneniy (variant metoda progonki) (Transfer of boundary conditions for systems of linear ordinary differential equations, (variant of the screw method)) // Zhurnal vychislitel'noy matematiki i matematicheskoy fiziki 1961. Vol. 1, № 3. P. 401-406. (in Russian)].
3. Бабаджанянц Л. К., Большаков А. И. Реализация метода рядов Тейлора для решения обыкновенных дифференциальных уравнений // Вычислительные методы и программирование. 2012. Т. 13. С. 497510. [Babadzhanyants L. K., Bol'shakov A. I. Realizat-siya metoda ryadov Teylora dlya resheniya obyknoven-nykh differentsial'nykh uravneniy (Implementation of the method of Taylor series for the solution of ordinary differential equations) // Vychislitel'nye metody i pro-grammirovanie. 2012. Vol. 13. P. 497-510. (in Russian)].
Поступила в редакцию 30.06.2014 г.
