The evaluation of the order of approximation of the matrix method for numerical integration of the boundary value problems for systems of linear non-homogeneous ordinary differential equations of the second order with variable coefficients. Message 1. Boundary value problems with boundary conditions of the first kind

Full Text

Введение. Использование конечных разностей для аппроксимации производных в классическом методе сеток численного интегрирования краевых задач для неоднородных обыкновенных дифференциальных уравнений второго порядка (ОДУ2) с переменными коэффициентами с граничными условиями первого рода x + p(t)x + q(t)x = f (t), x(a) = x0 , x(b) = xn , (1) где p(t), q(t), f (t) - заданные функции, дифференцируемые нужное число раз, [a, b] - область интегрирования, x0 , xn - заданные числа, приводит ко второму порядку аппроксимации [1-6]. Такой же (второй) порядок аппроксимации имеют методы численного интегрирования ряда краевых задач для уравнений в частных производных [4-10]. Второй порядок обусловлен тем, что при аппроксимации производных конечными разностями было удержано всего три члена разложения в ряд Тейлора искомого решения задачи. Метод, использующий средства матричного исчисления, численного интегрирования разностных краевых задач для неоднородных ОДУ2, позволяющий увеличить количество удержанных членов до произвольного натурального числа в разложении искомой функции в ряд Тейлора, предложен в [11], при этом аппроксимация производных конечными разностями не использовалась. Оценка порядка аппроксимации предложенного в [11] матричного метода интегрирования краевых задач для неоднородных линейных ОДУ2 с граничными условиями первого, второго и третьего рода при различных значениях удержанного числа членов в рядах Тейлора дана в [12]. Исследуем возможность использования матричного метода для численного интегрирования систем неоднородных линейных ОДУ2 с переменными коэффициентами с различными граничными условиями на примере следующей системы для неизвестных функций x(t), y(t): u1 x + p1 x + r1 x + v1 y + q1 y + s1 y = f1 , u2 x + p2 x + r2 x + v2 y + q2 y + s2 y = f2 , (2) где uj , pj , rj , vj , qj , sj , fj - заданные функции аргумента t, дифференцируемые нужное число раз, j = 1, 2 - номер уравнения в системе (2), и поставим целью вычисление порядка аппроксимации. Далее будем придерживаться принятых в [4] обозначений: 1) D - область интегрирования, ограниченная отрезком [a, b], Dh - узлы сетки, определяемые значениями ti = t0 + ih, i = 1, 2, . . . , n, t0 = a, tn = b, h = (b - a)/n, n + 1 - число узлов сетки; 2) x(t), y(t) - непрерывные функции, являющиеся точным решением системы (2) с теми или иными граничными условиями; 390 Оценка порядка аппроксимации матричного метода численного интегрирования краевых . . . 3) [x]h , [y]h - сеточные функции, совпадающие с точным решением в узлах сетки Dh ; 4) x(h) , y (h) - искомые сеточные функции. Для краткости примем для любой функции обозначение ϕ(ti ) = ϕi , где ti - узел сетки Dh . В дальнейшем опустим индекс h в наименованиях сеточных функций [x]h , [y]h , x(h) , y (h) . 1. Краевые задачи с граничными условиями первого рода. Рассмотрим систему ОДУ2 (2) с граничными условиями первого рода x0 = x0 , y0 = y0 , xn = xn , yn = yn , (3) где x0 , y0 , xn , yn - заданные числа. В соответствии с матричным методом численного интегрирования выполним следующие преобразования. Для некоторых фиксированных k 2, k - степень используемого многочлена Тейлора, и номера узла сетки i, i = 1, 2, . . . , n - 1, запишем следующие многочлены Тейлора для функций x(t), y(t):  k h3 h2  k h x(k) ,  x - x + · · · + (-1) x = x - hx +  i-1 i i   2! i 3! i k! i   2 3  hk h h   yi-1 = yi - hyi + yi - yi + · · · + (-1)k yi(k) , 2! 3! k! (4)  h2 h3 hk (k)   xi+1 = xi + hxi + xi + xi + · · · + xi ,   2! 3! k!    2 3 k  h h h  yi+1 = yi + hy + y + y + · · · + y (k) , i i i 2! 3! k! i вычислим производные по аргументу t от обеих частей уравнений системы (2) и запишем их в узле ti в виде  (u1i x + p1i x + r1i xi + v1i y + q1i y + s1i yi )(r) = (f1i )(r) , i i i i (5) (u x + p x + r x + v y + q y + s y )(r) = (f )(r) , 2i i 2i i 2i i 2i i 2i i 2i i 2i где r = 1, 2, . . . , k - 2; здесь и далее в наименованиях функций первый из пары нижних индексов означает номер уравнения в системе (2), второй - номер узла сетки. Из многочленов Тейлора (4), дифференциальных уравнений системы (2), записанных в узле ti , и производных (5) составим следующую систему из 2k + 2 уравнений:  k h2 h3  k h x(k) = x  x - hx + x - x + · · · + (-1)  i i-1 , i   2! i 3! i k! i   k  h3 h2  k h y (k) = y   y - y + · · · + (-1) y - hy + i i-1 , i   2! i 3! i k! i    h2 h3 hk (k)    xi + hxi + xi + xi + · · · + x = xi+1 , (6)   2! 3! k! i    h2 h3 hk (k)    y + y + · · · + y = yi+1 , y + hy + i  i i i  2! 3! k! i    r1i xi + p1i xi + u1i xi + s1i yi + q1i yi + v1i yi = f1i ,     r2i xi + p2i xi + u2i xi + s2i yi + q2i yi + v2i yi = f2i , 391 М а к л а к о в В. Н.  r1i xi + (r1i + p1i )xi + (p1i + u1i )xi + u1i xi +    +s1i yi + (s1i + q1i )yi + (q1i + v1i )yi + v1i yi = f1i ,       r2i xi + (r2i + p2i )xi + (p2i + u2i )xi + u2i xi +    +s2i yi + (s2i + q2i )yi + (q2i + v2i )yi + v2i yi = f2i ,      ..    .    (k-2) (k) (k-2) (k) (k-2)   r1i xi + · · · + u1i xi + s1i yi + · · · + v1i yi = f1i ,     (k-2) (k) (k-2) (k) (k-2) r2i xi + · · · + u2i xi + s2i yi + · · · + v2i yi = f2i . Система (6) - система линейных алгебраических уравнений (СЛАУ), имеющая в матричной форме вид Aki W ki = Gki , где  Aki 1 0 1 0 r1i r2i ... 0 1 0 1 s1i s2i ...        =      (k-2) (k-2) r1i s1i (k-2) (k-2) r2i s2i 2 h 0 -h 0 2! h2 0 -h 0 2! h2 h 0 0 2! 2 0 h 0 h2! p1i q1i u1i v1i p2i q2i u2i v2i ... ... ... ... ... ... ... ... ... ... ... ... k . . . (-1)k hk! ... 0 hk ... k! ... 0 ... 0 ... 0 ... ... ... u1i ... u2i (k) W ki = xi yi xi yi xi yi . . . xi Gki = xi-1 yi-1 xi+1 yi+1 f1i f2i . . . (k) yi (k-2) f1i  0 k (-1)k hk!     0  hk  k!  , 0   0  ...   v1i  v2i (7) , (k-2) f2i . Здесь и далее верхний индекс k означает степень используемого многочлена Тейлора, если речь не идет о показателях алгебраических степеней и степенях производных; второй из пары верхних индексов (если он существует) в наименованиях матриц и их элементов означает номер узла сетки. Матрицу Aki будем называть локальной матрицей. Предполагая существование обратной матрицы B ki = (Aki )-1 от локальной матрицы (7), запишем матричное равенство B ki Gki = W ki , или следующие соотношения в координатной форме: ki ki ki ki ki bki 11 xi-1 + b12 yi-1 + b13 xi+1 + b14 yi+1 + b15 f1i + b16 f2i + · · · + (k-2) + bki 1 (2k+1) f1i (k-2) = xi , (8) (k-2) = yi , (9) + bki 1 (2k+2) f2i ki ki ki ki ki bki 21 xi-1 + b22 yi-1 + b23 xi+1 + b24 yi+1 + b25 f1i + b26 f2i + · · · + (k-2) + bki 2 (2k+1) f1i 392 + bki 2 (2k+2) f2i Оценка порядка аппроксимации матричного метода численного интегрирования краевых . . . ki ki ki ki ki bki 31 xi-1 + b32 yi-1 + b33 xi+1 + b34 yi+1 + b35 f1i + b36 f2i + · · · + (k-2) + bki 3 (2k+1) f1i (k-2) = xi , (10) (k-2) = yi , (11) + bki 3 (2k+2) f2i ki ki ki ki ki bki 41 xi-1 + b42 yi-1 + b43 xi+1 + b44 yi+1 + b45 f1i + b46 f2i + · · · + (k-2) + bki 4 (2k+1) f1i + bki 4 (2k+2) f2i ... ki ki ki ki bki (2k+1) 1 xi-1 + b(2k+1) 2 yi-1 + b(2k+1) 3 xi+1 + b(2k+1) 4 yi+1 + b(2k+1) 5 f1i + (k-2) ki + bki (2k+1) 6 f2i + · · · + b(2k+1) (2k+1) f1i (k-2) + bki (2k+1) (2k+2) f2i (k) = xi , (12) ki ki ki ki bki (2k+2) 1 xi-1 + b(2k+2) 2 yi-1 + b(2k+2) 3 xi+1 + b(2k+2) 4 yi+1 + b(2k+2) 5 f1i + (k-2) ki + bki (2k+2) 6 f2i + · · · + b(2k+2) (2k+1) f1i (k-2) + bki (2k+2) (2k+2) f2i (k) = yi , (13) ki где bki jm - элементы матрицы B . Из соотношений (8), (9), являющихся разностными уравнениями второго порядка для трехточечного шаблона ti-1 , ti , ti+1 , i = 1, 2, . . . , n - 1, и связывающих значения xi-1 , yi-1 , xi , yi , xi+1 , yi+1 искомых функций, составим, учитывая граничные условия (3), систему  k1 x - bk1 y = bk1 x + bk1 y + x1 - b13  2 14 2 11 0 12 0   (k-2) (k-2)  k1 f + bk1 f + · · · + bk1  +b + bk1 ,  15 11 16 21 1 (2k+1) f11 1 (2k+2) f21    k1 k1 k1   y1 - bk1 23 x2 - b24 y2 = b21 x0 + b22 y0 +    (k-2) (k-2) k1 k1  +bk1 + bk1 ,  25 f11 + b26 f21 + · · · + b2 (2k+1) f11  2 (2k+2) f21   .   ..     ki ki ki  -bki  11 xi-1 - b12 yi-1 + xi - b13 xi+1 - b14 yi+1 =   (k-2) (k-2)  ki ki  = bki + bki , 15 f1i + b16 f2i + · · · + b1 (2k+1) f1i 1 (2k+2) f2i (14) ki ki ki ki  -b21 xi-1 - b22 yi-1 + yi - b23 xi+1 - b24 yi+1 =    (k-2) (k-2) ki ki  = bki + bki ,  25 f1i + b26 f2i + · · · + b2 (2k+1) f1i  2 (2k+2) f2i   .   ..     k (n-1) k (n-1) k (n-1) k (n-1)   -b11 xn-2 - b12 yn-2 + xn-1 = b13 xn + b14 yn +    k (n-1) k (n-1) k (n-1) (k-2) k (n-1) (k-2)   +b15 f1 (n-1) + b16 f2 (n-1) + · · · + b1 (2k+1) f1 (n-1) + b1 (2k+2) f2 (n-1) ,    k (n-1) k (n-1) k (n-1) k (n-1)   -b21 xn-2 - b22 yn-2 + yn-1 = b23 xn + b24 yn +    k (n-1) k (n-1) (k-2) k (n-1) (k-2)  +bk (n-1) f +b f + ··· + b f +b f . 25 1 (n-1) 26 2 (n-1) 2 (2k+1) 1 (n-1) 2 (2k+2) 2 (n-1) 393 М а к л а к о в В. Н. Отметим, что в центральном узле шаблона уравнения системы (14) с нечетными номерами содержат искомые значения xi , с четными номерами - значения yi . Система (14) является СЛАУ, имеющей блочную трехдиагональную матрицу, решение которой позволяет найти метод матричной прогонки [5, 13], что дает возможность вычислить искомые значения xi , yi во внутренних узлах (r) (r) сетки ti , i = 1, 2, . . . , n - 1. Производные xi , yi , r = 1, 2, . . . , k, во внутренних узлах сетки найдем по формулам (10)-(13) с использованием найденных значений xi , yi . Фактически СЛАУ (14) есть запись в развернутом виде разностной краевой задачи для системы двух разностных уравнений второго порядка с граничными условиями первого рода для отыскания вектор-функции zi = [xi yi ] , i = 1, 2, . . . , n - 1, удовлетворяющей следующим равенствам [4, 14]: 1 Bjki zi+j = g ki , i = 1, 2, . . . , n - 1, λ0 z0 = z0 , λn zn = zn , (15) j=-1 где ki B-1 = g ki = ki -bki 11 -b12 , ki -b21 -bki 22 g1i = C ki F ki , g2i 1 0 , 0 1 B0ki = λ0 = λn = B1ki = ki -bki 13 -b14 , ki -b23 -bki 24 ki ki ki ki ki bki 15 b16 b17 b18 . . . b1 (2k+1) b1 (2k+2) ki ki ki i ki bki 25 b26 b27 b28 . . . b2 (2k+1) b2 (2k+2) C ki = (k-2) F ki = f1i f2i f1i f2i . . . f1i (k-2) f2i , z0 = [x0 y0 ] , , zn = [xn yn ] . 2. Предварительные оценки. При некоторых фиксированных номере k 2 и номере i, i = 1, 2, . . . , n - 1, запишем следующие ряды Тейлора для функции x(t): xi-1 = xi - hxi + h3 h2 xi - x + ··· = 2! 3! i k = mh (-1) m=0 xi+1 = xi + hxi + m m! ∞ (m) xi (-1)m + m=k+1 hm (m) k k x = Px,i-1 + Rx,i-1 , m! i h2 h3 xi + x + ··· = 2! 3! i k = m=0 hm (m) x + m! i ∞ m=k+1 hm (m) k k x = Px,i+1 + Rx,i+1 , m! i где k k Rx,i-1 , Rx,i+1 = 394 hk+1 (k+1) x (ξi ) = O(hk+1 ), (k + 1)! ξi ∈ (ti , ti+1 ) Оценка порядка аппроксимации матричного метода численного интегрирования краевых . . . - дополнительные члены разложений в ряд Тейлора в форме Лагранжа [15]. Из двух последних точных равенств найдем ∞ k k Rx,i+1 - Rx,i-1 = m=k+1 hm (m) x - m! i ∞ (-1)m m=k+1 hm (m) x = m! i ∞ = m=k+1 hm (m) x (1 - (-1)m ) (16) m! i и, аналогично, ∞ k k Rx,i+1 + Rx,i-1 = m=k+1 hm (m) x + m! i ∞ (-1)m m=k+1 hm (m) x = m! i ∞ = m=k+1 hm (m) x (1 + (-1)m ) . (17) m! i Из равенств (16), (17) для четного k имеем k k Rx,i+1 - Rx,i-1 = k k Rx,i+1 + Rx,i-1 = hk+1 (k+1) x (1 + 1) + (k + 1)! i hk+2 (k+2) x (1 - 1) + · · · = O(hk+1 ), (18) + (k + 2)! i hk+1 (k+1) x (1 - 1) + (k + 1)! i hk+2 (k+2) + xi (1 + 1) + · · · = O(hk+2 ) (19) (k + 2)! и для нечетного k - k k Rx,i+1 - Rx,i-1 = k k Rx,i+1 + Rx,i-1 = hk+1 (k+1) x (1 - 1) + (k + 1)! i hk+2 (k+2) + xi (1 + 1) + · · · = O(hk+2 ), (20) (k + 2)! hk+1 (k+1) x (1 + 1) + (k + 1)! i hk+2 (k+2) + x (1 - 1) + · · · = O(hk+1 ). (21) (k + 2)! i Замена в равенствах (18)-(21) x на y даст аналогичные оценки для доk k полнительных членов Ry,i-1 , Ry,i+1 . 395 М а к л а к о в В. Н. 3. Оценка порядка аппроксимации разностной краевой задачи с граничными условиями первого рода. Соответствующую дифференциальной краевой задаче (1) разностную краевую задачу будем обозначать, по аналогии с [4], в компактной символической форме как Lkh x = fhk . При исследовании дифференциальной краевой задачи (1) для одного уравнения сеточная функция xi , i = 0, 1, . . . , n, являющаяся решением разностной краевой задачи, при подстановке в уравнения этой разностной краевой задачи обратит их в верные равенства. В [4] показано, что подстановка в уравнения задачи сеточной функции [xi ], отличающейся от xi , приведет к некоторому отличию от верных равенств. Эти отличия и характеризует невязка δfhk [4]. Иными словами, подстановка [x] в Lkh x = fhk приведет к Lkh [x] = fhk + δfhk . В работе [4] в качестве оценки величины невязки принята норма δfhk = max |δfhk0 |, |δfhkn |, |δfhki | , i = 1, 2, . . . , n - 1, где первые две компоненты характеризуют меру отличий в граничных узлах сетки Dh , оставшиеся - во внутренних узлах. Согласно [4, 6] разностная краевая задача (14) аппроксимирует дифференциальную краевую задачу (1) на точном решении x, если δfhk → 0 при h → 0. Если при этом имеет место неравенство δfhk Chk , где C > 0, k > 0 - некоторые постоянные, не зависящие от h, то говорят, что имеет место аппроксимация порядка k относительно величины h. Вычислим порядок аппроксимации рассматриваемой краевой задачи. СЛАУ (14), или, что то же самое, разностную краевую задачу (15) запишем, по аналогии с [4], в компактной символической форме как Lkh z = fhk , или, в дальнейшем для краткости, как Lkh , где Lkh z ≡ Lkh 396 x y  ki b   xi-1 - - 11   bki  15     bki  21  - ki xi-1 - b 26 =  x0 ,      y0 ,     xn ,    yn , bki xi bki bki 12 13 14 y + - x - yi+1 , i-1 i+1 ki ki ki b15 b15 b15 bki 15 bki yi bki bki 22 23 24 y + - x - y , i-1 i+1 ki ki ki i+1 bki b b b 26 26 26 26 (22) Оценка порядка аппроксимации матричного метода численного интегрирования краевых . . .   bki  16  f + f +  1i  ki 2i  b  15      bki   25 f1i + f2i + ki fhk = b26   x0 ,     y0 ,      xn ,    yn , bki bki bki bki 1 (2k+1) (k-2) 1 (2k+2) (k-2) 17 18 + , f f f f2i + + · · · + 1i 2i 1i ki ki ki b15 b15 b15 bki 15 bki bki bki bki 2 (2k+2) (k-2) 2 (2k+1) (k-2) 27 28 f1i + ki f2i + · · · + f1i + f2i , ki ki b26 b26 b26 bki 26 для всех i = 1, 2, . . . , n - 1. ki Как будет показано ниже, выбор того или иного знаменателя (bki 15 или b26 ) в двух первых равенствах задачи (22) не является существенным при оценке порядка аппроксимации. Для задачи (22) обобщение приведенного выше определения порядка аппроксимации запишем в форме Lkh [z] = fhk + δfhk , (23) где k k δfhk = δf1h δf2h ; k , δf k - величины невязок в узлах сетки D , появление которых обусловδf1h h 2h лено соответственно первым и вторым уравнениями задачи (22). Итоговую норму невязку задачи в соответствии с [4] определим как k k δfhk = max δf1h , δf2h . (24) При фиксированных k, i, i = 1, 2, . . . , n - 1, имеем точные равенства  k h2 h3  k h [x(k) ] = [x k  [x ] - h[x ] + [x ] - [x ] + · · · + (-1)  i i-1 ] - Rx,i-1 , i   2! i 3! i k! i    h2 h3 hk   k [yi ] - h[yi ] + [yi ] - [yi ] + · · · + (-1)k [yi(k) ] = [yi-1 ] - Ry,i-1 , 2! 3! k! (25)  h2 h3 hk (k)  k  [xi ] + h[xi ] + [xi ] + [xi ] + · · · + [xi ] = [xi+1 ] - Rx,i+1 ,   2! 3! k!    2 3 k  h h h  (k) k  [yi ] + h[yi ] + [yi ] + [yi ] + · · · + [yi ] = [yi+1 ] - Ry,i+1 . 2! 3! k! Составим СЛАУ, в которую внесем равенства (25), дифференциальные уравнения системы (2), записанные в узле ti , и точные равенства (5). В итоге получим 397 М а к л а к о в В. Н.  k h3 h2  k k h [x(k) ] = [x  [x [x ] - ] + · · · + (-1) [x ] - h[x ] +  i-1 ] - Rx,i-1 , i i i i   2! 3! k! i    h3 h2 hk (k)   k  , [yi ] - h[yi ] + [yi ] - [yi ] + · · · + (-1)k [yi ] = [yi-1 ] - Ry,i-1   2! 3! k!    2 3 k  h h h (k)  k  [x ] + h[x ] + [x ] + [x ] + · · · + [xi ] = [xi+1 ] - Rx,i+1 ,  i i i i  2! 3! k!     h3 hk (k) h2  k   [y [y [yi ] = [yi+1 ] - Ry,i+1 [y ] + h[y ] + ] + · · · + , ] + i  i i i  2! 3! k!     r1i [xi ] + p1i [xi ] + u1i [xi ] + s1i [yi ] + q1i [yi ] + v1i [y ] = f1i ,    r2i [xi ] + p2i [xi ] + u2i [xi ] + s2i [yi ] + q2i [yi ] + v2i [yi ] = f2i ,     r1i [xi ] + (r1i + p1i )[xi ] + (p1i + u1i )[xi ] + u1i [xi ]+      +s1i [yi ] + (s1i + q1i )[yi ] + (q1i + v1i )[yi ] + v1i [yi ] = f1i ,       r2i [xi ] + (r2i + p2i )[xi ] + (p2i + u2i )[xi ] + u2i [xi ]+     +s2i [yi ] + (s2i + q2i )[yi ] + (q2i + v2i )[yi ] + v2i [yi ] = f2i ,       ..   .    (k-2) (k) (k-2) (k) (k-2)   r1i [xi ] + · · · + u1i [xi ] + s1i [yi ] + · · · + v1i [yi ] = f1i ,      (k-2) (k) (k-2) (k) (k-2) r2i [xi ] + · · · + u2i [xi ] + s2i [yi ] + · · · + v2i [yi ] = f2i . (26) В матричной форме система (26) имеет вид Aki [W ki ] = [Gki ] , где  [xi ]  [yi ]   [x ]   i   [y ]   i   [x ]   i    [W ki ] =  [yi ]  ,  [xi ]     [yi ]     ...   (k)  [xi ]  (k) [yi ]   k [xi-1 ] - Rx,i-1 k   [yi-1 ] - Ry,i-1    [xi+1 ] - Rk x,i+1     [y ] - Rk  i+1 y,i+1    f1i   ki , [G ] =  f2i     f 1i     f 2i     ...   (k-2)   f 1i (k-2) f2i а матрица Aki определяется формулой (7). В предположении существования обратной матрицы B ki = (Aki )-1 от матрицы Aki найдем B ki [Gki ] = [W ki ]. Выпишем первые два уравнения последнего матричного равенства: 398 Оценка порядка аппроксимации матричного метода численного интегрирования краевых . . . k ki k ki k bki 11 ([xi-1 ] - Rx,i-1 ) + b12 ([yi-1 ] - Ry,i-1 ) + b13 ([xi+1 ] - Rx,i+1 )+ k ki ki ki ki + bki 14 ([yi+1 ] - Ry,i+1 ) + b15 f1i + b16 f2i + b17 f1i + b18 f2i + (k-2) + · · · + bki 1 (2k+1) f1i (k-2) + bki 1 (2k+2) f2i = [xi ], k ki k ki k bki 21 ([xi-1 ] - Rx,i-1 ) + b22 ([yi-1 ] - Ry,i-1 ) + b23 ([xi+1 ] - Rx,i+1 )+ k ki ki ki ki + bki 24 ([yi+1 ] - Ry,i+1 ) + b25 f1i + b26 f2i + b27 f1i + b28 f2i + (k-2) + · · · + bki 2 (2k+1) f1i (k-2) + bki 2 (2k+2) f2i = [yi ], ki в узле t , или, после преобразований, где bki i jm - элементы матрицы B - bki [xi ] bki bki bki 12 13 14 11 [x ] - [y ] + - [x ] - [yi+1 ] = f1i + i-1 i-1 i+1 ki ki ki ki b15 b15 b15 b15 bki 15 + bki bki bki bki bki 1 (2k+1) (k-2) 1 (2k+2) (k-2) 16 17 18 f + f + f + · · · + f + f2i - 2i 1i ki 1i ki 2i ki bki b b b bki 15 15 15 15 15 - - k ki k bki 11 Rx,i-1 + b13 Rx,i+1 bki 15 - k ki k bki 12 Ry,i-1 + b14 Ry,i+1 bki 15 , (27) bki [yi ] bki bki bki bki 22 23 24 25 21 [x ] - [y ] + - [x ] - [y ] = f1i + i-1 i-1 i+1 i+1 ki ki ki ki ki ki b26 b26 b26 b26 b26 b26 + f2i + bki bki bki bki 2 (2k+1) (k-2) 2 (2k+2) (k-2) 27 28 f + f + · · · + f + f2i - 1i 1i ki 2i ki bki b b bki 26 26 26 26 - k ki k bki 21 Rx,i-1 + b23 Rx,i+1 bki 26 - k ki k bki 22 Ry,i-1 + b24 Ry,i+1 bki 26 . (28) Отбрасывание двух последних дробей в составленной из уравнений (27), (28) во внутренних узлах сетки Dh системе, что равносильно переходу от точного решения [xi ], [yi ] к искомому приближенному xi , yi , приведет эту систему к разностной краевой задаче (22). Следовательно, в соответствии с (23) две последние дроби в уравнениях (27), (28) характеризуют величину невязки в узлах ti , i = 1, 2, . . . , n - 1; в итоге для рассматриваемой задачи имеем k ki k k ki k bki bki 11 Rx,i-1 + b13 Rx,i+1 12 Ry,i-1 + b14 Ry,i+1 ki δf1h =- - , (29) bki bki 15 15 ki δf2h =- k ki k bki 21 Rx,i-1 + b23 Rx,i+1 bki 26 - k ki k bki 22 Ry,i-1 + b24 Ry,i+1 bki 26 . (30) Первое слагаемое в равенствах (29), (30) характеризует величину невязки, появление которой обусловлено функцией x, второе - функцией y. Для краевой задачи с граничными условиями первого рода величины невязок в граничных узлах сетки обращаются в нуль в силу того, что два 399 М а к л а к о в В. Н. первых уравнения задачи (22) при i = 1, n - 1 содержат заданные значения x0 , y0 , xn , yn искомых функций. ki и δf ki во внутренних узлах сетки D . Исследуем невязки δf1h h 2h В узле ti введем компактные обозначения для определителей второго порядка, например u v (31) Ui Vi = 1i 1i , u2i v2i где, очевидно, Ui Vi = -Vi Ui . В дальнейшем будем опускать индекс i в обозначениях определителей вида (31). Непосредственными вычислениями можно для любого k 3 убедиться в справедливости оценки ki 2i M11 ≈ (U V )k-2 M11 , (32) ki - алгебраическое дополнение элемента aki транспонированной логде M11 11 кальной матрицы Aki . Отсутствие в настоящей работе полного вывода формулы (32) обусловлено лишь громоздкостью выкладок. Для краевых задач для одного ОДУ2 полный вывод формулы (32) дан в [12]. Тем не менее рассмотрим частный случай при k = 3, опуская индекс i - номер узла сетки. Воспользуемся в дальнейшем известными свойствами определителя [16]. Имеем, пренебрегая старшими степенями, следующее: 3 M11 1 0 -h = 02 h 2 0 h 0 1 0 h 0 h2 2 0 h3 3! 0 3 - h3! 0 1 0 0 h -h 0 + u1 0 h2 2 2 h 0 23 - h3! 0 s1 s2 s1 s2 p1 p2 p1 p2 q1 q2 q1 q2 u1 u2 u1 u2 = h2 v1 v2 v1 v2 2 0 0 0 u1 u2 h3 0 0 v1 v2 3! 1 s1 s2 s2 0 p1 p2 p2 h q1 q2 q2 0 u1 u2 u2 - u2 h2 v1 v2 v2 2 h3 0 0 v2 3! 1 0 3 -h h 3! 02 h 23 - h3! 1 0 -h 0 h2 23 - h3! 0 h 0 h2 2 0 0 1 0 h 0 h2 2 h3 3! 1 0 h 0 h2 2 h3 3! s1 s2 s1 s2 p1 p2 p1 p2 q1 q2 q1 q2 u1 u2 u1 u2 + v1 v2 v1 v2 0 0 v1 v2 s1 s2 s1 p1 p2 p1 q1 q2 q1 u1 u2 u1 = v1 v2 v1 0 0 v1 2 = c3 h6 + O h7 + u1 d3 h4 + O h5 + v2 M11 - 2 2 2 - u2 g3 h4 + O h5 + v1 M11 ≈ (u1 v2 - u2 v1 ) M11 = U V M11 , где pj , qj , uj , vj , j = 1, 2, - функции от pj , qj , uj , vj и их первых производных, c3 , d3 , g3 - не зависящие от h величины. Формулы, аналогичные (32), имеют место, по крайней мере, для первых шести элементов первой и второй строк матрицы (Aki ) ; на основании чего и очевидных равенств bki 1j bki 15 400 = k M1j k M15 , bki 2j bki 26 = k M2j k M26 , j = 1, 2, 3, 4, Оценка порядка аппроксимации матричного метода численного интегрирования краевых . . . следуют оценки bki 1j bki 15 ≈ 2 M1j 2 M15 , bki 2j bki 26 ≈ 2 M2j 2 M26 , j = 1, 2, 3, 4. (33) Невязка (29) с учетом соотношений (33) примет вид ki δf1h ≈- 2 Rk 2 k M11 x,i-1 + M13 Rx,i+1 2 M15 - 2 Rk 2 k M12 y,i-1 + M14 Ry,i+1 2 M15 . (34) Вычислим оценки первых пяти алгебраических дополнений первой строки матрицы A2i , i = 1, 2, . . . , n - 1. Имеем, пренебрегая старшими степенями, следующее: 2 M11 1 0 = -h 0 0 h 0 1 0 h 0 h2 2 h2 h2 2 0 2 s1 p1 q1 u1 v1  s2 0 p2 1 q2 = -2 h2 h 2 u2 0 v2 1 2 = -2 h - 0 h2 2 1 0 0 h2 2 s1 s2 p1 p2 u1 u2 = v1 v2 1 s1 s2 h 0 u1 u2 + 2 h2 v1 v2 2  s1 s2 p1 p2  = v1 v2 h h2 h2 SU + PV + SP = 2 2 2 h h2 h3 U V - P V + SU - SP ≈ h2 (2 U V - h P V ) , (35) 2 2 4 = -2 h2 -U V - = 2 h2 где были использованы обозначения определителей, аналогичные введенным формулой (31), 2 M13 0 1 -h 0 = 02 -h h 0 2 h2 0 2 1 0 h 0 s1 s2 p1 p2 q 1 q2 = u1 u2 v1 v2 h2 2 = 2 h2 U V + 2 M12 0 -h = - 02 0 h 0 h 2 h2 2 0 0 1 0 h 0 h2 2 h h2 h3 P V + SU + SP 2 2 4 s1 s2 p1 p2 q1 q2 = -h3 u1 u2 v1 v2 0 0 1 0 1 h 0 h2 2 ≈ h2 (2 U V + h P V ) , (36) s1 s2 q1 q2 u1 u2 = v1 v2 401 М а к л а к о в В. Н. = -h 3 1 h h2 2 s1 s2 h2 q1 q2 = h3 -QV + h SV - SQ ≈ h3 (-QV + h SV ) , (37) 2 v1 v2 2 M14 2 M15 0 1 -h 0 = 02 -h h 0 2 h2 0 2 0 1 -h 0 = - 02 -h h 0 2 h2 0 2 0 h 0 h2 2 0 1 0 h 0 h2 2 0 h 0 h2 2 0 s1 s2 p 1 p2 q1 q2 ≈ h3 (QV + h SV ) , u1 u2 v1 v2 s2 p2 q2 = -4 h3 u2 v2 0 = 2 h4 1 0 1 0 h2 2 0 0 1 0 0 h2 0 0 1 h 0 h2 2 (38) s2 q2 u2 = v2 s2 h2 u2 = 2 h4 -v2 + s2 2 v2 ≈ -2 h4 v2 . (39) С учетом оценок (35)-(39) величину невязки (34), появление которой обусловлено первым уравнением задачи (22), на точном решении [xi ], [yi ] во внутренних узлах ti сетки Dh задачи Lkh запишем как ki δf1h ≈- k k h2 (2 U V - h P V ) Rx,i-1 + (2 U V + h P V ) Rx,i+1 - -2 h4 v2 k k h3 (-QV + h SV ) Ry,i-1 + (QV + h SV ) Ry,i+1 - = -2 h4 v2 k k k k 2 U V (Rx,i+1 + Rx,i-1 ) + h P V (Rx,i+1 - Rx,i-1 ) = + 2 h2 v2 k k k k QV (Ry,i+1 - Ry,i-1 ) + h SV (Ry,i+1 + Ry,i-1 ) + 2 h v2 (40) для произвольного k 2. Из оценки (40) для четного k с учетом (18), (19) имеем ki δf1h ≈ 2 U V O(hk+2 ) + h P V O(hk+1 ) QV O(hk+1 ) + h SV O(hk+2 ) + = 2 h2 v2 2 h v2 = O(hk ) + O(hk ) + O(hk ) + O(hk+2 ) = O(hk ), (41) а для нечетного k с учетом (20), (21) имеем ki δf1h ≈ 402 2 U V O(hk+1 ) + h P V O(hk+2 ) QV O(hk+2 ) + h SV O(hk+1 ) + = 2 h2 v2 2 h v2 Оценка порядка аппроксимации матричного метода численного интегрирования краевых . . . = O(hk-1 ) + O(hk+1 ) + O(hk+1 ) + O(hk+1 ) = O(hk-1 ) (42) для всех i = 1, 2, . . . , n - 1. Исследование невязки (30) дало схожий результат: ki δf2h ≈ O(hk ) (43) ki δf2h ≈ O(hk-1 ) (44) для четного k и для нечетного k. Из оценки (41) следует, что при четном k вклады функций x, y в величиki сравнимы между собой, тогда как оценка (42) показывает, ну невязки δf1h что при нечетном k обусловленная функцией x составляющая невязки как имеющая более низкий порядок играет главенствующую роль по сравнению с составляющей, обусловленной функцией y. Практически аналогичной окаki с тем лишь отличием, что при нечетном k залась ситуация с невязкой δf2h главенствующую роль во вкладе в величину невязки играет функция y. В итоге в соответствии с (24) норму невязки в равенстве (23) запишем для четного k с учетом (41), (43) как δfhk = max k δf1h k δf2h = max k0 |, |δf k1 |, . . . , |δf kn | max |δf1h 1h 1h k0 k1 |, . . . , |δf kn | max |δf2h |, |δf2h 2h = max O(hk ) O(hk ) = = O(hk ), (45) и для нечетного k с учетом (42), (44) как δfhk = O(hk-1 ). (46) 2m+1 , где Анализ норм невязок (45), (46) показывает, что задачи L2m h , Lh m - натуральное число, имеют одинаковый порядок аппроксимации. Следовательно, на практике для уменьшения объема вычислений следует использовать задачи Lkh с четными значениями k. Действительно, число арифметических операций только для нахождения обратной матрицы от матрицы 2 размерности n методом Гаусса вычисляется по формуле 38 n3 - n2 - n6 - 1 [17]. Вывод о правомерности выбора четного значения k ниже будет подтвержден численным экспериментом. 2 ≈ -2h4 u имеет, как следует Значение алгебраического дополнения M26 1 2 , из (39), одинаковый порядок малости с алгебраическим дополнением M15 2 2 что свидетельствует о несущественной роли в использовании M15 или M26 ki или, что то же самое, в использовании bki 15 или b26 при вычислении невязок по формулам (29), (30). 4. Оценка погрешностей. При выполнении численного эксперимента в качестве суммарной оценки относительных погрешностей использованы следующие нормы: Dxk = n 2 i=0 (xi - [xi ]) n i=0 [xi ] · 100 %, Dyk = n 2 i=0 (yi - [yi ]) n i=0 [yi ] · 100 %, (47) 403 М а к л а к о в В. Н. которые можно трактовать как некий аналог коэффициента вариации в статистике, характеризующий меру разброса в процентах [18], а в качестве максимальной оценки абсолютных погрешностей: Exk = max xi - [xi ] , Eyk = max yi - [yi ] , i = 0, 1, . . . , n. (48) В качестве примера использована имеющая аналитическое решение система нелинейных ОДУ2 вида (2)  t2 + 2   x - tx - x - ty = 2t cos t, t2 (49)  1 x t2 - 4   x + +y + y = -2t sin t 2 t 2t2 с граничными условиями x(2π) = 0, x(3π) = 0, y(2π) = 4π 2 , y(3π) = -9π 2 . (50) В вычислениях использовались следующие параметры сетки: n = 15, h = 0.20944. Результаты численного эксперимента для краевой задачи (49) с граничными условиями (50) приведены в табл. 1, 2. В табл. 2 нормы Dxk , Dyk , Exk , Eyk для производных x (t), y (t) характеризуют суммарные оценки относительных погрешностей. Максимальные оценки абсолютных погрешностей, соответственно, вычислены по формулам (47), (48), в которых значения функций заменены на значения своих первых производных, найденных по формулам (10), (11), во внутренних узлах обла(j) (j) сти интегрирования; найденных при вычисленных x1 , y1 , j = 1, 2, . . . , k, по (10)-(13) при i = 1 по формулам h2 hk-1 (k) x1 - · · · + (-1)k-1 x , 2! (k - 1)! 1 h2 hk-1 (k) y0 = y1 - hy1 + y1 - · · · + (-1)k-1 y 2! (k - 1)! 1 x0 = x1 - hx1 + (51) (52) в левой границе t0 = a; найденных xn , yn по аналогичным (51), (52) разложениям в правой границе tn = b. Таблица 1 Значения погрешностей для решения краевой задачи (49), (50) [The values of the errors for the solution of the boundary value problem (49), (50)] k 2 3 4 5 6 7 Dxk , % Dyk , % Exk Eyk 404 1.62 · 10-1 1.21 · 10-1 3.59 · 10-1 2.19 · 10-1 1.47 · 10-1 2.13 · 10-1 3.47 · 10-1 3.56 · 10-1 4.21 · 10-4 6.96 · 10-4 1.01 · 10-3 1.15 · 10-3 1.70 · 10-4 4.49 · 10-4 4.34 · 10-4 7.21 · 10-4 1.57 · 10-5 9.85 · 10-6 3.45 · 10-5 1.77 · 10-5 1.60 · 10-5 1.09 · 10-5 3.52 · 10-5 1.97 · 10-5 Оценка порядка аппроксимации матричного метода численного интегрирования краевых . . . Таблица 2 Значения погрешностей для первых производных решения краевой задачи (49), (50) [The values of the errors for the first derivatives of the boundary value problem (49), (50)] k 2 3 4 5 6 7 Dxk , % Dyk , % Exk Eyk 4.53 · 10-1 3.43 · 10-1 7.32 · 10-1 8.20 · 10-1 1.98 · 10-1 1.50 · 10-1 4.07 · 10-1 3.90 · 10-1 1.31 · 10-3 1.09 · 10-3 2.27 · 10-3 2.57 · 10-3 2.80 · 10-4 2.90 · 10-4 5.09 · 10-4 6.94 · 10-4 3.78 · 10-5 2.34 · 10-5 7.09 · 10-5 5.66 · 10-5 3.86 · 10-5 2.41 · 10-5 7.08 · 10-5 5.88 · 10-5 Практически аналогичный характер изменения погрешностей (динамика и абсолютные значения) имел место для ряда систем ОДУ2, в частности для системы (1 + t)x + 2x + ty - 2y = 2 sin 2t, x + 2x - 2ty = 2(1 + t2 ) sin 2t с граничными условиями x(2π) = 0, x(3π) = 0, y(2π) = 2π, y(3π) = 3π. Анализ данных табл. 1, 2 свидетельствует, что суммарные относительные Dxk , Dyk , Dxk , Dyk и максимальные абсолютные Exk , Eyk , Exk , Eyk погрешности 2m+1 , имеющих одинаковый порядок аппроксимации, различазадач L2m h и Lh ются незначительно для любого натурального m = 1, 2, 3. Выводы. 1. Теоретически выявлены закономерности между порядком аппроксимации матричного метода и степенью k используемого многочлена Тейлора в разностных краевых задачах для систем линейных ОДУ2 с граничными условиями первого рода. Установлено, что порядок аппроксима2m+1 совпадает и равен 2m для ции разностных краевых задач L2m h и Lh натурального числа m. ki , i = 1, 2, . . . , 2. При четном k вклады функций x, y в величины невязок δf1h n - 1, сравнимы между собой; при нечетном k обусловленные функцией x составляющие невязок как имеющие более низкий порядок играют главенствующую роль по сравнению с составляющими, обусловленными функцией y. Практически аналогичной оказалась ситуация с невязками ki с тем лишь отличием, что при нечетном k главенствующую роль δf2h во вклады в величины невязок играет функция y. Декларация о финансовых и других взаимоотношениях. Исследование не имело спонсорской поддержки. Автор несет полную ответственность за предоставление окончательной версии рукописи в печать. Окончательная версия рукописи была одобрена автором. Автор не получал гонорар за статью.

About the authors

Vladimir N Maklakov

Samara State Technical University

Email: makvo63@yandex.ru
(Cand. Phys. & Math. Sci.; makvo63@yandex.ru), Associate Professor, Dept. of Higher Mathematics and Applied Informatics 244, Molodogvardeyskaya st., Samara, 443100, Russian Federation

References

Supplementary files