Вычисление коэффициентов разложения решения задачи Коши в ряд по многочленам Чебышева Текст научной статьи по специальности «Математика»

5. Михалев А.А., Михалев А.В., Чеповский А.А., Шампаньер К. Примитивные элементы свободных неассоциативных алгебр // Фунд. и прикл. матем. 2007. 13, № 5. 171-192.

6. Mikhalev A.A., Yu J.-T. Primitive, almost primitive, test, and Д-primitive elements of free algebras with the NielsenSchreier property //J. Algebra. 2000. 228. 603-623.

7. Климаков А.В., Михалев А.А. Почти примитивные элементы свободных неассоциативных алгебр малых рангов // Фунд. и прикл. матем. 2012. 17, № 1. 127-141.

Поступила в редакцию 17.10.2011 После доработки 02.11.2011

УДК 519.622

ВЫЧИСЛЕНИЕ КОЭФФИЦИЕНТОВ РАЗЛОЖЕНИЯ РЕШЕНИЯ ЗАДАЧИ КОШИ В РЯД ПО МНОГОЧЛЕНАМ ЧЕБЫШЕВА

О. Б. Арушанян1, Н. И. Волченскова2, С. Ф. Залеткин3

Предложены два способа определения начального приближения коэффициентов разложения решения задачи Коши для обыкновенных дифференциальных уравнений в ряды по смещенным многочленам Чебышева первого рода. Данное приближение предназначено для применения в аналитическом методе интегрирования обыкновенных дифференциальных уравнений с использованием ортогональных разложений.

Ключевые слова: обыкновенные дифференциальные уравнения, приближенные аналитические методы интегрирования, численные методы интегрирования, ортогональные разложения, смещенные ряды Чебышева, квадратурные формулы Маркова.

Two approaches are proposed to determine an initial approximation for the coefficients of an expansion of the solution to a Cauchy problem for ordinary differential equations in the form of series in shifted Chebyshev polynomials of the first kind. This approximation is used in an analytical method to solve ordinary differential equations using orthogonal expansions.

Key words: ordinary differential equations, approximate analytical methods of integration, numerical methods of integration, orthogonal expansions, shifted Chebyshev polynomials, Markov quadrature formulas.

Рассматривается задача Коши для нормальной системы

y'(x) = f (x,y(x)), y(xо) = yo, xo ^ x ^ xo + X, (1)

и задача Коши для канонической системы M обыкновенных дифференциальных уравнений второго порядка

y"(x) = f (x, y(x),y'(x)), y(xо) = yo, y'(xo) = y0, xo < x < xo + X. (2)

В [1] предложен приближенный аналитический метод интегрирования уравнений (1), (2), основанный на разложении решения y(x) в ряд Фурье по смещенным многочленам Чебышева первого рода. Частичная сумма такого ряда принимается в качестве приближения к решению задачи Коши (1), (2). В [1] приведены соотношения, которые связывают коэффициенты Чебышева решения y(x) (и его производной y'(x) для задачи (2)) с коэффициентами Чебышева правой части системы (1) (или (2)), взятой на решении

1 Арушанян Олег Багратович — доктор техн. наук, проф., зав. лаб. Научно-исслед. вычисл. центра МГУ, e-mail: arush@srcc.msu.ru.

2 Волченскова Надежда Ивановна — канд. физ.-мат. наук, ст. науч. сотр. Научно-исслед. вычисл. центра МГУ, e-mail: nad1946@mail.ru.

3 Залеткин Сергей Федорович — канд. физ.-мат. наук, ст. науч. сотр. Научно-исслед. вычисл. центра МГУ, e-mail: iraz@srcc.msu.ru.

задачи Коши. Для отыскания коэффициентов разложения правой части в ряд Фурье-Чебышева на основе квадратурных формул Маркова [2-5] в [1] построены уравнения и описан итерационный процесс их решения.

Настоящая статья посвящена изучению вопроса о выборе начального приближения коэффициентов Чебышева правой части системы для данного итерационного процесса. Здесь мы рассмотрим два способа определения начального приближения. Будем предполагать существование у правой части дифференциального уравнения непрерывных частных производных такого порядка, который гарантировал бы справедливость приводимых оценок и правомерность используемых преобразований.

1. Первый способ определения начального приближения. Опишем первый подход в применении к частичному сегменту [xo,Xo + h], 0 < h ^ X, на котором приближенное решение строится в виде конечной суммы ряда.

Обозначим k-ю частичную сумму смещенного ряда Чебышева для правой части системы (1) или (2) на сегменте [xo ,Xo + h] через Jk. Пусть U — частичная сумма смещенного ряда Чебышева для решения задачи (1) или (2), U' — частичная сумма смещенного ряда Чебышева для производной y'(x) решения задачи (2) на том же сегменте, а* — значения коэффициентов смещенного ряда Чебышева и Т*(а) — смещенные многочлены Чебышева первого рода.

В [1, 5] для вычисления коэффициентов Чебышева была использована квадратурная формула Маркова с одним наперед заданным узлом. Здесь мы приведем изложение на примере формулы Маркова с двумя наперед заданными узлами [6].

Рассмотрим сначала нормальную систему (1). Положим

т^о Ш = ^a*0[f(x0,yo)} = f(x0,yo)

и определим приближенные значения первых двух коэффициентов Чебышева для решения:

h 1 h

a*i[U] = -a*0[Jk], -a*0[U]=y0 + -a*0[Jk}.

Далее вычислим в k + 2 узлах квадратурной формулы

1 ( 3 п \

Xj = х0 + atjh, ао = 0, aj = - (Ч + cos -j—pyj, j = 1, ak+i = 1, (3)

значения приближенного решения и правой части:

U(x°) = i а*0[и]П(a,) + at[Umaj) = ± a*0[U] + a\[U] cos

Ф(ао) = / (хо,Уо), Ф(а ) = / ; = 1 ,...,к + 1.

В качестве нулевого приближения для а*[,1к], % = 0,1, ... ,к, возьмем значения, определяемые по квадратурной формуле Маркова с двумя наперед заданными узлами ао = 0 и ак+\ = 1, с к нефиксированными узлами <х/, = 1, ... , к, и весовой функцией 1/^/а(1 — а):

О к + 1 /1^1 о к

- = ^ 4(0) + ш *<и + ш £«. ш *«*>■ <4>

3=0 3=1

,, т ,, 1 1 Здесь символ определен формулой = — щ + + ... + ат-\ -\— ат, т > I. Эти значения

3=1

имеют погрешность 0(Н2), Н ^ 0.

Теперь рассмотрим каноническую систему (2). В этом случае последовательность вычислений будет следующей. Положим

2 а*оШ = 2 ао[/(жо,Уо,Уо)] = /(%о,Уо,Уо)

и определим приближенные значения первых двух коэффициентов Чебышева для производной от решения:

Н 1 Н

а* Р'] = ^а*0 [Л], - а*0 [II'} =у'0 + ^а*0 [Л]

и приближенные значения первых трех коэффициентов Чебышева для решения:

Ь? ~ Н Н 1 Н 3Н2

а*2[и} = —а*0Щ, а1[и} =-(у'0 + -а*0Щ), -а*0[и}=уо + -у'о + — а*0Ш.

По этим коэффициентам вычисляем в к + 2 узлах (3) значения приближенного решения и его производной:

и'{$) = I аЪ\и'Щ{аэ) + а1[и'}Т?(ау) = I а*0[и'} + аЦи'} сов ^ 1

2^ ^^ , ^ ^2^ > ' —к + 1'

и(х°3) = - аЪ[иЩ{а3) + аЦЩЩа,) + а*2[и]П(а3) = = I а*0[и] + аЦЩ со8 + а*2р] со8 ■

2 ^ ^ к + 1 ^ к + 1 ' а затем и правую часть системы (2):

Ф(ао ) = / (хо,Уо,Уо), ) = / (х0)), ; = 1, ...,к + 1.

В качестве нулевого приближения для а*[.к], г = 0,1,... ,к, возьмем значения, определяемые по квадратурной формуле Маркова (4). Эти значения, как и в случае нормальной системы (1), имеют погрешность 0(Н2), Н ^ 0.

Дальнейшее уточнение приближенных значений коэффициентов выполняется с помощью итерационного процесса так, как это описано в [1]. Для погрешности полученных в результате итерационного процесса приближенных значений коэффициентов Чебышева а* [.к] справедлива асимптотическая оценка

а*[Ф] - а* ]= 0(Нк+2), Н ^ 0, (5)

где для системы (1)

Ф(а) = ¥(х0 + аН) = /(х0 + аН, у(х0 + аН)), 0 ^ а ^ 1, (6)

а для системы (2)

Ф(а) = ¥(х0 + аН) = /(х0 + аН,у(х0 + аН),у'(х0 + аН)), 0 ^ а ^ 1. (7)

В рассмотренном способе исходными данными являются значения решения и его производной на левом конце промежутка интегрирования, поэтому приведенный алгоритм может быть применен на любом частичном сегменте [хп , хп + Н].

2. Второй способ определения начального приближения. Описанный ниже способ опирается на экстраполяцию коэффициентов Чебышева с предыдущего сегмента на следующий. Покажем, как коэффициенты Чебышева разложения функции Ф(а), заданной формулами (6) и (7), могут быть использованы для получения коэффициентов Чебышева разложения функции Ф:(а) = ¥(х1 + аН*), где 0 ^ а ^ 1, х1 = х0 + Н, Н* — длина следующего частичного сегмента [х1,х1 + Н*] и ¥(х) — правая часть системы (1) (или (2)), взятая на решении задачи Коши на сегменте [х1,х1 + Н*].

Обратимся к частичной сумме смещенного ряда Чебышева функции Ф(а) с приближенными коэффициентами. Пусть а = — (£ + 1) и Ф(£) = Фтогда а*[Ф(а)] = , где ец — коэффициенты разложения функции в ряд по несмещенным многочленам Чебышева. С учетом (5) имеем

к к

Бк(а, Ф) = .(а) = ^'а*.]7?(а) = Е'а#(*)№(*) + 0(Нк+2), (8)

г=0 г=0

т 1

где символ ^^ определен формулой ^^ а^ = — щ + + ... + ат, т ^ I. Используя выражение

з =1

коэффициента Чебышева разложения функции через производную этой функции

1 Ф М(6) ъ[т\ = ^ТТ —р2 , 0 е [-1,1],

и упорядочивая сумму по степеням переменной t, получим

к - 1 Ф<flfe)™ , А - 1

АЛ*) = E'jPT^W, + O(h^) = ±(£)»<*w> + o(ft«)

i=0 ' i=0 '

i= i' 2

где Xi/2 = Xo + h/2. Аналогично можно представить частичную сумму ряда Чебышева функции Ф^а):

к

Sk(t, Ф0 = ¿) (i (у JV%3/2) + + 0(h*k+2), (10)

1 / h\*

Ti(t) + O(hk+2) =

где Хз/2 = Ж1 + Н*/2.

Для дальнейшего определим следующие матрицы: Р — верхняя треугольная матрица Паскаля с элементами = С3—1, 1 ^ % ^ к + 1, 1 ^ j ^ к + 1, % ^ j; Т и I — диагональные матрицы (к + 1)-го

порядка: Т = {1,1 + £, (1 + £)2, ... , (1 + £)к} и I = {1, ш, ш2, ... , шк}, где £ = Н*/Н и ш = {/(1 + £). Введем векторы длиной к + 1:

Н .. 1 /Н\2 ,, 1 /Н\к 4 Т

(h 1 / Ъ\2 1 //} \к \

F{xl/2),-F'{xl/2),-[-) F"{x 1/2),...,-(-) f(%1/2)J

Ж^з/2) = (^з/2),у F\x3/2),l(^)2F"(x3/2), ... kF(k\x3/2)y.

Последовательно применяя линейные преобразования, задаваемые матрицами Т, Р и I, к вектору 2(Х1/2), приходим к соотношениям

2 'к^ 2

Н + Н* , , 1 /Н + Н^к 4 т

— iГ(*3/2),...,й( —

IPTZ(Жl/2)= 2 (хз/2)+ 0(^

Точные компоненты вектора 2(Ж1/2) нам неизвестны, но приближения к ним могут быть получены из (9). Применяя указанные выше преобразования к вектору, составленному из этих приближений, мы найдем компоненты вектора 2(Ж3/2) с погрешностями 0(Н*),0(Н*2), ... ,0(Н*к+1), а именно получим величины

* = ^(у)^(<)(*з/2 ) + 0(Г+1).

Далее составим полином переменной Ь с коэффициентами, равными числам и выполним с ним те же преобразования, что были использованы при выводе формул (8)—(10), но осуществим их в обратном порядке:

к к -1 7 ^ к л

= (y)>W(-3/2) + 0(r+1))f = £ (¿Ф(1°(0) +о(нпу =

i=0 ' i=0

¿'(¿Г I + 0{к*г+1))т) = ¿' (a#i(i)] + 0{к*г+1))т)

,i'.\2) v °/2/ v V ^ w

i=0 i=0 i=0

2 1 , ^г)(и) + 04/V ) ]'Ji(t) = > la.

i=0 2 i=0 к

:£'( а*[Фх(а)]+ O(h*i+1 )) Т*(а).

i=0

Таким образом, в итоге мы имеем коэффициенты Чебышева функции с погрешностями

O(h*),O(h* ), ... ,O(h* + ), которые и принимаем за нулевое приближение для неизвестных a*[Jk] на сегменте [xi,xi + h*]. Дальнейшее уточнение значений коэффициентов a*[Jk] выполняется с помощью итерационного процесса так, как это описано в [1].

3. Примеры. 1) Интегрируется система обыкновенных дифференциальных уравнений

8 2 2 8 1 2/1 = —g 2/1 — з г/2, г/2 = з г/1-д г/2, ш(о) = i, у{( о) = -, у2(о) = о, у'2( о) = i. (11)

У системы колеблющееся решение с быстро возрастающей амплитудой: yi(x) = ex/3 cos ж, У2(ж) = ex/3 sin ж. Задача решалась методом рядов Чебышева на сегменте [0,Xp], Xp = 25. Задавалось разбиение интервала на несколько частичных сегментов длиной h; на каждом таком сегменте решение представлялось в виде (k + 2)-й частичной суммы ряда Чебышева, а производная — в виде (k + 1)-й частичной суммы ряда при k = 30. Вычисления проводились с 16 значащими цифрами. Значения h, способ выбора начального приближения для коэффициентов Чебышева правой части системы (11), количество верных цифр в приближенных значениях обеих компонент решения y(Xp) и производной y' (Xp), вычисленных в конце интервала, а также количество вычислений правой части системы на [0, Xp] приведены в табл. 1.

Таблица 1

h Количество верных цифр для Í/1 (Xf), у2(Xf) Количество верных цифр для vÍ(Xf), V2(Xf) Количество вычислений правой части

Первый способ Второй способ Первый способ Второй способ Первый способ Второй способ

1,0 15 15 15 15 16 15 16 15 24775 24055

2,0 15 14 15 14 15 15 15 15 12883 12523

3,0 15 14 14 14 15 15 16 14 8919 8679

4,0 16 14 14 14 15 15 15 15 6937 6757

5,0 14 14 15 13 14 15 14 14 4955 4835

6,0 14 13 14 15 14 14 15 14 4955 4835

7,0 14 13 13 12 14 14 13 13 3964 3874

8,0 14 12 12 12 13 13 13 13 3964 3874

10,0 13 11 12 11 12 13 13 12 2973 2913

15,0 10 10 10 10 10 10 10 10 1982 1952

Для (11) выполнены расчеты одношаговыми и многошаговыми методами [2-4, 7, 8] с постоянным шагом Н, принимающим указанные в табл. 1 значения. Результаты, полученные методами Адамса и

Штермера пятого порядка типа

Таблица 2

Номер коэффициента Коэффициенты Чебышева для yi{x)

Первый способ Второй способ

0 0,5534836873275461D+04 0,5534836873275460D+04

1 0,1450726099625513D+04 0,1450726099625513D+04

2 — 0,7289200606349429D+02 — 0,7289200606349432D+02

3 — 0,2118248114775613D+02 — 0,2118248114775612D+02

4 — 0,4675386269120383D+00 — 0,4675386269120360D+00

5 0,5821224677364755D-01 0,5821224677364673D-01

6 0,2710917034657765D-02 0,2710917034657737D-02

7 — 0,3170047963506896D-04 — 0,3170047963482767D-04

8 — 0,4030145656278271D-05 — 0,4030145656548678D-05

9 — 0,4417775897784090D-07 — 0,4417775908785338D-07

10 0,2375651196801003D-08 0,2375651387409987D-08

11 0,6396422013194907D-10 0,6396423072230448D-10

12 — 0,3612562281282079D-12 — 0,3613642944633837D-12

13 — 0,3317039480993550D-13 — 0,3315759487000173D-13

14 — 0,2973936773353770D-15 — 0,2488375043259698D-15

15 0,4090733518293968D-16 0,2593140546042205D-16

16 0,1684459544377387D-16 0,2535239366786601D-16

32 0,3648143580180812D-17 0,3714266797353567D-17

предиктор-корректор для решения у(Хр) и производной у'(Хр) в конце интервала интегрирования для всех значений шага Н из табл. 1, не имеют ни одной верной цифры. При интегрировании с постоянным шагом трехстадийным неявным методом Рунге-Кутты шестого порядка результаты для у(Хр), у ' (Хр) имеют от трех до одной верной цифры при Н = 1,2,3, а при интегрировании классическим методом Рунге-Кутты четвертого порядка — только одну верную цифру при Н = 1. Для остальных значений шага Н из табл. 1 результаты, полученные явным и неявным методами Рунге-Кутты, вообще не содержат ни одной верной цифры.

В табл. 2 даны с удвоенной точностью коэффициенты разложения

первой компоненты решения у 1(х) на сегменте [24, 25] в виде частичной суммы смещенного ряда Чебы-шева при к = 30:

32

У1(24 + а) 'а*МТ», 0 < а < 1.

г=0

Эти коэффициенты получены с использованием описанных выше способов определения начального приближения. Видно, что последовательность коэффициентов достаточно регулярно стремится к нулю и ряд Чебышева на данном сегменте быстро сходится.

Кроме того, выполнено решение системы (11) одношаговыми и многошаговыми методами с автоматическим выбором шага интегрирования. Количество верных цифр, соответствующее наилучшей фактически достигнутой точности решения У(Хр) в конце интервала интегрирования, и количество вычислений правой части системы, которое потребовалось для достижения этой точности, представлены для каждого метода в табл. 3.

Таблица 3

Количество верных цифр для Количество вычислений

Метод компонент решения и производной правой части

УЛХР) У2(ХР) у'ЛХр) У2 (.Хр)

Классический метод 11 9 и 11 228424

Рунге-Кутты

Метод Мерсона 11 10 и 11 228460

Метод Адамса 11 10 и 11 32913

Метод Штермера 12 11 12 12 18163

Неявный метод 12 11 12 12 6951

Рунге-Кутты

Заметим, что фактически достигнутая точность и количество вычислений правой части дифференциального уравнения при решении задачи методами Адамса и Штермера зависят от допустимого наибольшего значения шага интегрирования НМАХ; приведенные в табл. 3 данные получены при НМАХ = 0,1. Следует также отметить, что для примера (11) на каждом шаге интегрирования соответствующая система алгебраических уравнений неявного метода Рунге-Кутты, которая является линейной, решалась прямым, а не итерационным методом.

2) Интегрируется дифференциальное уравнение

У = —У

У(-9) =

1

0,3025

-9 < х < 1.

Решением уравнения является рациональная функция

(12)

Таблица 4

У(х) =

1

-, с = 3,05,

Метод Абсолютная погрешность для у(Хр) Количество вычислений правой части

Метод рядов 0,278 • 10~1Ь 7810

Чебышева

Классический метод 0,694 • 10~1Ь 143348

Рунге-Кутты

Метод Мерсона 0,153 • 10~1Ь 116955

Метод Адамса 0,290 • 10~14 6839

х + с2 ;

график которой — ветвь гиперболы. Задача решалась на сегменте [—9, Хр], Хр = 1, методом рядов Чебышева, при этом задавалось разбиение интервала интегрирования на десять частичных сегментов длиной Н = 1 и на каждом из них решение представлялось в виде (к + 1)-й частичной суммы смещенного ряда Чебышева при к = 30. Вычисления проводились с 16 значащими цифрами. Абсолютная погрешность для У(Хр) и количество вычислений правой части уравнения (12) приведены в табл. 4. Кроме того, указанное уравнение решалось одношаговыми и многошаговым методами с автоматическим выбором шага интегрирования. Наилучшая фактически достигнутая точность для У(Хр) и количество вычислений правой части уравнения, потребовавшееся при получении этой точности, приведены для каждого метода в табл. 4.

В табл. 5 содержатся с удвоенной точностью коэффициенты разложения решения у(х) на сегменте [0,1] в виде частичной суммы смещенного ряда Чебышева при к = 30:

31

У(х) Ч*[У]Т*(х),

0 < х < 1.

1=0

2

Второй столбец табл. 5 включает в себя приближенные значения коэффициентов, полученные с помощью метода рядов, а третий столбец — точные значения коэффициентов, которые определялись из известного ряда Чебышева рациональной функции [9]

1

8р2

(-1)У*ТТ(*),

X + с2

1 - р4

1 г=0

р = л/1 + С2 — С,

С > 0.

Таблица 5

Приближенные значения коэффициентов а0 [у], а* [у] имеют абсолютные погрешности 0,250 ■ 10-15 и 0,130 ■ 10-16. Абсолютные погрешности остальных коэффициентов существенно меньше 10-17. Из табл. 5 видно, что последовательность коэффициентов достаточно регулярно стремится к нулю и смещенный ряд Чебы-шева на данном сегменте быстро сходится.

Таблица 5 наглядно иллюстрирует действительную способность рассматриваемого аналитического метода правильно воспроизводить ряды Чебышева для решений обыкновенных дифференциальных уравнений.

Из приведенных результатов можно сделать следующий вывод. Представление решения обыкновенного дифференциального уравнения в виде частичной суммы ряда Чебышева позволяет вычислять приближение к решению с такой высокой точностью, которая на практике может оказаться недостижимой для одношаговых методов типа Рунге-Кутты и многошаговых методов Адамса и Штермера при той же разрядной сетке, поскольку эта точность требует для указанных методов столь малых размеров шага интегрирования, что эти шаги выходят за границу их реальной области асимптотики [8]. Как следует из табл. 1, 3, 4, более высокая точность может быть достигнута с помощью метода рядов Чебышева даже с меньшим числом вычислений правой части, чем в других рассмотренных методах.

Работа выполнена при финансовой поддержке РФФИ, грант № 10-01-00297.

Номер коэффициента Коэффициенты Чебышева

приближенные точные

0 0,2042955207682877Б+00 0,2042955207682874Б+00

1 — 0,52136846622785230-02 -0,5213684662278510Б-02

2 0,1330548396531354Б-03 0,1330548396531349Б-03

3 -0,3395600520915374Б-05 -0,3395600520915382Б-05

4 0,866567719573306Ш-07 0,8665677195732996Б-07

5 — 0,221150753155051Ш-08 -0,2211507531527976Б-08

6 0,5643835391680775Б-10 0,5643835388206232Б-10

7 -0,1440324170350716Б-11 -0,1440324187689344Б-11

8 0,3675750896458252Б-13 0,3675751723691420Б-13

9 -0,9380651486359762Б-15 -0,9380631700628355Б-15

10 0,2395198144517287Б-16 0,2393966124960754Б-16

11 -0,6140044377730617Б-18 -0,6109475342770058Б-18

12 0,4739744909928849Б-20 0,1559156939387654Б-19

31 0,8041262883838078Б-21 -0,8389232458235244Б-50

СПИСОК ЛИТЕРАТУРЫ

1. Арушанян О.Б., Волченскова Н.И., Залеткин С.Ф. О применении ортогональных разложений для приближенного интегрирования обыкновенных дифференциальных уравнений // Вестн. Моск. ун-та. Матем. Механ. 2010. № 4. 40-43.

2. Бабенко К.И. Основы численного анализа. М.: Наука, 1986.

3. Мысовских И.П. Лекции по методам вычислений. СПб.: Изд-во СПб. ун-та, 1998.

4. Ильин В.П., Кузнецов Ю.И. Алгебраические основы численного анализа. Новосибирск: Наука, 1986.

5. Арушанян О.Б., Залеткин С.Ф. О применении формулы численного интегрирования Маркова в ортогональных разложениях // Вестн. Моск. ун-та. Матем. Механ. 2009. № 6. 18-22.

6. Залеткин С.Ф. Формула численного интегрирования Маркова с двумя фиксированными узлами и ее применение в ортогональных разложениях // Вычислительные методы и программирование. 2005. 6, № 1. 141-157.

7. Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы. М.: БИНОМ, 2007.

8. Арушанян О.Б., Залеткин С.Ф. Численное решение обыкновенных дифференциальных уравнений на Фортране. М.: Изд-во МГУ, 1990.

9. Пашковский С. Вычислительные применения многочленов и рядов Чебышева. М.: Наука, 1983.

Поступила в редакцию 08.11.2011