Boundary value problems solving method with the implicit use of the Taylor expansions

Full Text

В работе [1] предложен неклассический метод численного решения краевых задач, основанный на применении рядов Тейлора со старшими производными, порядок которых выше порядка старшей производной заданного дифференциального уравнения. В работе [2] этот метод был применён к решению краевых задач для уравнения Лапласа, а в [3, 4] он развит для решения обыкновенных дифференциальных уравнений, где исследованы вопросы сходимости и устойчивости метода. Данная работа является дальнейшим развитием указанных работ. Пусть в области D с кусочно-гладкой границей Γ задано уравнение L[u(x, y)] ≡ q0 (x, y)u + q1 (x, y)ux + q2 (x, y)uy + q3 (x, y)uxx + + q4 (x, y)uxy + q5 (x, y)uyy = F (x, y), (1) где L — линейный дифференциальный оператор, u(x, y) — искомая неизвестная функция, qi (i = 0, 1, . . . , 5), F — дифференцируемые необходимое число раз заданные функции. Не касаясь всей широты и строгости задачи, рассмотрим предлагаемый метод на примере краевой задачи для уравнения (1) в частных производных, при этом граничные условия пока формулировать не будем. Выполним дискретизацию области по переменным x и y. Обозначим ∆xik = xi − xk (i, k = 1, n) и ∆yjs = yj − ys (j, s = 1, m) соответственно. Рассмотрим шаблон, представленный на рис. 1. 1◦ . Разложим в ряд Тейлора функцию u(x, y) в точке (xi , yj ), ограничившись членами со второй производной: (∆xki )2 uxx (xi , yj )+ 2 (∆ysj )2 uyy (xi , yi ) = u(xk , ys ), (2) + ∆xki ∆ysj uxy (xi , yj ) + 2 u(xi , yj ) + ∆xki ux (xi , yj ) + ∆ysj uy (xi , yj ) + 227 У с о в А. А. Рис. 1 Рис. 2 где пары индексов (k, s) в соответствии с шаблоном принимают следующие значения: (k, s) = {(i, j + 1), (i + 1, j + 1), (i + 1, j), (i, j − 1), (i − 1, j)}. (3) С учётом изменения индексов соотношение (2) представляет систему пяти уравнений. Запишем (1) в точке (xi , yj ): q0 (xi , yj )u(xi , yj ) + q1 (xi , yj )ux (xi , yj ) + q2 (xi , yj )uy (xi , yj )+ + q3 (xi , yj )uxx (xi , yj ) + q4 (xi , yj )uxy (xi , yj ) + q5 (xi , yj )uyy (xi , yj ) = F (xi , yj ). Объединяя пять уравнений (2) и (1) и считая значения функции и её производных в точке (xi , yj ) неизвестными, получим систему шести линейных уравнений (2) относительно шести неизвестных. Обозначим через pksij вектор-строку коэффициентов при неизвестной функции u и её производных в разложении (2) (верхний индекс здесь означает максимальный порядок частных производных, используемых в разложении Тейлора): (2) pksij = 1, ∆xki , ∆ysj , (∆ysj )2 (∆xki )2 , , ∆xki ∆sj , 2 2 (4) где индексы k и s принимают значения (3). Через l обозначим вектор-строку, составленную из функций qm (xi , yj ) (m = = 0, 1, . . . , 5): l = (q0 , q1 , . . . , q5 ). (5) Введём в рассмотрение матрицу M6×6 , строками которой являются пять строк (4), а шестая её строка состоит из элементов вектор-строки (5), а также матрицы X6×1 и B6×1 следующим образом:   (2)  u(x , y )   u(x , y )  pi i+1 ij i j+1 i j   (2) pi+1 j+1 ij  u(xi+1 , yj+1 )  ux (xi , yj )    (2)  u(xi+1 , yj )   u (x , y )     M6×6 =  pi+1 j ij  , X =  y i j  , B =   u(xi , yj−1 )  . uxx (xi , yj )   (2)      pi j−1 ij  u(xi−1 , yj−1 ) uxy (xi , yj )  p i−1 j−1 ij F (xi , yj ) uyy (xi , yj ) l Тогда система уравнений (1), (2) может быть записана в матричном виде M X = B. 228 (6) Способ решения краевых задач с неявным применением рядов Тейлора Если матрица M является невырожденной, т.е. detM = 0, то решение (6) имеет вид X = M −1 B. (7) Обозначим коэффициенты обратной матрицы M −1 через m−1 (α, β = 1, 2, . . . , 6). αβ Тогда первая строка матричного равенства (7) даёт следующее соотношение: u(xi , yj ) − m−1 u(xi , yj+1 ) − m−1 u(xi+1 , yj+1 )− 11 12 −1 − m13 u(xi+1 , yj ) − m−1 u(xi , yj−1 ) − m−1 u(xi−1 , yj−1 ) = m−1 F (xi , yj ). (8) 14 15 16 Обозначим через N общее количество внутренних узлов области D, ограниченной границей Γ, и пусть известно значение функции u(x, y) на этой границе: u(x, y) Γ = ϕ(x, y). (9) Тогда, записывая соотношение (8) в каждом внутреннем узле дискретизации, получим систему N уравнений относительно стольких же неизвестных искомой функции u = u(x, y) в этих узлах. При этом за счёт неравномерной сетки дискретизации по переменным x и y всегда можно добиться того, что граничные узлы точно будут принадлежать границе Γ и значения функции в этих узлах будут определяться равенством (9) (эти члены в соотношениях (8) необходимо перенести в их правую часть). 2◦ . Увеличить точность численного решения краевой задачи (1), (9) (не изменяя количества узлов сетки) можно увеличением количества членов в разложении Тейлора, учитывая частные производные более высокого (чем второй) порядка. При этом, если порядок старшей частной производной равен n, то число узлов шаблона равно 2n + 1. Проиллюстрируем этот вариант на конкретном примере при n = 4. Рассмотрим шаблон из 9 узлов, представленный на рис. 2. Запишем разложение функции u(x, y) в точке (xi , yi ) до частных производных четвёртого порядка: (∆xki )2 uxx (xi , yj )+ 2 (∆xki )3 (∆ysj )2 uyy (xi , yj ) + uxxx (xi , yj )+ + ∆xki ∆ysj uxy (xi , yj ) + 2 3! (∆ysj )4 uyyyy (xi , yj ) = u(xk , ys ), (10) + ...+ 4! u(xi , yj ) + ∆xki ux (xi , yj ) + ∆ysj uy (xi , yj ) + где пары индексов (k, s) в соответствии с шаблоном принимают следующие значения: (k, s) = {(i + 2, j + 1), (i + 1, j − 1), (i + 1, j), (i + 1, j + 1), (i, j − 1), (i, j + 1), (i − 1, j − 1), (i − 1, j), (i − 1, j + 1)}. Таким образом, (10) представляют из себя 9 уравнений. Примем в качестве неизвестных функцию u(xi , yj ) и все её частные производные до четвертого порядка включительно в этой же точке. В итоге число неизвестных равно 15. Присоединим к системе девяти уравнений (10) шесть уравнений, полученных из (1) следующим образом: уравнение (1), а также уравнения, полученные дифференцированием (1) по переменным x, y, xx, xy и yy, записываются в точке (xi , yi ): q0 (xi , yj )u(xi , yi ) + q1 (xi , yj )ux (xi , yj ) + . . . + +q5 (xi , yj )uyy (xi , yj ) = F (xi , yj ), ... ρ0 (xi , yj )u(xi , yj ) + ρ1 (xi , yj )ux (xi , yj ) + . . . + +ρ14 (xi , yj )uyyyy (xi , yj ) = Fyy (xi , yj ), (11) где ρα (α = 0, 1, . . . , 14) — функции, зависящие от qβ (β = 0,1, . . . ,5) и их частных производных. 229 У с о в А. А. Объединим девять уравнений (10) и 6 уравнений (11) в одну систему и будем её рассмтаривать относительно 15 неизвестных u(xi , yj ) и всех её частных производных в точке (xi , yj ) до четвёртого порядка включительно. Тогда эта система также может быть записана в матричном виде (7), где элементами матрицы M15×15 являются коэффициенты при функции u и её производных в системе (10), (11), а матрицы X15×1 B15×1 имеют вид X = u(xi , yj ), ux (xi , yj ), uy (xi , yj ), uxx (xi , yj ), uxy (xi , yj ), uyy (xi , yj ), uxxx (xi , yj ), uxxy (xi , yj ), uxyy (xi , yj ), uyyy (xi , yj ), uxxxx (xi , yj ), uxxxy (xi , yj ), uxxyy (xi , yj ), uxyyy (xi , yj ), uyyyy (xi , yj ) , B = u(xi+2 , yj+1 ), u(xi+1 , yj−1 ), u(xi+1 , yj ), u(xi+1 , yj+1 ), u(xi , yj−1 ), u(xi , yj+1 ), u(xi−1 , yj−1 ), u(xi−1 , yj ), u(xi−1 , yj+1 ), F (xi , yj ), Fx (xi , yj ), Fy (xi , yj ), Fxx (xi , yj ), Fxy (xi , yj ), Fyy (xi , yj ) . Обозначим коэффициенты обратной матрицы M −1 также через m−1 (α, β = 1, 15). αβ Тогда первая строка матричного равенства (7) даёт следующее соотношение: u(xi , yj ) − m−1 u(xi+2 , yj+1 ) − m−1 u(xi+1 , yj−1 ) − . . . 11 12 − m−1 u(xi−1 , yj+1 ) = m−1 F (xi , yj ) − . . . − m−1 Fyy (xi , yj ). (12) 10 1 10 1 15 Записывая (12) для каждой внутренней точки области дискретизации, получим замкнутую систему уравнений для всех узлов сетки, причём учёт граничных условий осуществляется аналогично рассмотренному случаю разложения Тейлора до производных второго порядка. Аналогично можно получить расчётную схему для любого порядка n > 4 старшей производной в разложении Тейлора. 3◦ . Исследуем теперь на конкретных примерах погрешность метода для различных значений порядка n старшей производной в разложении Тейлора. Обозначим через u(n) (xi , yj ) численные решения заданной краевой задачи в узле (xi , yj ), а через ua (xi , yi ) — аналитическое решение этой же задачи в заданном узле сетки. В качестве оценки погрешности выберем следующую норму: δ = i, j u(n) (xi , yj ) − ua (xi , yj ) a 2 i,j (u (xi , yj )) 2 1/2 . (13) Пример 1. Найти решение уравнения ∂2u ∂2u + 2 =0 ∂x2 ∂y в квадрате 0,3 x 0,9, 0,3 y 0,9, если u(0,3; y) = ln(0,09 + y 2 ), u(x; 0,3) = = ln(x2 + 0,09), u(0,9; y) = ln(0,81 + y 2 ), u(x; 0,9) = ln(x2 + 0,81). Аналитическое решение данной задачи имеет вид ua (x, y) = ln(x2 + y 2 ). Для численного решения использовалось разложение до n = 4 с одинаковыми шагами h дискретизации по осям x и y. Расчёты выполнены при трёх значениях шага дискретизации: h = 0,2; h = 0,1; h = 0,05. В результате применения изложенной выше методики получены следующие оценки погрешности (13): при h = 0,2 величина δ h=0,2 = 6,55 · 10−5; при h = 0,1 — δ h=0,1 = 9,1 · 10−7 и при h = 0,05 — δ h=0,05 = = 1,25 · 10−8 . Нетрудно получить следующие результаты: δ h=0,2 / δ h=0,1 ≈ 72, 230 Способ решения краевых задач с неявным применением рядов Тейлора δ h=0,1 / δ h=0,05 ≈ 72,8, откуда следует, что при n = 4 погрешность метода имеет порядок o(h6 ). При использовании классического метода сеток [5] имеем порядок o(h2 ). Пример 2. Выполнить оценку погрешности метода для дифференциального уравнения (1) в области D = {(x, y) : 25 x2 + y 2 169}, если u(x, y) x2 +y2 =25 = = ln 25, u(x, y) x2 +y2 =169 = ln 169. В данном примере было выполнено исследование при двух равных шагах дискретизации по осям Ox и Oy: h = 1 и h = 0,5 и значениях n = 2, n = 4 и n = 6. При n = 2 получены следующие оценки нормы (13): δ h=1 = 1,47 · 10−4 , δ h=0,5 = 3,8 · 10−5 , поэтому δ h=1 / δ h=0,5 ≈ 3,87, откуда имеем погрешность метода близкой к o(h2 ). При n = 4 имеем δ h=1 = 1,148 · 10−5 , δ h=0,05 = 8,27 · 10−7 , δ h=1 / δ h=0,5 ≈ 13,85, откуда погрешность метода в этом случае составляет величину порядка o(h4 ). При n = 6 имеем δ h=1 = 3,6·10−7, δ h=0,5 = 5,28·10−9, δ h=1 / δ h=0,5 = 70,8, откуда погрешность метода o(h6 ). Из данного примера следует, что не изменяя числа узлов (h = const), за счёт увеличения порядка разложения Тейлора можно существенно увеличить точность расчётов (до нескольких порядков), в чём можно убедиться, сравнивая величины погрешности δ во втором примере, например, при h = 0,05 и различных значениях n.

About the authors

Alexandr A Usov

Samara State Technical University

Email: ktmz@samgtu.ru
Engineer, Dept. of Higher Mathematics and Applied Computer Science 244, Molodogvardeyskaya st., Samara, 443100, Russia

References

Supplementary files