Modification of finite differences method with use of Taylor expansions

Full Text

Введение. В настоящее время продолжается разработка и модификация существующих численных методов решения дифференциальных уравнений с точки зрения повышения точности и минимизации количества вычислительных процедур. Последняя проблема остается актуальной, несмотря на улучшение технических возможностей вычислительной техники. В методе сеток увеличение точности достигается за счет уменьшения шага разбиения необходимого интервала. Однако из-за этого происходит увеличение размерности системы алгебраических уравнений, что существенно затрудняет вычисления. Метод тейлоровских разложений, предложенный [1] и являющийся модификацией метода сеток, позволяет, не увеличивая количества узлов сетки, за счет увеличения количества членов в ряде Тейлора в разы повышать точность приближенного решения обыкновенного дифференциального уравнения. Таким образом, порядок системы алгебраических уравнений не увеличивается, что, безусловно, является серьезным преимуществом метода. Ключевым моментом является получение оценок аппроксимации и устойчивости, что является важной задачей любых численных методов. В качестве иллюстрации полученных теоретических выводов в работе рассматриваются конкретные примеры, подтверждающие эффективность метода тейлоровских разложений. 152 Модификация метода сеток с использованием разложений Тейлора Целью данного исследования является модификация метода сеток для обыкновенных дифференциальных уравнений второго порядка с помощью разложений Тейлора с тремя, четырьмя и пятью членами; оценка порядка аппроксимации, проверка устойчивости и сходимости метода тейлоровских разложений; сравнительная с классическим методом сеток оценка количества операций и точности. Метод тейлоровских разложений с тремя членами ряда. Рассмотрим краевую задачу первого рода Ly ≡ y + p(x)y + q(x)y = f (x), a x l0 y = y(a) = γ0 , l1 y = y(b) = γ1 , b; (1) (2) где p(x), q(x), f (x) ∈ C 2 [a, b], γ0 , γ1 — заданные числа. Введём тейлоровский оператор T h y [1]. Выполним дискретизацию отрезка [a, b] точками a = x0 < x1 < x2 < . . . < xn−1 < xn = b с шагом h. Через Yh и Y∗ обозначим линейные пространства сеточных функций, определённых h соответственно на множествах {x1 , x2 , . . . , xn−1 } и {x0 , xn }, в которых введём соответственно нормы y h = max |yi |, 1 i n−1 y ∗ h = max{|y0 |, |yn |}. Пусть осуществляется разложение функции в окрестности точки xi в ряд Тейлора до третьего члена включительно: y (xi )h2 ; 2 y (xi )h2 y(xi+1 ) ∼ y(xi ) + y (xi )h + ; i = 1, 2, . . . , n − 1. = 2 y(xi−1 ) ∼ y(xi ) + y (xi )(−h) + = (3) Из уравнения (1) имеем q(xi )y(xi ) + p(xi )y (xi ) + y (xi ) = f (xi ), i = 1, 2, . . . , n − 1. (4) Запишем уравнение (1) в каждой точке xi и возьмём за неизвестные y(xi ), y (xi ), y (xi ). Обозначим для краткости значения функции y(x) в узлах xi через yi и представим систему (3), (4) в матричной форме M X = B, где X = (yi , yi , yi ) , B = (yi−1 , yi+1 , fi ) ,   1 −h h2 /2 M =  1 h h2 /2 . q i pi 1 В предположении, что det M = 0, разрешим систему (3), (4): yi =m−1 (xi )yi−1 + m−1 (xi )yi+1 + m−1 (xi )fi , 11 12 13 yi =m−1 (xi )yi−1 + m−1 (xi )yi+1 + m−1 (xi )fi , 21 22 23 yi =m−1 (xi )yi−1 31 + m−1 (xi )yi+1 32 + m−1 (xi )fi , 33 (5) i = 1,2, . . . , n − 1, 153 П а в л о в а Г. А., Б е л я е в а И. В. где mkl (xi ) — элементы матрицы M −1 . С учётом краевых условий (2) перепишем первое уравнение (5) в виде m−1 (xi )yi−1 − yi + m−1 (xi )yi+1 = −m−1 (xi )fi , 11 12 13 i = 1, 2, . . . , n − 1. В результате получим систему n − 1 линейных уравнений с n − 1 переменными, которую можно разрешить различными методами (например, методом прогонки), после чего можно найти величины yi , yi из второго и третьего уравнений (5) чисто алгебраически. Аппроксимируем дифференциальный оператор Ly на сетке тейлоровским оператором T h y. Для этого составим оператор T h y: (T h y)i = yi + pi yi + qi yi , где yi , yi , yi определяются по формулам (5). Тогда в развернутой форме оператор T h y можно записать так: (T h y)i = m−1 (xi )yi−1 + m−1 (xi )yi+1 + m−1 (xi )fi + 31 32 33 + p(xi )(m−1 (xi )yi−1 + m−1 (xi )yi+1 + m−1 (xi )fi )+ 21 22 23 + q(xi )(m−1 (xi )yi−1 + m−1 (xi )yi+1 + m−1 (xi )fi ). (6) 11 12 13 Составим разностную краевую задачу (точнее говоря, семейство разностных краевых задач, зависящее от параметра h): T h y = f, th y = g, где T h y — оператор тейлоровских разложений, определённый по формуле (6), th y — оператор, имеющий вид (th y)i = yi , i = 0, n, f ∈ Yh — сеточная функция, порождаемая правой частью уравнения (1), g ∈ Y∗ — заданная сеточная функция со значениями g0 = γ0 , gn = γ1 [2]. h Утверждение 1. Пусть y ∈ C 4 [a, b] — произвольная функция, которая, в частности, может быть решением краевой задачи (1), (2). Тогда выполняется неравенство Ly − T h y h cu h2 , где cu — не зависящая от h постоянная. Д о к а з а т е л ь с т в о. Оценим Ly − T h y Ly − T h y h = y (xi ) + p(xi )y (xi ) + q(xi )y(xi ) − yi − pi yi − qi yi y (xi ) − yi 154 h: h+ p c y (xi ) − yi h+ q c h y (xi ) − yi h. (7) Модификация метода сеток с использованием разложений Тейлора По формулам (5) найдём yi , yi , yi . Для этого вычислим матрицу M −1 . При условии, что det M = 2h − h3 q = 0,   − h2 1 − hpi /2 1 + hpi /2  2 − q h2 2 − qi h2 2 − qi h2  i   −1 −1 −1 . M =  −2h 2h 0     pi − q i h pi + q i h 2 − 2h − h3 qi 2h − h3 qi 2 − qi h2 Отсюда получаем, что 1 − pi yi = h h yi−1 + 1 + pi yi+1 − h2 fi 2 2 = 2 − h2 qi h pi (yi+1 − yi−1 ) + yi+1 + yi−1 − h2 fi = = 2 2 − h2 qi yi+1 − yi−1 yi+1 − 2yi + yi−1 2yi pi + + 2 − fi 2h h2 h = h2 . 2 − h2 qi Аналогично методу сеток [3] используем формулы численного дифференцирования y(x1 ) − y(x−1 ) h2 − y (ξ), y (x0 ) = 2h 6 (8) y(x1 ) − 2y(x0 ) + y(x−1 ) h4 (4) y (x0 ) = − y (ξ), h2 12 где ξ ∈ [x−1 , x1 ]. Тогда yi = h2 h2 h2 y (xi ) + y (4) (ξ) + pi y (xi ) + y (ξ) + qi y(xi )− 2 − h2 qi 12 6 2y(xi ) 1 −qi y(xi ) − fi + , = y(xi ) + M h4 2 h 2 − h2 qi откуда y(xi ) − yi h = O(h4 ). Аналогично проведём рассуждения для yi : yi = yi+1 − yi−1 . 2h Полученное выражение в точности соответствует первой формуле из (8). Следовательно, в данном случае оценка Ly − T h y h для метода тейлоровских разложений эквивалентна оценке для метода сеток. Рассмотрим теперь 155 П а в л о в а Г. А., Б е л я е в а И. В. yi = (pi − hqi )yi−1 − (pi + hqi )yi+1 + 2hfi = 2h − h3 qi h3 (yi+1 − 2yi + yi−1 ) (yi+1 − yi−1 ) − qi −pi − qi yi + fi 2h 2h h2 = . h2 1 − qi 2 Подставляя формулы (8) в это выражение, получим h2 h2 y (ξ) − qi h2 y (xi ) − qi y (4) (ξ) − qi yi + fi 6 12 = yi = h2 1 − qi 2 h2 h2 −pi y (ξ) − qi h2 y (xi ) − qi y (4) (ξ) h2 6 12 = y (xi ) + = y (xi ) + M , h2 h2 1 − qi 1 − qi 2 2 откуда y (xi ) − yi h = O(h2 ). −pi y (xi ) − pi Таким образом, оператор T h y, как и разностный оператор, аппроксимирует дифференциальный оператор Ly со вторым порядком точности относительно шага h в сеточной норме · h . Устойчивость разностной схемы. Разностная схема представляет собой систему линейных алгебраических уравнений, которую всегда можно записать в матричной форме [4] Ay = Y, (9) где y — искомый вектор, Y — заданный вектор, определяемый правыми частями разностных уравнений и дополнительными (начальными) условиями. Уравнение (9) можно рассматривать так же, как операторное уравнение, где A — линейный оператор, действующий в пространстве Y, y — элемент этого пространства, Y ∈ Y — заданный элемент. Для разностных схем каждая схема определяет целое семейство уравнений (10) Ah yh = Yh , зависящих от шага сетки h. При каждом допустимом значении h оператор Ah действует в конечномерном пространстве Yh . Предположим, что в пространстве Yh заданы нормы · (1h) и · (2h) , в которых измеряются соответственно решение уравнения (10) и его правая часть. Определение [4]. Уравнение (10) называется корректным, если 1) решение уравнения (10) существует и единственно при любых Yh ∈ Yh ; 2) существует константа M > 0, не зависящая от h, такая, что для Yh ∈ Yh y 156 (1h) M Yh (2h) . (11) Модификация метода сеток с использованием разложений Тейлора При этом условие 1 эквивалентно существованию оператора A−1 , обратh ного к Ah , а условие 2 — равномерной по h ограниченности A−1 . Свойство, выh раженное оценкой (10), называется устойчивостью разностной схемы и означает, что при небольшом изменении входных данных решение изменяется мало. Теперь предположим, что Yh — вещественное конечномерное пространство, в котором введены скалярное произведение (y, ν)h и норма y h = = (y, y)h . Тогда справедливо следующее утверждение. Утверждение 2. Если существует постоянная δ > 0, не зависящая от h, такая, что ∀vh ∈ Yh (Ah vh , vh )h δ vh 2 h, то уравнение (10) корректно и для его решения выполняется оценка yh δ−1 Yh h h. Проверим выполнение этого утверждения в случае разностной схемы, порожденной тейлоровским оператором. Итак, разностная схема в рассматриваемом случае имеет вид m−1 (xi )yi−1 − yi + m−1 (xi )yi+1 = m−1 (xi )fi , 11 12 13 i = 1, 2, . . . , n − 1. (12) При y0 = a, yn = b систему (12) можно представить в следующем виде:  −1 −1 −1 y1 − m12 (x1 )y2 = −m13 (x1 )f1 + m11 (x1 )a;   −m−1 (x )y + y − m−1 (x )y = −m−1 (x )f ;  2 1 2 2 3 2 2  11 12 13  −m−1 (x )y + y − m−1 (x )y = −m−1 (x )f ; 3 2 3 3 4 3 3 11 12 13  . . .  −m−1 (x −1 −1   n−2 )yn−3 + yn−2 − m12 (xn−2 )yn−1 = −m13 (xn−2 )fn−2 ; 11   −m−1 (xn−1 )yn−2 + yn−1 = −m−1 (xn−1 )fn−1 + m−1 (xn−1 )b. 11 13 12 Тогда разностную схему (12) можно записать в операторной форме: T h y = f. Матрица оператора T h является трёхдиагональной:  . . . 0 0   .   −1 −1 . −m11 (x2 )  1 −m12 (x2 ) . 0 0   .   −1 . 0 −m11 (x3 ) 1 . 0 0   Th =  . .   . ... ... ... . ... ...     . −1  . 0 0 0 . 1 −m12 (xn−2 )   . . −m−1 (x . 1 0 0 0 n−1 ) 11  1 −m−1 (x1 ) 12 0 157 П а в л о в а Г. А., Б е л я е в а И. В. В пространстве Yh введём скалярное произведение n−1 y i vi h (y, v)h = i=1 и норму, порожденную этим скалярным произведением: 1/2 n−1 y h 2 yi = . i=1 Тогда ∀vh ∈ Yh 2 2 (Th vh , vh ) = h v1 − m−1 (x1 )v1 v2 − m−1 (x2 )v1 v2 + v1 − m−1 (x2 )v2 v3 − 12 11 12 2 2 − m−1 (x3 )v2 v3 + +v3 − m−1 (x3 )v3 v4 + . . . − m−1 (xn−1 )vn−2 vn−1 + vn−1 = 11 12 11 n−1 2 vi h − h v1 v2 m−1 (x1 ) + m−1 (x2 ) + v2 v3 m−1 (x2 ) + m−1 (x3 ) + . . . + 12 11 12 11 = i=1 +vn−1 vn−2 m−1 (xn−1 ) + m−1 (xn ) 12 11 . Воспользуемся очевидным неравенством −ab − a2 + b2 2 и получим оценку 2 2 v1 + v2 + m−1 (x2 )+ 12 2 2 2 v 2 + vn−2 v 2 + v3 + . . . + m−1 (xn−1 ) + m−1 (xn ) n−1 +m−1 (x3 ) 2 12 11 11 2 2 h 2 2 2 2 max m−1 (xi ) + m−1 (xi+1 ) v1 + 2v2 + 2v3 + . . . + 2vn−2 + vh 2 − h 12 11 2 i=1,...,n−1 2 +vn−1 vh 2 − max m−1 (xi ) + m−1 (xi+1 ) · vh 2 = vh 2 (1 − Mh ), h h h 12 11 (Th vh , vh )h vh 2 h m−1 (x1 ) + m−1 (x2 ) 12 11 −h i=1,...,n−1 где Mh = max i=1,...,n−1 |m−1 (xi ) + m−1 (xi+1 ) . Заметим, что 12 11 m−1 (xi ) = 12 1 + h pi 2 2 2(1 − qi h ) 2 , m−1 (xi+1 ) = 11 Если предположить, что p(x) ≡ 0, а q(x) m−1 (xi ) + m−1 (xi+1 ) = 12 11 158 2 2(1 − qi+1 h ) 2 . 0 ∀x ∈ [a, b], то 1 2(1 − 1 − h pi+1 2 2 qi h ) 2 + 1 2 2(1 − qi+1 h ) 2 1 Модификация метода сеток с использованием разложений Тейлора и 1 − Mh > 0 не зависит от шага h. Таким образом, при сформулированных выше условиях оператор T h y является корректным и схема, им порожденная, является устойчивой. Сходимость решения разностной схемы. Поскольку поставленная задача удовлетворяет основной теореме теории разностных схем [2, стр. 209], то решение, полученное методом тейлоровских разложений, сходится к решению исходной задачи со вторым порядком точности. Метод тейлоровских разложений с четырьмя членами ряда. В данном случае для решения краевой задачи (1), (2) использовались такие разложения: y (xi )h2 y (xi )h3 y(xi−1 ) ∼ y(xi ) + y (xi )(−h) + − ; = 2! 3! y (xi )h2 y (xi )h3 y(xi+1 ) ∼ y(xi ) + y (xi )h + + . = 2 3! После аналогичных предыдущему случаю преобразований получим следующую систему уравнений: m−1 (xi )yi−1 −yi +m−1 (xi )yi+1 = − m−1 (xi )fi + m−1 (xi )fi , 13 14 11 12 i = 1, . . . , n−1. Результат произведённых аналогично предыдущему случаю оценок порядка аппроксимации выглядит следующим образом: h4 2 −pi qi + pi qi + qi + O(h5 ) + c0 h4 ⇒ 12 yi − y(xi ) h M h4 ⇒ yi − y(xi ) h = O(h4 ); yi = y(xi ) 1 − h2 (−p2 + qi + pi ) + O(h3 ) + c1 h2 ⇒ i 6 yi − y (xi ) h M h2 ⇒ yi − y (xi ) h = O(h2 ); yi = y (xi ) 1 − yi = y (xi ) 1 + yi − y (xi ) h h2 2 p − qi − pi + O(h3 ) + c2 h2 ⇒ 6 i M h2 ⇒ yi − y (xi ) h = O(h2 ), где константы c0 , c1 , c2 находятся в процессе преобразований, зависят от qi , pi , qi , pi и не приводятся из-за громоздкости. Теперь можно сделать вывод, что тейлоровский оператор T h y в случае четырёх членов разложения ряда аппроксимирует дифференциальный оператор Ly со вторым порядком точности относительно шага h в сеточной норме · h . Вопросы устойчивости и сходимости данного метода рассматриваются аналогично случаю с тремя членами разложения. Метод тейлоровских разложений с пятью членами ряда. Аналогично предыдущим случаям рассматривается разложение искомой функции в ряд Тейлора с учётом пяти членов разложения: y (xi )h2 y (xi )h3 y (4) (xi )h4 y(xi−1 ) ∼ y(xi ) + y (xi )(−h) + − + ; = 2! 3! 4! 159 П а в л о в а Г. А., Б е л я е в а И. В. y (xi )h2 y (xi )h3 y (4) (xi )h4 y(xi+1 ) ∼ y(xi ) + y (xi )h + + + . = 2 3! 4! Найдена следующая система алгебраических уравнений: m−1 (xi )yi−1 − yi + m−1 (xi )yi+1 = − m−1 (xi )fi + m−1 (xi )fi + m−1 (xi )fi , 11 12 13 14 15 i = 1, . . . , n − 1. Для оценки аппроксимации использовалась формула (7) и были получены следующие оценки: yi = y(xi ) − yi − y(xi ) h h4 h4 y(xi ) + O(h5 ) + c0 ⇒ 2 2 M h4 ⇒ yi − y(xi ) h = O(h4 ); h2 qi + c1 h2 + O(h3 ) ⇒ 6 M h2 ⇒ y (xi ) − yi h = O(h2 ); yi = y (xi ) − yi − y (xi ) h yi = y (xi ) + y (xi ) yi − y (xi ) h h2 qi + pi − p2 + c2 h2 + O(h3 ) ⇒ i 6 M h2 ⇒ yi − y(xi ) h = O(h4 ), где константы c0 , c1 , c2 зависят от qi , pi , qi , pi , qi , pi и здесь не приводятся из-за громоздкости. Исследование метода показало, что результаты, полученные для оценок аппроксимации, сходимости и устойчивости, полностью аналогичны методу тейлоровских разложений с четырьмя членами ряда. Сравнительный анализ числа арифметических операций. Метод сеток и метод тейлоровских разложений, используемые для решения задачи (1), (2), имеют разное количество арифметических операций. Увеличение количество членов ряда в разложении ведёт к резкому увеличению числа арифметических операций в методе тейлоровских разложений. Кроме этого, метод тейлоровских разложений более высокого порядка предполагает вычисление значений производных на каждом шаге, что, в свою очередь, также увеличивает число арифметических операций. Для примера рассмотрим такую задачу: y − 2 2 (x − 1)2 + 1 y + 2y = , x x x2 y(2) = 8,545; y(10) = 287,31. При шаге сетки h = 0,8 на решение этой задачи методом сеток понадобится 360 действий, методом тейлоровских разложений с тремя, четырьмя и пятью членами ряда соответственно — 528, 3735 и 27087 действий. Если метод тейлоровских разложений с тремя членами ряда ещё сопоставим по числу операций с методом сеток, то в случае четырёх и пяти членов количество действий существенно увеличивается. Однако метод тейлоровских разложений с четырьмя и пятью членами ряда даёт более высокую точность, что немаловажно. 160 Модификация метода сеток с использованием разложений Тейлора Кроме того, данное уравнение решалось методом сеток с разным шагом и полученная погрешность сравнивалась с точностью, даваемой методом тейлоровских разложений. Перед нами стояла задача определить, каким должен быть шаг в методе сеток, чтобы оба метода давали сравнимую погрешность. Когда погрешность уравнивалась, проводился подсчет количества арифметических операций. Результаты для сравнения метода сеток и метода тейлоровских разложений с четырьмя и пятью членами ряда приведены в таблице. Здесь требовалось добиться приблизительно одинаковой погрешности вычислений при решении методом тейлоровских разложений и методом сеток и в каждом случае сравнить число вычислительных действий. Сравнение числа операций для метода сеток и метода тейлоровских разложений Метод тейлоровских разложений Шаг ПогрешЧисло ность операций Шаг Метод сеток ПогрешЧисло ность операций 4 члена ряда 0,8 0,4 0,2 1,288 · 10 3,104 · 10−5 9,120 · 10−6 0,8 0,4 0,2 2,251 · 10 2,090 · 10−6 2,060 · 10−6 −4 3735 7875 16155 0,4 0,1 0,05 2,084 · 10−4 1,400 · 10−5 5,037 · 10−6 819 3219 6419 2,203 · 10−6 2,160 · 10−6 2,089 · 10−6 32019 40019 160019 5 членов ряда −6 27087 57137 117237 0,01 0,008 0,002 Полученные результаты ясно показывают, что уже при пяти членах разложения метод сеток теряет свое преимущество относительно небольшого количества операций, поскольку для получения заданной погрешности необходимо брать большое число узлов разбиения. Выводы по работе. Проведенное исследование модификации метода сеток с использованием разложений Тейлора показало следующее: 1) решение, полученное методом тейлоровских разложений с тремя, четырьмя и пятью членами ряда, аппроксимирует точное решение со вторым порядком точности; 2) устойчивое решение, полученное методом тейлоровских разложений, сходится к точному решению со вторым порядком точности; 3) погрешность, даваемая методом тейлоровских разложений с тремя членами ряда, сопоставима с погрешностью метода сеток; однако при увеличении количества членов в разложении происходит резкое улучшение точности метода тейлоровских разложений по сравнению с методом сеток; 4) при применении метода сеток тейлоровских разложений с пятью и большим количеством членов ряда используется меньшее количество арифметических операций по сравнению с методом сеток для достижения высокой точности одного порядка.

About the authors

Galina A Pavlova

Samara State Technical University

Email: pavlova-ga@mail.ru
(Ph. D. (Phys. & Math.)), Associate Professor, Dept. of Applied Mathematics and Computer Science 244, Molodogvardeyskaya st., Samara, 443100, Russia

Irina V Belyaeva

Samara State Technical University

Email: irr.bel@gmail.com
student, Dept. of Applied Mathematics and Computer Science 244, Molodogvardeyskaya st., Samara, 443100, Russia

References

Supplementary files