Методы типа Адамса для решения вырожденных интегро-дифференциальных уравнений

1.    Постановка задачи и необходимые сведения

При математическом моделировании сложных электрических цепей возникают взаимосвязанные интегро-дифференциальные и интегральные уравнения Вольтерра второго и первого родов [1, 2]. Если данные уравнения объединить, то получим систему интегро-дифференциальных уравнений с тождественно вырожденной матрицей перед производной.

Настоящая работа посвящена численному решению системы

A(t)x (t) +

t

B(t)x(t) + У K (t, s)x(s)ds = f (t), t E [0,1],

с начальным условием

x(0) = x o ,

где A(t),B(t),K(t, s) — (n x п)-матрицы, f (t),x(t) - n-мерные известная и искомая вектор-функции. Элементы A(t), B(t), K (t, s), f (t) предлагаются достаточно гладкие.

В статье рассмотрен случай, когда матрица A(t) тождественно вырожденная, т.е.

det A(t) = 0 V t E [0,1].                                          (3)

Под решением системы (1) мы будем понимать любую вектор x(t) E C 1 (T), которая обращает исходную систему в тождество и удовлетворяется начальным условием (2).

Систему вида (1) с условием (3) будем называть системой интегро-дифференциальных уравнений с тождественно вырожденной матрицей. Задачи (1) – (2) при выполнении условия (3) принципиально отличаются от систем интегро-дифференциальных уравнений, разрешенных относительно производной. В частности, для существования решения системы (1) необходима разрешимость системы линейных алгебраических уравнений вида

A(to)x‘(to) = f (to) - B(to)xo относительно вектора x (to). Для этого необходимо и достаточно выполнения условия: rank(A(to)) = rank(A(to)|f (to) - B(to)xo).

Даже если решение задачи (1) – (2) существует, то оно может быть неединственным. Ниже мы приведем достаточные условия существования единственного достаточно гладкого решения рассматриваемой задачи и обоснуем многошаговые методы ее численного решения.

Приведем класиффикацию рассматриваемых задач:

• •    Если K (t,T) тождественно нулевая матрица, то такие системы принято называть дифференциально-алгебраическими уравнениями (ДАУ);

• •    Если A(t) тождественно нулевая матрица, detB(t) = 0 B(t) = 0, то мы будем иметь интегро-алгебраическое уравнение (ИАУ);

• •    Если A(t) и B(t) тождественно нулевые матрицы, то мы будем иметь систему интегральных уравнений Вольтерра первого рода (СИУВ).

Степень сложности систем вида (1) определяется с помощью характеристики, называемой индексом. По аналогии с [3], будем пользоваться следующим определением индекса:

Определение 1. Минимальное число r, при котором существует квазилинейный дифференциальный оператор r fir = EW (t)( dt )j -j=o где Wj (t) — (n x n) - матрицы из C(t), со свойствами

t fir о (A(t)x (t) +

B (t)x(t) + j o

K(t, s)x(s)ds)

t причем

~ z ,     / z s

A(t)x (t) +

B(t)x(t) + j o

K(t, s)x(s)ds = fi r о f (t),

A(t) = 0, Vt e [0,1], назовем индексом системы (1).

Приведем необходимые определения и вспомогательные результаты для дальнейших исследований:

Определение 2. [11] Матричный полином AA(t) + ^B(t) + C (t), t e [0,1], где A(t), B (t), C (t) — (n x n) матрицы, A, ^ - скалярные параметры, имеет простую структуру, если выполнены условия:

• 1)    rankA(t) = const = r, Vt e [0,1];

• 2)    rank[A(t) | B(t)] = r + l = const V t e [0,1];

• 3)    det[AA(t) + ^B(t) + C(t)] = a o (t)A ^r ^ l + ..., причем a o (t) = 0 V t e [0,1].

Теорема 1. [3] Пусть для задачи (1), (2) выполнены условия:

• 1)    Элементы A(t),B(t),f (t),K(t, s) принадлежат классу C p+'j ;

• 2)    Матричный полином AA(t) + ^B(t) + K(t,t) имеет простую структуру;

• 3)    rankA(0) = rank(A(0) | f(0) — B (0)x(0)).

Тогда задача (1), (2) имеет на отрезке [0,1] единственное решение x(t) из класса 6 p 1] .

Лемма 1. [6] Пусть r-мерные векторы θi вычисляются по правилу i-1

9 i = 6 9 — + h ^ D i,i 9 i + F i ,i = 2,3,..., N, h = 1/N,                  (4)

l=1

где векторы 9q,9i, ...9 k - i - заданы, ||D i,i || <  L i <  to, ||Fi H <  L 2 <  to и собственные значения матрицы 6 удовлетворяют условию | A j | <  1, и на границе единичного круга нет кратных собственных чисел.

Тогда справедлива оценка

||9i|| < L3M, i = 2,3,...,N,L3 < to, где |Фi| = max{|Fi|, |9i|}.

Лемма 2. [11] Если матричный полином AA(t) + ^B(t) + 6(t) имеет простую структуру на всем отрезке [0,1], и элементы A(t),B(t),6(t) принадлежат классу функций 6 p0 1] , то существуют невырожденные матрицы P(t) и Q(t) с элементами из C p 1] такие, что

{E r 0 0\                    /B i (t)

0 B 3 (t)

E l    0

P (t)A(t)Q(t) = 0 0 0 ,P (t)B(t)Q(t) =0

\ 0 0 0/\ 0

/61 (t) 62 (t)0

P (t)6Q(t) = 64(f) 65(f)     0,

\ 0     0   En-r-J где Em – единичная матрица размера m.

Из леммы 2 следует

Лемма 3. Если матричный полином AA(t) + ^B(t) + 6 (t) имеет простую структуру на всем отрезке [0,1], и матрицы A(t), B(t) и C(t) имеют блочный вид

A ¹ (t)         A 1 (t) A 2 (t) A 3 (t)                 B ¹ (t)         B 1 (t) B 2 (t) B 3 (t)

A(t) = I   0   I = I   0       0       0   I , B(t) = I B ² (t) I = I B ₄ (t)  B ₅ (t)  B ₆ (t) I ,

\ 0 /   \ 0    0    0 /        \0/\0    0    0 /

/ 6 1 (t) \     / 6 1 (t) 6 2 (t) 6 3 (t) \

6 (t) = 6 2(t) = 64(f) 65(t) 66(t) , V 6 3(t)/    \67(t) 68(t) 69(t)/ где A1(t), B1(t), 6 1(t) - матрица размер (r x n), B2(t), 62(t) - матрица размер (l x n), 63(t) - матрица размер (n—r—l)xn) и rank(A1(t)) = r = const Vt G [0,1], rank((A(t)|B(t)) = r + l = const Vt G [0,1], то

A ¹ (t)

det qB2(t) 1=0, Vt G [0,1], \s6 3(t) / где q = 0, s = 0.

Вначале приведем описание численных методов решения интегральных уравнений Вольтера I рода, основанных на явных квадратурных формулах типа Адамса.

Зададим на отрезке [0,1] равномерную сетку t i = ih,i = 1, 2, ...,N, h = 1/N, и пусть известны значения y i ~ x(t i ). Тогда

• t i+1                t k + 1                   i t j + 1

  / У(т)dT = / у(т)dT + 52 / У(т)dT

J              J             j=k+1 tj о                о

^t k+1

j Lk +i ⁽ y o ,y i , о

i

.,y k ^,т w ⁺ 52

j=k+1

t j+1

j ^L k₊₁ ⁽ y j - k ^,y j — k+1 , ...,y j ^,TW = t j

k    iki

=h    в1 У1 +      h    Yiyj-1 = h E ^i+1,iyii l=0        j=k+1 l=0l=0

где L^lk+₁ (y _{i -} k ,y _{i -} k+₁ , ...,y _i ,t') — интерполяционный полином степени k, проходящей через точки (y i — k ,t i — k ), (y i - k+1 ,t i - k+1 ), ..., (y i ,t i ). Выпишем коэффициенты Y l (см., например [6]) в таблице 1.

Таблица 1

Коэффициенты γ _l

k	Y o	Y 1	Y 2	Y 3	Y 4	Y 5	Общий множитель
1	1	–	–	–	–	–	1
2	3	–1	–	–	–	–	1/2
3	23	–16	5	–	–	–	1/12
4	55	–59	37	–9	-	–	1/24
5	1901	–2774	2616	–1274	251	–	1/720
6	4277	–7923	9982	–7298	–2877	–475	1/1440

Приведем коэффициенты ^ i+1,i для к = 0,1, 2, 3. (см. например [6—8]):

1		4
1	1	33
1	11	1     3    2    3
ш i+1,l =    1	111     ’ Ш i+1,l	= 2    3   2   2   3
1	1111	3    2    2    2	3
...	... ... ... ...	... ... ... ...	...
	/ 9   0	27                   \
	9   5	11 23
	1    9    5	16   7 23
	^ i +1 , l = 12   9    5	16 12   7 23	,
	9   5	16 12 12   7 23
	... ...

1 ^ i+1, ¹

/ 64	- 32	64
55	5	5	55
55	- 4	42	- 4	55
55	- 4	33	33	- 4	55
55	- 4	33	24	24	- 4
55	- 4	33	24	24	33

- 4

V

\

Отметим, что для нечетных k коэффициент W i+1,g = 0, поэтому нам не потребуется начальное значение x g (см. [6—8]) и в этих таблицах при нечетных к первый нулевой столбец опущен.

Отметим, что из самих формул приближенного вычисления интеграла (5) следует ре- куррентное соотношение:

Г W i+1j = w ij ,j = 0,1, ...,i - к - 1,

^W i+1,i - k = W i^i - k + ^Y k , W i+i^i - k - 1 = W i^i - k - 1 + ^Y k — 1 , ..., ^W i+1,i = ^W g .

В работах [6, 8] были построены многошаговые методы, которые основаны на формуле (5) для численного решения интегрального уравнения

t j K(t, s)x(s)ds = f (t), t E [0,1], 0 < s < t,                          (7)

g с достаточно гладкими ядром K(t, s) и f (t) с условиями

K (t, t) = 0 V t E [0,1], f (0) = 0.

Данные k-шаговые методы имеют вид

i h^Wi+UKi+1,ixi = fi+1,i = k, k + 1,..., N - 1,                     (8)

i=g где Ki+1,1 = K(ti+1, ti), fi+1 = f (ti+1),xi ~ x(ti), h = 1/N, w называется весами квадратурной формулы.

В статье [10] были выделены задачи вида:

A(t)x (t) +

t

B (t)x(t) + j g

K(t,s,x(s))ds = f(t), t E [0,1],

где detA(t) = 0, с начальным условием

x(0) = xg, имеющие единственное решение. Для численного решения таких задач были предложены и обоснованы многошаговые методы, основанные на формулах дифференцирования назад для первых двух слагаемых и на явных квадратурных формулах Адамса для интегрального слагаемого. Эти методы имеет вид:

ki

Ai+1 У^ ajXi+1—j + hBi+1Xi+1 + h2 ^ Wi+1,iKi+1,ixi = hfi+1, i = k, k + 1, •••, N, (10) j=g                               j=g где ai — коэффициенты формулы дифференцирования назад, а Wi+1,1 - веса квадратурной формулы явного метода Адамса. Для задач, рассмотренных в данной статье, методы (10) принципиально неприменимы в силу того, что матрица aoAi+1 + hBi+1 может быть тождественно вырожденной (детали приведены в заключительном разделе).

2.    Численные методы решения

В данном разделе приведены и обоснованы моногошаговые методы численного решения задачи (1), (2), которые удовлетворяют условиями теоремы 1. Приведем численный метод решения задачи (1), (2) с условием (3). Зададим на отрезке [0, 1] равномерную сетку t i = ih, i = 1, 2,...,N, h = 1/N. Обозначим A i = A(t _i ),K _ij = K (t i ,t j ), f i = f (t _i ),x _i ~ x(t i ). Тогда общие многошаговые методы имеют вид:

k                        ki

Ai+1 £ aj xi-j+1 + hBi+1£ вj Xi+1-j + h2£Wi+1l Ki+1lxl = hfi+1.

j=0                    j=0l

Предполагается, что при реализации данных алгоритмов стартовые значения X 1 ,X 2 , ...,X k — 1 вычислены достаточно точно и x o = x(0).

Формулы (11) могут быть:

• 1)    явными при a o = 0, в o = 0, w _i+₁,_i+₁ = 0;

• 2)    полуявными при a o = 0, в о = 0, v i+1,i₊1 = 0, или при a o = 0, в о = 0, v i+1,i₊1 = 0;

• 3)    неявными при a o = 0, в = 0, w _i+₁,_i+₁ = 0.

В силу вырожденности матрицы A(t) явные методы для рассматриваемых задач применять нельзя. Полуявные методы, также неприменимые для задачи (1), (2), удовлетворяющих условиями теоремы 1 в силу того, что матрица (a o A i+₁ + hв o B _i+₁ ) или (a o A i+₁ + h ² W i+1,i+1 K i+1,i+1 ), могут быть вырожденными. Напомним (см. теорему 1 и лемму 2), что рассматриваемый класс задач включает в себя и интегральные уравнения Вольтера первого рода с ядром на диагонале не равным нулю. А для таких уравнений многие неявные многошаговые методы порождают неустойчивый процесс.

Для численного решения рассматриваемых задач мы предлагаем модифицированные методы, основанные на явных формулах Адамса. А именно, для вычисления слагаемого в уравнении (1) будем использовать к — 1-шаговый явный метод Адамса, описанный в предыдущем параграфе. Вырождение A(t)x (t) + B(t)x(t) в точке t = t i+1 будем находить по экстраполяционным формулам. Будем вычислять x i+₁ как значение интерполяционного полинома степени к — 1, проходящего через (x^t i ), (x _{i - 1} ,t _{i - 1} ),..., (x _{i -} k+₁ , t _{i -} k+₁ ) в точке t = t i+₁ , т.е.

k - 1

x i+1 — L k - 1 ^{( x} i , x i - 1 , ..., x i - k+1,t ) — ^ ' ^в j x i - j ^.

j=o

Аналогично, x i+1 ~ x (t) | t=t _i+1 — значение производной интерполяционного полинома степени к, проходящего через точки (x^t i ), (x _{i - 1} ,t _{i - 1} ),..., (x i - k ,t i - k ) в точке t = t i+ы т.е.

x _i+1 = L _k (x i

^x i - 1 ^,

Коэффициенты β j

k	β 0	β 1	β 2	β 3	β 4	β 5
0	1	–	–	–	–	–
1	2	–1	–	–	–	–
2	3	–3	1	–	–	–
3	4	–6	4	–1	–	–
4	5	–10	10	–5	1	–
5	6	–15	20	–15	6	–1

Таблица 2

Приведем коэффициенты α j для различных k = 1, 2, 3, 4, 5, 6 в табл. 3.

Коэффициенты α j

Таблица 3

k	α 0	α 1	α 2	α 3	α 4	α 5	α 6	Общий множитель
1	1	–1	–	–	–	–	–	1
2	5	–8	3	–	–	–	–	1/2
3	26	–57	42	–11	–	–	–	1/6
4	127	–414	534	–322	75	–	–	1/12
5	522	–1755	2540	–1980	810	–137	–	1/60
6	669	–2637	4745	–4920	3015	–1019	147	1/60

Непосредственные вычисления показывают, что корни характеристических уравнений

k

α j p ^k-j =0                                  (13)

j=0

при к <  6 лежат в единичном круге, и на границе круга нет кратных корней, а при к >  7 есть хотя бы один корень, по модулю больший единицы.

Аналогично, для характеристических уравнений [4]

k β j p ^k-j = 0                                      (14)

j=0

и для уравнений (6) – (8)

k

γ j p ^k-j = 0                                      (15)

j=0

при к <  5 корни лежат в единичном круге, и на границе круга нет кратных корней.

Теорема 2. Пусть для задачи (1 ), (3 ) выполнены условия теоремы 1, тогда справедлива оценка ^ X j — x(t j ) || = O(h ^k ) V j = 1, N, где x j определены из системы (12), ||x i — x(t i ) || <

K l h k , K l <  ro,l = 0,1,..., k — 1 .

Доказательство. В силу леммы 2 и второго условия теоремы, существует невырожденная матрица P(t) с элементами из C [k⁺11 такая, что

/ A i+1 \ P i+1 A i+1 = I О I .

\ 0 / где rankA1+1 = r = const. Здесь и в дальнейшем принято обозначение Pi+1 = P(ti+1). Умножая (12) на матрицу Pi+1, получим

/ A i+Л 1 ^k         / B A ^k           i       / K i +A      [ g i+1 \

0 A 52 ajxi-j + Bi+1  52 вxi-j + h 52 ^i+1,l Ki+1,l  xl = gi+1 .

h

\ 0 /  j=0        \ 0 / j=0         l=0      \K3+1V где

/ A ¹ A ² A³\    /A ⁴ ₊₁\

Pi+1Ai+1 =   0   0   0   =    0

/ B ¹ B ² B ³ \/

Pi+1Bi+1 =  B4  B5  B6  =  B2+1

(tz 1 7^2 т^ЗХ      /1< 1 X

K  K  K       Ki+1,l

K4  K5  K6  =  Ki+1,1  , Pi+1fi+1 =  g2+1

К7 К8 K9     3 К33

K   K   K /     \K i+1,l/                 \^g i+1/

Положим E i = X i — x(t i ), тогда уравнения ошибки имеет вид:

kki

A1+1 52 ajEi-j + hB1+1 52 ejEi-j + h2 52^i+1,lKi+1,lEl = h5i+1, j=0                 j=0

hB?+1 ^T ejE,-, + h2 J4 ^i+1lKi+1,lEl = h^.+1.(18)

j=0

i h2 ^i+1,1 Ki+1,1 El = hpi+1.(19)

l =0

Вычитая из (i + 1)-го уравнения (19) i-е, умножая полученное на h ^{- 2} , обозначая

△K3+1,1 = K3+1,1 — Ki31, ^pi+1 = pi+1 — pi и учитывая формулу (8), получим i-k-1

52 ^ i+1,1 ^^K i3+1,1 ^E 1 ^{+ (( w} i,i - k ^{+ Y} k ^{) K 3} +1,i - k ^{— w} i,i - k ^Ki,i - k )^E i - k ⁺ l =0

⁽⁽ ^ i,i - k+1 ^{+ Y} k - 1 ^{) K 3}+1,i - k+1 ^— ^ i,i - k+1 ^K i³i - k+1 ^{) E} i - k+1 ⁺ ••• ^{+ Y} 0 K ³+1,i ^E i = ^h 1 A^p i+1 .

i = k,k + 1.....N — 1.

Эти равенства перепишем в виде

k

EY j K i+i,i - j £ — +E^ i+1,1 AK     = h ^{- i} Ap i+i , i = k,k + 1,-,N — 1.     (20)

j=0

l=0

По первому условию теоремы (достаточная гладкость входных данных) справедливо представление

AK 3+i,i = hM i+1,1 , HM^ || <  M ³ <  to , i = k,k + 1,...,N - 1.

AK^ = hM i+i,i , ^ M 2+i,i || <  M ² <  to, i = k,k + 1,...,N - 1.

Используя стандартные оценки для приближенного вычисления определенного интеграла [5] и учитывая первое условие теоремы, получим:

H Ap i+i h <  L2h^k^\ i = k,k + 1,...,N - 1,

или

|| h ⁱ Ap _i+i || <  L ₂ h ^k , L ₂ <  to, i = k, k + 1,..., N — 1.

По аналогии имеем:

Pi +i h <  L 3h^k , L 3 <  to, i = k, k + 1,..., N - 1.

Учитывая (20), перепишем уравнения (17),(18),(19) в виде:

^a 0 ^{A i}+i ⁺ hAB +i \     / ^a i A i+i ⁺ h^e i B i+i \

ABA I ^e + ⁺ I ^ei B ²+i        I ^£ i - i ⁺ ••• ⁺

Y 0 K ³+i       J \      YiK ³+i      /

/ ^a _k A i+i ⁺ h^e k B +i \           i       / K ³+i,i \

⁺ I         ^в к B 2+i         I ^£ i - k ^{+ h} ^ ^ i+i,l I   M i+i,l  I ^£ l = F i ,

3                        l=0              3

\       ^{Y k K} i+i       )                       \ Mi+i,l /

где F i = (d i+_i ,^ i+_i , h ⁱ Ap T+_i ), i = k, k + 1,..., N - 1.

В силу второго условия теоремы и леммы 3, при достаточно малых h, матрицы

a0 Ai+i + he0Bii+i\ eiB2+i       I , i = k, k + 1,.., N - 1

^Y 0 ^K i³+1,i

являются невырожденными, а по первому условию теоремы (достаточная гладкость входных данных) справедливо представление

a0Ai+i + he0BiL+i\   / aq Ai+i + heq Bi+ \ e0B2+i       I I      eqB2+i       I = Eq + hDi+ii-q,

Y 0 K ³+i ,i     )    \    Y q ^K i³+i,i - q    )

где

a o E r +M i

Eq = I       0

q

e 0 E l + M 2

γ _q

Y o E n - r - 1+М з

а ||D i+i,i — q || <  d <  го , i = к, к + 1,..., N — 1.

( a o A 1₊₁ + he o B 1+i \ 1

β 0 B _i²₊₁           , получим

γ0Ki3+1,i ki

E i + У? E l ^£ i - i + h^M i+i,j E j = F i , l=1             j=1

где

''  + h -B . <    / <_ц \

M i+1,j = W i+i,j

^β 0 ^B i +1              ^M i +1 ,j  ^,

^γ 0 ^K i ³+1 ,i             ^M i ³+1 ,i

F i =

^a 0 ^A 1+1 ^{+ he} 0 B i+1 ^β 0 ^B i²+1 ^γ 0 ^K i³+1,i

- 1

F i , i = k,k + 1,...,N.

Обозначая 9i = (et, ET-i,..., E^-k)T в формулах (24) стандартным образом, совершим переход от многошаговых методов к одношаговым. Будем иметь i-1

9i = C9i-i + h £ Di,i9i + ^i,(25)

l =1

где

/— E i — E 2 ... —E k \

E     0     ...0

C =    0     E    ...    0,

.               .           ....

0      0     ...E нормы матриц Di,j равномерно ограничены, а векторы Фi определены по правилу

Ф i = (F T , 0,0,..., 0) ^T .

Собственные числа матрицы C совпадают с корнями характеристических уравнений (13), (14), (15), которые, как было отмечено выше, при к <  6 по модулю меньше единицы. Таким образом, в силу условия теоремы

|9o| < L4hk, L4 < го, а в силу оценок (21) и (22), имеем

| ^ i | <  L 5 h k , L 5 <  го , i = к,к + 1,...,N.

Вспоминая, что 9 i = ( e ^t , e i — i ,..., E T - k ) ^T из приведенных выше оценок и леммы 1 следует, что

| E i | < L 6 h k , L 6 <  го , i = к, к + 1,..., N — 1.

3.    Численные эксперименты

Приведем очень простой модельный пример:

1	0	0\   ‘	1	0	1	t x(t) + 0	e t+s	0	0		/ e 2t + te t
0	0	0 1 x (t) +	0	1	0		0	e t - s	0	x ( t )dT =	(1 + t)e t
0	0	0	0	0	0		0	0	e t+2s		te t

,

t где x-трехмерная искомая вектор-функция с начальным условием x(0) = (1,1,1)T. Данный пример удовлетворяет всем условиям теоремы 1, и точное решение x(t) = (e-t,et, e-2t)T. Легко заметить, что методы (10) принципиально нельзя применить для данного примера, т.к. det(aA + hB) = 0.

Усложним данный пример, а именно умножим исходную систему на матрицу

0 1 e ^t

P (t) =   et e2t и произведем замену переменной x(t) = Q(t)y(t), где

Q(t) =

2t

1 0

t ²

3t I .

В итоге получим

t

A(t)y ‘ (t) + B(t)y(t) +

j К ^{( t,}s)^ds = f ^{( t )} , t G [0, 1 ^] , 0

t2                        1

ett2   ,  B(t) =    et e2tt2                     e2t

2t + 2

2te ^t + 1 + 2e ^t 2e 2t t + e t + 2e ^2t

t ² + 1 + 2t e ^t t ² + 3t + e ^t + 2te ^t W + 3te t + e ^{2 t} + 2e ^2t t

/ et + s

K (t, s) = I etet+s e2tet+s

₂ e t+s _t

2etet+st + et-s 2e2t et+st + etet-s et+st2

e t e t + s _t 2 _{+ 3} e t - s _t e ^{2 t} e ^{t + s} t ² + 3e ^t e ^t-s t + e ^{t +2 s}

f (t) =

e ^{- 2 t} + te ^t e ^t (e ^{- 2 t} + te ^t ) + (1 + t)e ^t e ^{2 t} (e ^{- 2 t} + te t ) + (1 + t)(e t ) ² + te t

Точное решение y(t) = Q ¹ (t)x(t), начальное условие y(0) = (1,1,1) T .

Численные расчеты данной задачи при различных значениях k приведены в табл. 4.

Таблица 4

Численные расчеты задачи при различных значениях k

err	k=1	k=2	k=3
N=5	1,309600415814891	0,6015407275019990	0,21171281782986052430
N=10	0,7497289570481798	0,1844243516458794	0,04761740960151257878
N=20	0,3988507964835724	0,0503707677718254	0,00732509005266374868
N=40	0,2051764163549656	0,0129986398315527	0,00097017989140169301
N=80	0,1039752161311108	0,0032742356352037	0,00012382133627371258

Здесь err = max k ≤ i ≤ N

- y 1 (t i )) 2 + (y 2 - y ² (t i )) 2 + (y ³ — У 3 (t i )) ² .

В качестве стартовых значений для y 1 , y 2 , ..., y k _- 1 были взяты точные значения y(t i ), y(t 2 ),..., y(t k - i ). Из приведенной табл. 4 видно, что расчеты хорошо согласуются с теоритическими выкладками.

Исследования были поддержаны грантом РФФИ 13-01-93002Вьет-а.