CC BY
16
Izvestiya of Saratov University. Mathematics. Mechanics. Informatics, 2023, vol. 23, iss. 1, pp. 36-47
Известия Саратовского университета. Новая серия. Серия: Математика. Механика. Информатика. 2023. Т. 23, вып. 1. С. 36-47
mmi.sgu.ru https://doi.org/10.18500/1816-9791-2023-23-1-36-47, EDN: BFDVVG
Научная статья УДК 517.98
Новый подход к формированию систем линейных алгебраических уравнений для решения обыкновенных дифференциальных уравнений методом коллокаций
Л. А. Севастьянов1'2, К. П. Ловецкий1, Д. С. Кулябов1'20
Российский университет дружбы народов (РУДН), Россия, 117198, г. Москва, ул. Миклухо-Маклая,
2Объединенный институт ядерных исследований, Россия, Московская область, 141980, г. Дубна, ул. Жолио-Кюри, д. 6
Севастьянов Леонид Антонович, доктор физико-математических наук, 1профессор кафедры прикладной информатики и теории вероятностей, 2сотрудник лаборатории теоретической физики, 8еуа8Шпоу-la@pfur.ru, https://orcid.org/0000-0002-1856-4643, АиЙогГО: 6741
Ловецкий Константин Петрович, кандидат физико-математических наук, доцент кафедры прикладной информатики и теории вероятностей, lovetskiy-kp@rudn.ru, https://orcid.org/0000-0002-3645-1060, Аи^ГО: 584101
Кулябов Дмитрий Сергеевич, доктор физико-математических наук, 1профессор кафедры прикладной информатики и теории вероятностей, 2сотрудник лаборатории информационных технологий, kulyabov-ds@rudn.ru, https://orcid.org/0000-0002-0877-7063, AuthorID: 360953
Аннотация. Реализован новый алгоритм численного решения одномерных задач Коши и уравнений Пуассона, основанный на методе коллокации и представлении решения в виде разложения по полиномам Чебышева. Предлагается вместо обычного подхода, заключающегося в слиянии всех известных условий — дифференциальных (само уравнение) и начальных/ граничных — в одну систему приближенных линейных алгебраических уравнений (СЛАУ), перейти к методике решения задачи в несколько отдельных этапов. Вначале выделяются спектральные коэффициенты, определяющие «общее» решение исходной задачи. По методу коллокации определяются интерполяционные коэффициенты производной решения, а тем самым и коэффициенты разложения самого решения (кроме начальных). На этом этапе выбор удачного базиса, обладающего дискретной ортогональностью, дает возможность применения весьма эффективных алгоритмов поиска искомых коэффициентов. Трудоемкость приведения матрицы СЛАУ к диагональной форме становится эквивалентной сложности умножения чебышевской матрицы коэффициентов на вектор правой части системы. Затем коэффициенты разложения самого решения (кроме первых одного-двух) получаются с помощью умножения известной трехдиагональной матрицы интегрирования (обратной по отношению к матрице дифференцирования Чебышева) на вектор интерполяционных коэффициентов производной. На последнем этапе учет начальных/граничных условий выделяет «частное» искомое решение, однозначно доопределяя недостающие коэффициенты искомого разложения. Ключевые слова: начально-краевые задачи, метод коллокации, многочлены Чебышева, множества Гаусса - Лобатто, численная устойчивость, дискретная ортогональность
д. 6
Благодарности: Работа выполнена при финансовой поддержке Программы стратегического академического лидерства РУДН.
Для цитирования: Севастьянов Л. А., Ловецкий К. П., Кулябов Д. С. Новый подход к формированию систем линейных алгебраических уравнений для решения обыкновенных дифференциальных уравнений методом коллокаций // Известия Саратовского университета. Новая серия. Серия: Математика. Механика. Информатика. 2023. Т. 23, вып. 1. С. 36-47. https://doi.org/10.18500/1816-9791-2023-23-1-36-47, EDN: BFDVVG
Статья опубликована на условиях лицензии Creative Commons Attribution 4.0 International (CC-BY 4.0)
Article
A new approach to the formation of systems of linear algebraic
equations for solving ordinary differential equations by the collocation method
L. A. Sevastianov12, K. P. Lovetskiy1, D. S. Kulyabov120
1Peoples' Friendship University of Russia (RUDN University), 6 Miklukho-Maklaya St., Moscow 117198, Russia
Joint Institute for Nuclear Research, 6 Joliot-Curie St., Dubna 141980, Moscow region, Russia
Leonid A. Sevastianov, sevastianov-la@pfur.ru, https://orcid.org/0000-0002-1856-4643, AuthorlD: 6741 Konstantin P. Lovetskiy, lovetskiy-kp@rudn.ru, https://orcid.org/0000-0002-3645-1060, AuthorlD: 584101
Dmitry S. Kulyabov, kulyabov-ds@rudn.ru, https://orcid.org/0000-0002-0877-7063, AuthorlD: 360953
Abstract. A new algorithm for the numerical solution of one-dimensional Cauchy problems and Poisson equations is implemented. The algorithm is based on the collocation method and representation of the solution as an expansion in Chebyshev polynomials. It is proposed instead of the usual approach, which consists in combining all known conditions — differential (the equation itself) and initial / boundary — into one system of approximate linear algebraic equations, to go to the method of solving the problem in several separate stages. First, spectral coefficients are identified that determine the "general" solution of the original problem. The collocation method determines the interpolation coefficients of the derivative of the solution, and thus the expansion coefficients of the solution itself (except for the initial ones). At this stage, the choice of a good basis with discrete orthogonality makes it possible to use very efficient algorithms for finding the desired coefficients. The complexity of reducing the matrix of a system of linear algebraic equations to a diagonal form becomes equivalent to the complexity of multiplying the Chebyshev matrix of coefficients by the vector of the right side of the system. Then the expansion coefficients of the solution itself (except for the first one or two) are obtained by multiplying the known tridiagonal integration matrix (inverse to the Chebyshev differentiation matrix) by the vector of interpolation coefficients of the derivative. At the last stage, considering the initial/boundary conditions select a "particular" desired solution, unambiguously redefining the missing coefficients of the desired expansion.
Keywords: initial boundary value problems, collocation method, Chebyshev polynomials, Gauss-Lobatto sets, numerical stability, discrete orthogonality
Acknowledgements: The work was supported by the RUDN University Strategic Academic Leadership Program.
For citation: Sevastianov L. A., Lovetskiy K. P., Kulyabov D. S. A new approach to the formation of systems of linear algebraic equations for solving ordinary differential equations by the collocation method. Izvestiya of Saratov University. Mathematics. Mechanics. Informatics, 2023, vol. 23, iss. 1, pp. 36-47 (in Russian). https://doi.org/10.18500/1816-9791-2023-23-1-36-47, EDN: BFDVVG
This is an open access article distributed under the terms of Creative Commons Attribution 4.0 International License (CC-BY 4.0)
Введение
Спектральные методы — это класс методов, активно используемых для численного решения различных дифференциальных уравнений [1-4]. Основная идея спектральных методов состоит в том, чтобы представить искомое решение дифференциального уравнения в виде взвешенной суммы определенных «базисных функций» [5] (например, в виде разложения по степенным функциям — ряд Тейлора, или по синусам и косинусам — тригонометрический ряд Фурье), а затем вычислить коэффициенты разложения и тем самым получить приближенное решение.
Спектральные методы обеспечивают более высокую по сравнению с методами конечных элементов сходимость сходимости. При этом их «экспоненциальная сходимость» наиболее быстрая из возможных, когда решение является гладким.
Спектральные методы численного решения обыкновенных дифференциальных уравнений с заданными начальными условиями обычно сводятся к решению СЛАУ, в которую включены как начальные условия, так и условия, обеспечивающие выполнение дифференциальных соотношений [6]. Однако априорное встраивание начальных (граничных) условий в систему линейных уравнений приводит к существенному увеличению заполненности матриц и, следовательно, усложнению алгоритма и метода решения [7] задачи.
Новизна же предлагаемого авторами данной статьи подхода заключается в том, чтобы сначала с помощью устойчивого и простого с вычислительной точки зрения метода выделить класс функций, удовлетворяющих дифференциальному уравнению, тем самым выделив общее решение дифференциального уравнения. Затем рассчитать старшие коэффициенты разложения будущего решения (по вычисленным интерполяционным коэффициентам производной), используя явное представление матрицы интегрирования. И лишь после этого выделить из множества общих решений те частные решения, которые удовлетворяют заданным начальным или граничным условиям.
В работе решение тривиального обыкновенного дифференциального уравнения предлагается разбить на три этапа. Первый этап — вычисление интерполяционных коэффициентов первой либо второй производной умножением матрицы чебышевских коэффициентов на вектор правой части системы. На втором этапе умножение трехдиа-гональной матрицы интегрирования на вектор, полученный в результате выполнения первого этапа, дает коэффициенты разложения «общего» решения задачи. На третьем этапе для определения недостающих коэффициентов, определяющих «частное» решение, привлекаются начальные либо граничные условия.
Каждая подзадача решается устойчиво и просто. Решение первой задачи сводится к умножению вектора правой части на транспонированную матрицу значений функций Чебышева на сетке Гаусса - Лобатто. На следующем шаге решаем СЛАУ с
диагональной положительно определенной матрицей и затем умножаем полученный вектор слева на трехдиагональную чебышевскую матрицу интегрирования [8]. Число начальных (доопределяемых на следующем этапе) коэффициентов соответствует порядку дифференциального уравнения. Они являются решениями элементарных линейных алгебраических уравнений первого (задача Коши) либо второго (уравнение Пуассона) порядка.
Точное решение простейшего обыкновенного дифференциального уравнения при заданном начальном (граничном) условии
правая часть которого не зависит от у и непрерывна по х, имеет единственное решение.
1. Аппроксимация производной. Задача Коши
Вначале рассмотрим задачу определения (восстановления) функции по значениям ее производной на сетке и (одному) дополнительному условию. В такой формулировке задача распадается на три подзадачи:
• полиномиальную интерполяцию производной (вычисление коэффициентов разложения производной в конечный ряд по базисным функциям);
• вычисление коэффициентов решения (кроме первого), которые определяются дифференциальными условиями задачи - умножение обратной по отношению к матрице дифференцирования на вектор интерполяционных коэффициентов производной;
• доопределение коэффициентов решения на основе граничных (или других независимых) условий.
Не ограничивая общности, будем далее считать, что областью определения решения является интервал [—1,1].
Часто аппроксимация непрерывных функций ограничивается отбрасыванием членов ряда Чебышева, величина которых мала [9,10]. В отличие от приближений, получаемых при использовании других степенных рядов, приближение по полиномам Чебышева (обладает свойством почти оптимальности) минимизирует количество слагаемых, необходимых для аппроксимации функции многочленами с заданной точностью. С этим связано также и то свойство, что приближение на основе ряда Чебышева оказывается довольно близко к наилучшему равномерному приближению (среди многочленов той же степени), но вычисляется проще. Кроме того, выбор интерполяционной сетки Гаусса - Лобатто позволяет избавиться от эффекта Гиббса.
Рассмотрим подробнее проблему отыскания коэффициентов полинома р(х), аппроксимирующего решение уравнения (1) в заданном конечном числе точек на интервале [— 1,1]. Спектральный метод коллокации [11] решения задачи основывается на представлении искомой аппроксимирующей функции в виде разложения в конечный ряд
к=о
по полиномам Чебышева первого рода (ж)}§°, которые образуют базис в гильбертовом пространстве функций на отрезке [— 1,1].
у' = f (ж), х ^ х0, х Е [-1,1], У(х о) = Уо,
(1) (2)
(3)
Устойчивое определение коэффициентов разложения решения по известным коэффициентам разложения производной
Продифференцируем равенство (3). Выражение для производной имеет вид
п— 1
р'(х) = Ск Тк (Х) = Ьк Тк (х), X € [-1, 1].
(4)
к=0
к=0
Используя рекуррентные соотношения, которым удовлетворяют чебышевские полиномы первого рода и их производные [2,8], и приравнивая коэффициенты при одинаковых полиномах в (4), приходим [3] к следующей зависимости коэффициентов Ск от Ьк:
Б+Ъ
0 0 1 0 0 4 00 00
00 00 00 00
0
0 0
0
1 6 0
0 0 0 0
0 0 0
0
1
0 0 0 0
0 0 0 0
0 0 0 0 0
0 0 0 0 0
0 0 0 0 0
0 0 0 0 0
0 1/2 п— 3 0 0 0
1/2 п-2 0 1/2 п—2 0 0
1/2 п— 1 0 1/2 - п— 1 0
0 1/2 п
" Ьо - ' Со
Ь1 С1
Ь2 С2
3 -о С3
Ь4 = С4
Ьп—3 С-п—3
Ьп—2 Сп—2
Ьп—1 С-п— 1
_ Ьп . _ С-п _
(5)
Здесь через Б+ обозначена усеченная трехдиагональная матрица с нулевыми элементами на главной диагонали, являющаяся обратной по отношению к матрице дифференцирования Чебышева в спектральном пространстве [8]. Первая строка заполнена нулевыми элементами, поскольку компонента не может быть определена по коэффициентам производной, для ее однозначного определения требуется дополнительное условие, например начальное.
Таким образом, вектор коэффициентов {с0,с\,... ,сп} есть результат умножения простой трехдиагональной матрицы интегрирования Б+ на вектор Ъ. Алгоритм вычисления коэффициентов разложения решения по полиномам Чебышева реализуется по следующим явным формулам:
С1 = Ьо - 62/2, Ск = (Ьк - Ьк+\)/2к, Ск = Ьк-1/2к,
к = 1,
1 < к < п - 1, к = п - 1,п.
(6)
Таким образом, зная коэффициенты разложения функции задачи (1) по полиномам Чебышева 1-го рода, мы можем однозначно восстановить последние п коэффициентов разложения искомой функции по тем же базисным функциям.
Устойчивое определение коэффициентов разложения производной
Итак, первая часть задачи сводится к вычислению коэффициентов {Ъ0,Ъ1, разложения правой части (1) по полиномам Чебышева на интервале [-1,1]
. ,Ьп}
У^Ьк Тк (х) к=0
/ (х).
Метод коллокации заключается в подборе таких коэффициентов {60,&1,... ,ЬП} разложения интерполяционного полинома р'(х), что удовлетворяются (выполняются) следующие равенства при искомых коэффициентах Ьд., к = 0,1,... ,п:
п
ьк Тк ) = f ), 3 = 0,...n, (7)
к=0
в (п +1) точках коллокации {х0, х\,..., хп}.
Последнее утверждение эквивалентно тому, что коэффициенты должны быть решением системы линейных алгебраических уравнений (7) метода коллокаций. В матричной форме это можно записать так:
ТЬ = {. (8)
Выбор точек коллокации может и должен обеспечивать невырожденность системы уравнений (7), для этого достаточно, чтобы все точки сетки были различными, а в остальном их выбор произволен, т. е. решение системы (7) на произвольной сетке интервала [— 1,1] определяет нужную аппроксимацию. Для произвольной сетки матрица Т является полностью заполненной, и численное решение такой системы достаточно трудоемко. Чтобы упростить вид матрицы и ускорить отыскание вектора Ь, воспользуемся свойством дискретной ортогональности чебышевской матрицы Т на сетке Гаусса - Лобатто. В качестве точек коллокации рассмотрим множество х^ = сов('П]/п), ] = 0,... п. Для дальнейшего улучшения свойств системы линейных уравнений, решением которой будет вектор {Ъ0,Ъ\,... ,Ъп}, умножим первое и последнее уравнения из (7) на множитель 1/л/2 . Получим эквивалентную «модифицированную» систему с новой матрицей Т (вместо Т) и вектором f вместо f. Новая система хороша тем, что она обладает свойством «дискретной ортогональности», и умножение ее слева на транспонированную Тт дает диагональную матрицу
~т~ , / п п \ ТТТ = ^ [п,-,...,-,п).
Преобразуем систему (8), умножив ее слева на транспонированную матрицу Тт, получим простое матричное уравнение с диагональной матрицей для определения искомых коэффициентов разложения {Ь0, Ъ\,..., Ьп}:
ТТТ Ь = Тт?. (9)
Обозначая через / произведение Тт на вектор (п, ^,..., § ,п) в (9), выписываем искомые коэффициенты разложения производной решения — функции /(х) — в явном виде
Ъ0 = -, Ъг = 2 ^, Ъ2 = 2 ^, ..., Ьп = ^. (10)
п п п п
В предлагаемом алгоритме дискретная ортогональность матриц ТТ обеспечивает непосредственное вычисление интерполяционных коэффициентов без использования громоздких методов решения СЛАУ.
Вычисление недостающего коэффициента разложения по заданному дополнительному условию — начальному или граничному
Для однозначного определения недостающего коэффициента с0 необходимо привлечение, по крайней мере, еще одного дополнительного условия. Это может быть как граничное условие на левом или на правом конце интервала рассмотрения функции, так и условие прохождения искомой функции через любую заданную точку внутри интервала задания функции.
Таким обазом, рассматриваемый метод позволяет решать в зависимости от вида дополнительного условия как задачу Коши с начальными условиями, так и задачи с краевыми условиями общего вида, требующие, например, применения итерационного метода стрельбы [4] (Shooting method). Рассмотрим подробнее различные частные случаи задания дополнительного (граничного) условия, позволяющего выделить единственное решение задачи (1), (2).
В том случае, когда граничное условие задано на левом конце интервала интегрирования, нулевой коэффициент определяется из уравнения
п
Со + Ckтк (-1) = Уо (11)
к=1
(учитывая, Тк(-1) = (—1)к для к ^ 0) по формуле
п п
со = Уо - ^ СкТк(-1) = уо - ^ Ск(-1)к. (12)
к=1 к=1
Если же дополнительное «граничное» условие задано в произвольной точке интервала интегрирования уь = у(хь), хь Е [-1,1], коэффициент со определяется по формуле
п
со = Уь СкТк (хь). (13)
к=1
На правом конце интервала интегрирования уг = у(1), хг = 1, полиномы Чебышева любого порядка принимают значение, равное 1 (Тк(1) = 1). Поэтому коэффициент со в этом случае определяется по формуле
п п
Со = Уг СкТк (Хг) = Уг Ск. (14)
к=1 к=1
В предлагаемом методе в любом из вариантов задания дополнительных условий недостающий коэффициент вычисляется непосредственно, без необходимости решения сложных систем уравнений или использования итерационных процессов.
2. Аппроксимация второй производной. Уравнение Пуассона
Продемонстрируем практичный и простой устойчивый метод решения одномерных дифференциальных уравнений второго порядка с различными граничными условиями на примере решения простейшего уравнения Пуассона, ограничившись уравнением с постоянными коэффициентами к(х) = 1 и q(x) = 0:
u"(x) = f (х), a ^ ж ^ Ь. (15)
Рассмотрим решение задачи Пуассона с различными граничными условиями, такими как:
• условия Дирихле на обоих концах интервала
и(а) = а, и(Ь) = ¡3;
(16)
• условия Неймана - Дирихле
и' (а) = а, и(Ь) = ¡3;
(17)
• или условия Дирихле - Неймана
и(а) = а, и' (Ь) = ¡3.
(18)
Метод решения задачи (15)-(18) для ОДУ второго порядка состоит из последовательного решения нескольких подзадач. Рассмотрим, как и в случае задачи Коши для ОДУ первого порядка, алгоритм поиска решения, состоящий из трех этапов:
1) вычисление коэффициентов полиномиальной интерполяции второй производной решения - вектора правой части (15) на сетке Гаусса - Лобатто;
2) вычисление тех коэффициентов решения (кроме двух первых), которые определяются из дифференциальных условий задачи (позволяющих решению удовлетворять дифференциальным условиям) — умножение обратной матрицы [8] на вектор интерполяционных коэффициентов;
3) доопределение коэффициентов решения на основе граничных (или иных независимых дополнительных) условий.
Как и в предыдущем разделе, будем представлять приближенное решение в виде конечного ряда (3) — разложения по базису из ортогональных полиномов Чебышева. И первый этап решения задачи Пуассона практически полностью совпадает с первым этапом решения задачи Коши. Изменяется лишь число вычисляемых на этом этапе коэффициентов разложения решения, обеспечивающих выполнение дифференциальной части (15) решения задачи Пуассона.
Продифференцируем дважды функцию (3). Выражение для второй производной имеет вид
Используя рекуррентные соотношения, которым удовлетворяют чебышевские полиномы первого рода и их производные [2,8], и приравнивая коэффициенты при одинаковых полиномах в (19), приходим [2] к следующей зависимости коэффициентов
Сг, г = 2, 3,... ,п от Ьк, г = 0,... ,п:
где Б+ — обобщенная обратная матрица по отношению к матрице дифференцирования Чебышева в спектральном пространстве (см. уравнение (5)). Иными словами, вектор коэффициентов со,с\,... ,сп является результатом двукратного умножения простой трехдиагональной матрицы интегрирования на вектор Ь0,Ъ\,...,ЪП. Алгоритм вычисления коэффициентов разложения решения по полиномам Чебышева реализуется по дважды примененным явным формулам (6) в случае вычисления второй производной.
Отличие от метода решения задачи Коши появляется на третьем этапе решения задачи Пуассона — определении двух первых коэффициентов разложения искомого
(19)
Б+Б+Ъ = с,
(20)
решения по полиномам Чебышева. Рассматриваемый метод позволяет решать в зависимости от вида дополнительных условий как задачи Коши с начальными условиями, так и задачи с краевыми условиями общего вида.
В случае граничных условий Дирихле (граничные условия первого рода): р(а) = а, р(Ь) = /3, определение остающихся пока еще неизвестными коэффициентов со, с1 сводится к решению системы двух уравнений, которыми могут быть, например, уравнения, определяющие поведение решения в граничных точках: х = ±1
п п
Со + С1 Т1,о(-1) + ^ скТм(-1) = а, Со + сгТх,п(1) + ^ скТм(1) = (21)
к=2 к=2
Если дополнительно учесть тот факт, что полиномы Чебышева первого рода принимают на границе интервала значения Тк^(±1) = ±1, ],к = 0,..., то решение системы (21) можно записать в явном виде
Со = 1 (а + р - £ сИ , С1 = 1 и - а - £ сИ . (22)
\ к=2, к — четное / \ к=3, к — нечетное /
В случае смешанных условий Неймана - Дирихле (граничные условия второго и первого рода): р'(а) = а, р(Ь) = коэффициенты с0,сг определяются по формулам
п п
сг = а - £(-1)^*, со = $ - С1 - Е ск, (23)
к=2 к=2
а в случае Дирихле - Неймана —
С1 = /3 - ^ к2Ск, со = а - С1 - ^(-1)к+1^2Ск. (24)
к=2 к=2
Третий тип граничных условий встречается при моделировании процессов охлаждения и носит название условий Робена или Ньютона [12]. Их использование в качестве граничных условий приводит лишь к несколько более сложным соотношениям для определения недостающих коэффициентов по сравнению с условиями (21)-(24).
Решение 1-0 уравнения Пуассона. Сравнение с точным решением
Разработанный выше метод решения уравнения Пуассона (15), (16) рассмотрим на конкретном примере. Считаем, что наше уравнение имеет вид
и"(х) = -(2х + 5)ех, -1 ^ ж ^ 2, (25)
а с помощью «случайно» известного точного решения иехаа = (2х + 1)ех граничные условия Дирихле определяются выражениями
Ц-1) = -1/е, и(2) = 5е2. (26)
Сравним отклонения решения, полученного с помощью изложенного выше алгоритма, от точного решения. На рисунке приводятся график решения задачи и вычисленные поточечные значения невязки в зависимости от числа точек аппроксимации п на всем интервале интегрирования.
argument
в / c
argument
г / d
Точность восстановления второй производной (решение уравнения Пуассона) в зависимости от количества п точек аппроксимации: а) ueXact = (2х + 1)еж; б) residual — невязка при п = 10; в) residual — невязка при п = 12; г) Residual — невязка при п = 15
Figure. The accuracy of the recovery of the second derivative (solution of the Poisson equation) depending on the number of n approximation points: a) uexact = (2x + 1)еж; b) residual at n = 10;
c) residual at n = 12; d) residual at n = 15
Заключение
Среди методов решения начальных и граничных задач для линейных ОДУ первого и второго порядка есть методы, основанные в основном на локальной аппроксимации решения, которые сразу используют начальное приближение (граничные условия) при решении дифференциальных уравнений. Это методы типа Эйлера, Рунге - Кутты и т.д. Другие методы, основанные на аппроксимации решения с помощью глобальных функций [1-3,6,10], основываются на построении таких систем уравнений, которые одновременно включают в себя как начальные (граничные) условия, так и условия, задающие поведение производных искомого решения.
В нашей работе в рамках спектрального метода коллокаций решение основной задачи представляется в виде последовательного решения нескольких подзадач. Вначале выделяется множество решений, которое удовлетворяет дифференциальному уравнению. Однако оно не обязательно удовлетворяет начальным (граничным) условиям.
Учет начальных (граничных) условий осуществляется на последнем этапе решения исходной задачи и фактически сводится к решению линейного уравнения с одним неизвестным коэффициентом для ОДУ первого порядка или с двумя для уравнения второго порядка.
Решение первой задачи сводится к умножению матрицы значений функций Чебы-шева на сетке Гаусса - Лобатто на вектор значений функции, задающей правую часть
исходного дифференциального уравнения для определения интерполяционных коэффициентов разложения производной решения. Далее умножение двухдиагональной матрицы интегрирования [3] на вектор этих коэффициентов приводит к получению всех коэффициентов искомого решения, кроме первого (двух первых). На последнем этапе использование начального (граничного) условия дает возможность определить и первый коэффициент (два первых коэффициента) разложения решения.
Новый подход, основанный на разделении задачи решения ОДУ первого и второго порядков на отдельные подзадачи, является весьма перспективным. Авторы продолжат разработку подходов и алгоритмов применения многостадийного спектрального метода для решения начальных и краевых задач с дифференциальными уравнениями более высоких порядков.
Список литературы
1. Boyd J. P. Chebyshev and Fourier Spectral Methods: Second Revised Edition. Dover Books on Mathematics, 2013. 668 p.
2. Mason J. C., Handscomb D. C. Chebyshev Polynomials. Chapman and Hall/CRC Press, 2002. 360 p. https://doi.org/10.1201/9781420036114
3. Fornberg B. A Practical Guide to Pseudospectral Methods. New York : Cambridge University Press, 1996. 231 p. https://doi.org/10.1017/CB09780511626357
4. Press W. H., Teukolsky S. A., Vetterling W. T., Flannery B. P. Numerical Recipes: The Art of Scientific Computing. 3rd ed. New York : Cambridge University Press, 2007. 1235 p.
5. Shen J., Tang T., Wang L.-L. Spectral Methods: Algorithms, Analysis and Applications. Berlin ; Heidelberg : Springer, 2011. 472 p. (Springer Series in Computational Mathematics, vol. 41). https://doi.org/10.1007/978-3-540-71041-7
6. Olver S., Townsend A. A Fast and Well-Conditioned Spectral Method // SIAM Review. 2013. Vol. 55, iss. 3. P. 462-489. https://doi.org/10.1137/120865458
7. Chandrasekaran S., Gu M. Fast and Stable Algorithms for Banded Plus Semiseparable Systems of Linear Equations // SIAM Journal on Matrix Analysis and Applications. 2003. Vol. 25, iss. 2. P. 373-384. https://doi.org/10.1137/S0895479899353373
8. Amiraslani A., Corless R. M., Gunasingam M. Differentiation matrices for univariate polynomials // Numerical Algorithms. 2020. Vol. 83, iss. 1. P. 1-31. https://doi.org/10. 1007/s11075-019-00668-z
9. Zhang X., Boyd J. P. Asymptotic coefficients and errors for Chebyshev polynomial approximations with weak endpoint singularities: Effects of different bases // Science China Mathematics. 2023. Vol. 66, iss. 1. P. 191-220. https://doi.org/10.1007/s11425-021-1974-x
10. Boyd J. P., Gally D. H. Numerical experiments on the accuracy of the Chebyshev-Frobenius companion matrix method for finding the zeros of a truncated series of Chebyshev polynomials // Journal of Computational and Applied Mathematics. 2007. Vol. 205, iss. 1. P. 281-295. https://doi.org/10.1016/j-.cam.2006.05.006
11. Dutykh D. A Brief Introduction to Pseudo-spectral Methods: Application to Diffusion Problems. 2019. 55 p. URL: https://arxiv.org/pdf/1606.05432 (дата обращения: 30.05.2022).
12. Dawkins P. Differential Equations. 2018. 524 p. URL: https://tutorial.math.lamar.edu/ Classes/DE/DE.aspx (дата обращения: 30.05.2022).
References
1. Boyd J. P. Chebyshev and Fourier Spectral Methods: Second Revised Edition. Dover Books on Mathematics, 2013. 668 p.
2. Mason J. C., Handscomb D. C. Chebyshev Polynomials. Chapman and Hall/CRC Press, 2002. 360 p. https://doi.org/10.1201/9781420036114
3. Fornberg B. A Practical Guide to Pseudospectral Methods. New York, Cambridge University Press, 1996. 231 p. https://doi.org/10.1017/CB09780511626357
4. Press W. H., Teukolsky S. A., Vetterling W. T., Flannery B. P. Numerical Recipes: The
Art of Scientific Computing. 3rd ed. New York, Cambridge University Press, 2007. 1235 p.
5. Shen J., Tang T., Wang L.-L. Spectral Methods: Algorithms, Analysis and Applications. Springer Series in Computational Mathematics, vol. 41. Berlin, Heidelberg, Springer, 2011. 472 p. https://doi.org/10.1007/978-3-540-71041-7
6. Olver S., Townsend A. A Fast and well-conditioned spectral method. SIAM Review, 2013, vol. 55, iss. 3, pp. 462-489. https://doi.org/10.1137/120865458
7. Chandrasekaran S., Gu M. Fast and Stable Algorithms for Banded Plus Semiseparable Systems of Linear Equations. SIAM Journal on Matrix Analysis and Applications, 2003, vol. 25, iss. 2, pp. 373-384. https://doi.org/10.1137/S0895479899353373
8. Amiraslani A., Corless R. M., Gunasingam M. Differentiation matrices for univariate polynomials. Numerical Algorithms, 2020, vol. 83, iss. 1, pp. 1-31. https://doi.org/10. 1007/s11075-019-00668-z
9. Zhang X., Boyd J. P. Asymptotic coefficients and errors for Chebyshev polynomial approximations with weak endpoint singularities: Effects of different bases. Science China Mathematics, 2023, vol. 66, iss. 1, pp. 191-220. https://doi.org/10.1007/s11425-021-1974-x
10. Boyd J. P., Gally D. H. Numerical experiments on the accuracy of the Chebyshev -Frobenius companion matrix method for finding the zeros of a truncated series of Chebyshev polynomials. Journal of Computational and Applied Mathematics, 2007, vol. 205, iss. 1, pp. 281-295. https://doi.org/10.1016/j-.cam.2006.05.006
11. Dutykh D. A Brief Introduction to Pseudo-spectral Methods: Application to Diffusion Problems, 2019. 55 p. Available at: https://arxiv.org/pdf/1606.05432 (accessed 30 May 2022).
12. Dawkins P. Differential Equations, 2018. 524 p. Available at: https://tutorial.math.lamar. edu/Classes/DE/DE.aspx (accessed 30 May 2022).
Поступила в редакцию / Received 14.06.2022 Принята к публикации / Accepted 26.09.2022 Опубликована / Published 01.03.2023
