Numerical integration by the matrix method of boundary value problems for linear inhomogeneous ordinary differential equations of the third order with variable coefficients

Full Text

Предложенный в работе [1] метод, использующий средства матричного исчисления и численного интегрирования краевых задач для обыкновенных дифференциальных уравнений второго порядка (ОДУ2) с переменными коэффициентами позволяет удерживать произвольное число членов в разложении в ряд Тейлора искомого решения задачи, отказавшись при этом от аппроксимации производных конечными разностями. Известно, что использование конечных разностей приводит ко второму порядку аппроксимации при численном интегрировании краевых задач как для ОДУ2 [2-7], так и для ряда краевых задач для дифференциальных уравнений в частных производных [5-11]. Последнее обусловлено тем, что при аппроксимации производных было удержано всего три члена в разложении в ряд Тейлора искомого решения задачи, или, что то же самое, был использован многочлен Тейлора второй степени. Оценка порядка аппроксимации матричного метода при интегрировании краевых задач для ОДУ2 дана в [12], где показано, что для граничных условий первого, второго и третьего рода порядок аппроксимации пропорционален степени используемого многочлена Тейлора и меньше него на единицу во всех случаях, кроме граничных условий первого рода при четной степени используемого многочлена Тейлора - в этом случае порядок аппроксимации совпадает со степенью многочлена Тейлора. Поставим целью исследовать возможность использования матричного метода при численном интегрировании краевых задач в области интегрирования [a, b] для ОДУ3 s(t)x + r(t)x + p(t)x + q(t)x = f (t) (1) с граничными условиями x(a) = x0 , 154 x (a) = x0 , x(b) = xn , (2) Численное интегрирование матричным методом краевых задач. . . где x(t) - искомое решение; s(t), r(t), p(t), q(t), f (t) - заданные функции, дифференцируемые нужное число раз; x0 , x0 , xn - заданные числа, с последующим вычислением порядка аппроксимации разностной краевой задачи в зависимости от используемой степени многочлена Тейлора и числа узлов сетки в шаблоне. Далее будем придерживаться принятых в [5] обозначений: 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) - функция, являющаяся точным решением краевой задачи (1), (2); 3) [x]h - сеточная функция, совпадающая с точным решением в узлах сетки Dh ; 4) x(h) - искомая сеточная функция. Для краткости примем для любой функции обозначение ϕ(ti ) = ϕi , где ti - узел сетки Dh . В дальнейшем опустим индекс h в наименованиях сеточных функций [x]h , x(h) и будем особо оговаривать случаи, в которых используется функция x(t), являющаяся точным решением задачи. 1. Выбор узлов сетки при использовании четырехточечного шаблона. В матричном методе численного интегрирования краевых задач минимальная степень используемого многочлена Тейлора, что соответствует простейшему случаю, совпадает со старшей степенью производной в исследуемом ОДУ [1]; поэтому в начале остановимся на использовании многочленов Тейлора степени k = 3. В соответствии с матричным методом [1] при исследовании краевых задач для ОДУ2 в простейшем случае (при k = 2) необходимо составить систему линейных алгебраических уравнений (СЛАУ), в которую следует внести: а) два многочлена Тейлора степени k = 2, полученных из двух разложений в ряд Тейлора искомого точного решения x(t) в окрестностях слева и справа от некоторого внутреннего узла ti , i = 1, 2, . . . , n - 1, сетки Dh - последнее приводит к появлению трехточечного шаблона ti-1 , ti , ti+1 , где, как нетрудно заметить, узел ti является центральным узлом шаблона; б) исследуемое ОДУ2, записанное в центральном узле ti шаблона. Преобразование полученной замкнутой СЛАУ, как показано в [13], приводит к разностному уравнению традиционного метода сеток [5]. Отметим, что в простейшем случае число уравнений в замкнутой СЛАУ на единицу превышает число k. Наличие старшей производной третьей степени в левой части уравнения (1) указывает на то, что в простейшем случае СЛАУ краевой задачи (1), (2) при построении разностного уравнения должна содержать четыре уравнения с четырьмя неизвестными значениями искомой сеточной функции x, записанными или в узлах ti-2 , ti-1 , ti , ti+1 , i = 2, 3, . . . , n - 1, или в узлах ti-1 , ti , ti+1 , ti+2 , i = 1, 2, . . . , n - 2. В указанных двух шаблонах будем использовать три многочлена Тейлора, записанных в окрестностях слева и справа от узла ti , который в дальнейшем будем называть ведущим узлом шаблона. Первый из указанных выше четырехточечных шаблонов будем далее называть правым (ведущий узел ti расположен ближе к правой границе 155 М а к л а к о в В. Н., С т е л ь м а х Я. Г. шаблона) четырехточечным шаблоном или просто правым шаблоном, второй - левым шаблоном. Добавление ОДУ3 (1), записанного в ведущем узле ti , к трем многочленам Тейлора замкнет СЛАУ, позволяющую построить разностное уравнение. Выясним, следует ли отдавать предпочтение тому или иному из двух шаблонов. Для ведущих узлов ti , i = 1, 2, . . . , n-2, используем левые шаблоны. После выписывания многочленов Тейлора добавим к ним ОДУ3 (1), записанное в ведущем узле ti . В итоге получим следующую СЛАУ:  h3 h2   x - x = xi-1 , x - hx +  i i i   2! 3!3 i  2  h h  xi + hxi + xi + xi = xi+1 , (3) 2! 2 3!  3  h h x + 2hx + 22 x + 23 x = x ,   i i+2 i  2! i 3! i   qi xi + pi xi + ri xi + si xi = fi . В матричной форме система уравнений (3) в обозначениях   2 3     1 -h h2! - h3! xi xi-1   2 3 h h  1 h x  x  2! 3!  , A3i =  W 3i =  i  , G3i =  i+1  2  x xi+2 h h3  2 3 i  1 2h 2 2! 2 3!  xi fi qi p i ri si имеет вид A3i W 3i = G3i . Здесь и далее первый верхний индекс k (сейчас k = 3) означает степень используемого многочлена Тейлора P k , если речь не идет о показателях алгебраических степеней, порядках производных, символах обратных матриц и транспонировании; второй из пары верхних индексов в наименованиях матриц и их элементов означает номер ведущего узла сетки. Матрицы Aki , как и ранее в [12, 13], будем называть локальными матрицами. Предполагая существование обратной матрицы B 3i = (A3i )-1 от локальной матрицы A3i , найдем W 3i = (A3i )-1 G3i , или в координатной форме 3i 3i 3i xi = b3i 11 xi-1 + b12 xi+1 + b13 xi+2 + b14 fi , (4) 3i 3i 3i xi = b3i 21 xi-1 + b22 xi+1 + b23 xi+2 + b24 fi , 3i 3i 3i xi = b3i 31 xi-1 + b32 xi+1 + b33 xi+2 + b34 fi , 3i 3i 3i xi = b3i 41 xi-1 + b42 xi+1 + b43 xi+2 + b44 fi , (5) (6) (7) 3i где b3i lm - элементы матрицы B в ведущем узле ti . Выполненные выше действия для узла ti назовем для краткости процедурой построения (ПП) с использованием левого шаблона. Из равенств (4) составим систему разностных уравнений 3i 3i 3i -b3i 11 xi-1 + xi - b12 xi+1 - b13 xi+2 = b14 fi , 156 i = 1, 2, . . . , n - 2, (8) Численное интегрирование матричным методом краевых задач. . . для левого шаблона. Отметим, что в системе (8) число неизвестных (совпадающее с числом внутренних узлов сетки Dh и поэтому равное n - 1) на единицу превышает число разностных уравнений (совпадающее с числом ведущих узлов сетки). Однако пока не использовано второе граничное условие x (a) = x0 в (2). Для замыкания системы (8) в ведущем узле t1 выполним ПП с использованием левого шаблона, в которой вместо первого по счету многочлена Тейлора для функции x, как это было сделано в (3), запишем многочлен Тейлора для производной x . Тогда полученная таким образом система  h2   x1 - hx1 + x1 = x0 ,    2!   h3 h2  x1 + hx1 + x1 + x1 = x2 , (9) 2! 2 3!   h h3  2 3  x1 + 2hx1 + 2 x +2 x = x3 ,   2! 1 3! 1   q1 x1 + p1 x1 + r1 x1 + s1 x1 = f1 позволит построить, с учетом (2), недостающее в (8) разностное уравнение 31 31 31 -c31 11 x0 + x1 - c12 x2 - c13 x3 = c14 f1 , 31 = (A31 )-1 ; A31 - локальная матрица, при где c31 lm - элементы матрицы C построении которой использована система (9). В итоге получим следующую разностную краевую задачу, которую запишем в компактной символической форме в обозначениях, заимствованных из [5]: 3 , (10) L3h,L x = fh,L где принято  x1 c31 c31  12 13  - x -  2 31 31 31 x3 ,  c c c  14 14 14    31 31  x b b 1  12 13    b31 - b31 x2 - b31 x3 , 14 14 14 L3h,L x = 3i b x b3i b3i i  11 12 13  - x + - x - x , i = 2, 3, . . . , n - 3,  3i i-1 3i 3i i+1 3i i+2  b b b b  14 14 14 14    3(n-2) 3(n-2)  xn-2  b12 b11   x + - x , -  3(n-2) n-3 3(n-2) 3(n-2) n-1 b14 b14 b14  c31  11  f + 1  31 x0 ,  c  14     b31 11  f1 + 31 x0 , 3 b14 fh,L =   fi , i = 2, 3, . . . , n - 3,     3(n-2)   b13   fn-2 + 3(n-2) xn , b14 (11) (12) 157 М а к л а к о в В. Н., С т е л ь м а х Я. Г. а нижний индекс L означает, что были использованы только левые шаблоны. Решение задачи (10) дает значения искомой сеточной функции xi во внутренних узлах ti , i = 1, 2, . . . , n - 1, сетки Dh ; значения в граничных узлах t0 , tn определены граничными условиями (2). Отметим недостатки полученных расчетных формул для левых шаблонов. Во-первых, после нахождения решения задачи (10) равенства (5)-(7) позволяют вычислить производные вплоть до третьего порядка во всех внутренних узлах сетки, кроме узла tn-1 , который не являлся ведущим узлом ни разу. Во-вторых, узел t1 был ведущим дважды, а именно при построении матриц B 31 и C 31 , что дублирует формулы для вычисления производных в этом узле. При использовании в ПП правого шаблона вместо (10)-(12) будет получено 3 , (13) L3h,R x = fh,R где принято  32 c x c32  - 12 x1 + 2 - 13 x3 ,   c32 c32 c32  14 14 14    b32  x2 b32  12 13   - b32 x1 + b32 - b32 x3 , 14 14 14 L3h,R x = 3i 3i xi b b b3i - 11 x 12 13  - x + - x , i = 3, 4, . . . , n - 2,  i-2 i-1 3i 3i 3i 3i i+1  b b b b  14 14 14  14   3(n-1) 3(n-1)   b11 b12 xn-1   - 3(n-1) xn-3 - 3(n-1) xn-2 + 3(n-1) , b14 b b14  14 32 c   x0 , f2 + 11   c32  14     b32 11  f2 + 32 x0 , 3 b 14 fh,R =  fi , i = 3, 4, . . . , n - 2,     3(n-1)   b13   fn-1 + 3(n-1) xn , b14 а нижний индекс R означает, что были использованы только правые шаблоны. Решение задачи (13) дает значения искомой сеточной функции xi во внутренних узлах ti , i = 1, 2, . . . , n - 1, сетки Dh ; значения в граничных узлах t0 , tn определены граничными условиями (2). Полученные расчетные формулы для правых шаблонов обладают теми же недостатками, что и для левых шаблонов. С одной стороны, после нахождения решения задачи (13) равенства (5)-(7) позволяют вычислить производные вплоть до третьего порядка во всех внутренних узлах сетки, кроме узла t1 , который не являлся ведущим узлом ни разу. С другой стороны, узел t2 был ведущим дважды, а именно при построении матриц B 32 и C 32 , что дублирует формулы для вычисления производных в этом узле. Однако ситуацию довольно просто исправить - для этого достаточно вместо C 32 вычислить матрицу C 31 в узле t1 с помощью системы (9), при построении 158 Численное интегрирование матричным методом краевых задач. . . которой выполнена ПП с левым шаблоном, а затем использовать только правые шаблоны. В итоге получим следующую задачу, в которой использован смешанный вариант ПП: L3h x = fh3 , (14) где принято  c31 c31 x1  12 13  - x -  2 31 31 31 x3 ,  c c c  14 14 14     b32 x2 b32  12 13  - x + -   b32 1 b32 b32 x3 , 14 14 14 L3h x = b3i b3i xi b3i  11 12 13  - x - x + - xi+1 , i = 3, 4, . . . , n - 2,  i-2 i-1 3i 3i 3i 3i  b14 b14 b14  b14    3(n-1) 3(n-1)   b11 b12 xn-1   - 3(n-1) xn-3 - 3(n-1) xn-2 + 3(n-1) , b14 b14 b14  c31  11  f + 1  31 x0 ,  c  14     b32 11  f2 + 32 x0 , b14 fh3 =   fi , i = 3, 4, . . . , n - 2,     3(n-1)   b13   fn-1 + 3(n-1) xn . b14 Очевидно, что для граничных условий вида x(a) = x0 , x(b) = xn , (15) (16) x (b) = xn сначала следует (n - 2) раза использовать левый шаблон, а затем один раз использовать правый шаблон. Вопрос о влиянии вида шаблона на порядок аппроксимации разностной краевой задачи обсуждается ниже. 2. Численное интегрирование краевых задач для ОДУ3 при использовании четырехточечного шаблона. Исследуем разностную краевую задачу при различных значениях k. Для некоторого фиксированного k > 3, что приведет к увеличению числа слагаемых в многочленах Тейлора, в узле сетки t1 выполним ПП с использованием левого шаблона лишь с тем отличием от приведенной выше процедуры, что СЛАУ (9) дополним уравнениями q1 x1 + p1 x1 + r1 x1 + s1 x1 (r) (r) = f1 , r = 1, 2, . . . , k - 3, (17) полученными дифференцированием обеих частей ОДУ3 (1) и записанными в узле t1 ; в узлах сетки ti , i = 2, 3, . . . , n - 1 выполним ПП с использованием правого шаблона, дополнив систему вида (3) уравнениями (17), записанными в узлах ti . В итоге вместо (14)-(16) получим разностную краевую задачу Lkh x = fhk , (18) 159 М а к л а к о в В. Н., С т е л ь м а х Я. Г. где принято  ck1 ck1 x1  12 13  - x - x ,  2 31 k1 k1 3  c c c  14 14 14   k2   b12 x2 bk2  13   - bk2 x1 + bk2 - bk2 x3 , 14 14 14 Lkh x = ki ki b b xi bki  11 12 13  - x - x + - x , i = 3, 4, . . . , n - 2,  i-2 i-1  ki ki ki i+1  bki b b b  14 14 14 14   k(n-1) k(n-1)   b12 xn-1 - b11   k(n-1) xn-3 - k(n-1) xn-2 + k(n-1) , b14 b14 b14  k+1 k1  c1m (m-4) ck1  11   f + k1 x0 , f + 1  k1 1  c c  14 14  m=5   k+1 k2   b1m (m-4) bk2   f + f + 11 x ,  2  k2 2 k2 0 b b 14 14 m=5 fhk = k+1 ki  b1m (m-4)    f + fi , i = 3, 4, . . . , n - 2, i   bki  14  m=5   k+1 k(n-1) k(n-1)   b1m b13 (m-4)   f + k(n-1) xn ,  fn-1 + k(n-1) n-1 b14 m=5 b14 (19) (20) решение которой дает значения искомой сеточной функции xi во внутренних узлах ti , i = 1, 2, . . . , n - 1, сетки Dh ; значения в граничных узлах t0 , tn определены граничными условиями (2); производные во внутренних узлах ti могут быть вычислены по формулам k+1 (l-1) x1 (m-4) k1 k1 k1 = ck1 l1 x0 + cl2 x2 + cl3 x3 + cl4 f1 + ck1 lm f1 , (21) m=5 k+1 (l-1) xi (m-4) ki ki ki = bki l1 xi-2 + bl2 xi-1 + bl3 xi+1 + bl4 fi + bki lm fi , (22) m=5 где l = 2, 3, . . . , k + 1, i = 2, 3, . . . , n - 1. 3. Вычисление порядка аппроксимации разностной краевой задачи при использовании четырехточечного шаблона. Сеточная функция xi , i = 0, 1, . . . , n, являющаяся решением некоторой разностной краевой задачи, при подстановке в уравнения этой разностной краевой задачи обратит их в верные равенства. В [5] показано, что подстановка в уравнения задачи сеточной функции [xi ], отличающейся от xi , приведет к некоторому отличию от верных равенств. Эти отличия характеризуются невязкой δfhk [5]. Иными словами, подстановка [x] в Lkh x = fhk приведет к равенству Lkh [x] = fhk + δfhk . 160 (23) Численное интегрирование матричным методом краевых задач. . . В соответствии с [5] в качестве оценки величины невязки примем норму k k k k k δfhk = max |δfh0 |, |δfh1 |, |δfh2 |, . . . , |δfh,n-1 |, |δfhn | , (24) где первые две и последняя компоненты характеризуют меры отличий, появление которых обусловлено граничными условиями (2), остальные - оставшимися разностными уравнениями задачи (18). Отметим, что величины невяk |, |δf k |, появление которых обусловлено первым и третьим граничзок |δfh0 hn ными условиями в (2), как показано в [5], обращаются в нуль. Согласно [5, 7], разностная краевая задача аппроксимирует дифференциальную краевую задачу на точном решении x(t), если δfhk → 0 при h → 0. Chk , где C > 0, k > 0 - Если при этом имеет место неравенство |δfhk некоторые постоянные, не зависящие от h, то говорят, что имеет место аппроксимация порядка k относительно величины h. Вычислим порядок аппроксимации разностной краевой задачи (18). Исследование рассматриваемой задачи, в которой используем смешанный вариант ПП, начнем с ПП, использующей левый шаблон, т.е. исследуем ведущий узел t1 для граничного условия x (a) = x0 . При фиксированном k 3 вместо многочленов Тейлора используем точные равенства  k-1 h2  (k) k-1 h   [x ] - h[x ] + [x ] + · · · + (-1) [x1 ] = [x0 ] - R0k-1 , 1 1  1  2! (k - 1)!   h3 hk (k) h2 ] + ] + · · · + [x ] + h[x ] + [x [x [x1 ] = [x2 ] - R2k , 1 1 1 1  2! 3! k!    2 3 k   [x1 ] + 2h[x ] + 22 h [x ] + 23 h [x ] + · · · + 2k h [x(k) ] = [x3 ] - Rk , 1 1 1 3 2! 3! k! 1 где R0k-1 , R2k , R3k - дополнительные члены разложений в ряд Тейлора в форме Лагранжа [14]: Rik-1 = hk+1 (k+1) x (ξi ) = O(hk+1 ), (k + 1)! ξi ∈ (t1 , t3 ), i = 0, 2, 3. Тогда получим  k-1  h2 (k)  k-1 h  [x ] - h[x ] + [x ] + · · · + (-1) [x1 ] = [x0 ] - R0k-1 ,  1 1 1  2! (k - 1)!    2 3  h h hk (k)   [x ] + h[x ] + [x ] + [x ] + · · · + [x ] = [x2 ] - R2k ,  1 1 1 1   2! 3! k! 1    h2 h3 hk (k)  [x1 ] + 2h[x1 ] + 22 [x1 ] + 23 [x1 ] + · · · + 2k [x1 ] = [x3 ] - R3k , (25) 2! 3! k!    q1 [x1 ] + p1 [x1 ] + r1 [x1 ] + s1 [x1 ] = f1 ,      (4)   q1 [x1 ] + (q1 + p1 ) [x1 ] + (p1 + r1 ) [x1 ] + (r1 + s1 ) [x1 ] + s1 [x1 ] = f1 ,     ···    (k-3) (k) (k-3)  q1 [x1 ] + · · · + s1 [x1 ] = f1 . 161 М а к л а к о в В. Н., С т е л ь м а х Я. Г. В матричной форме система (25) в обозначениях  Ak1 0 -h 1 h2 2! h3 3! 3 23 h3! s1  h2  1 h 2!  2  1 2h 22 h2!  =  q1 p1 r1   q1 q + p p + r1 r1 + s1 1 1 1   ... ... ... ... (k-3) q1 ... ... ... k-1 h . . . (-1)k-1 (k-1)! hk k! k 2k hk! ... ... ... ... ... ...      ,     0 0 ... s1 (4) (k) [W k1 ] = [x1 ] [x1 ] [x1 ] [x1 ] [x1 ] . . . [x1 ]  , (k-3) [Gk1 ] = [x0 ]-R0k-1 [x2 ]-R2k [x3 ]-R3k f1 f1 . . . f1 имеет вид Ak1 [W k1 ] = [Gk1 ]. В предположении существования обратной матрицы C k1 = (Ak1 )-1 от локальной матрицы Ak1 найдем C k1 [Gk1 ] = [W k1 ]. Выпишем первое уравнение последнего матричного равенства: k k1 k k1 c11 [x0 ] - R0k-1 + ck1 12 [x2 ] - R2 + c13 [x3 ] - R3 + k+1 (m-4) + ck1 14 f1 + ck1 1m f1 = [x1 ], (26) m=5 k1 в узле t . где ck1 1 1m - элементы матрицы C Преобразуем равенство (26) к виду - k1 [x1 ] c12 ck1 ck1 13 11 [x ] + - [x ] - [x3 ] = 2 0 k1 k1 k1 ck1 c c c 14 14 14 14 k+1 = f1 + k-1 k k1 k ck1 ck1 + ck1 1m (m-4) 11 R0 12 R2 + c13 R3 f - . (27) ck1 1 ck1 14 m=5 14 Выполняя ПП с использованием правых шаблонов в оставшихся ведущих узлах сетки, получим k+1 - bki bki [xi ] bki bki 11 1m (m-4) [xi-2 ] - 12 [xi-1 ] + ki - 13 [xi+1 ] = fi + f - ki ki ki ki i b14 b14 b14 b14 b 14 m=5 - k ki k ki k bki 11 Ri-2 + b12 Ri-1 + b13 Ri+1 , bki 14 i = 2, 3, . . . , n - 1. (28) Равенства (27), (28) с учетом граничных условий (2) преобразуем к виду 162 Численное интегрирование матричным методом краевых задач. . . ck1 [x1 ] ck1 12 - k1 [x2 ] - 13 [x3 ] = k1 c14 c14 ck1 14 k+1 = f1 + - k-1 k k1 k ck1 ck1 ck1 + ck1 1m (m-4) 11 11 R0 12 R2 + c13 R3 f + x - , (29) 0 ck1 1 ck1 ck1 14 14 m=5 14 [x2 ] bk2 bk2 12 [x ] + - 13 [x3 ] = 1 bk2 bk2 bk2 14 14 14 k+1 = f2 + k k2 k k2 k bk2 bk2 bk2 1m (m-4) 11 11 R0 + b12 R1 + b13 R3 f + x - , (30) 0 2 bk2 bk2 bk2 14 14 m=5 14 k+1 bki bki [xi ] bki bki 11 13 1m (m-4) - ki [xi-2 ] - 12 [x ] + - [x ] = f + f - i-1 i+1 i ki ki ki b14 b14 b14 b14 bki i m=5 14 - k(n-1) - k ki k ki k bki 11 Ri-2 + b12 Ri-1 + b13 Ri+1 , bki 14 k+1 k(n-1) b11 i = 3, 4, . . . , n - 2, (31) k(n-1) b12 [xn-1 ] b1m (m-4) [x ] - [x ] + = f + f + n-3 n-2 n-1 k(n-1) k(n-1) k(n-1) k(n-1) n-1 b14 b14 b14 b m=5 14 k(n-1) + b13 x - k(n-1) n b14 k(n-1) b11 k(n-1) k Rn-3 + b12 k(n-1) k Rn-2 + b13 k(n-1) Rnk . (32) b14 Отбрасывание последних дробей в равенствах (29)-(32), что равносильно переходу от точного решения [xi ] к искомому приближенному xi , приведет эти равенства к уравнениям разностной краевой задачи (18)-(20). Следовательно, в соответствии с (23), последние дроби в равенствах (29)-(32) характеризуют величину невязки в уравнениях задачи, при построении которых были использованы ведущие узлы ti , i = 1, 2, . . . , n - 1. В итоге для рассматриваемой задачи имеем  k1 k-1 k k1 k c11 R0 + ck1  12 R2 + c13 R3   - ,   ck1  14    k k2 k k2 k  bk2  11 R0 + b12 R1 + b13 R3   - ,   bk2 14 δfhk = k ki k ki k bki  11 Ri-2 + b12 Ri-1 + b13 Ri+1   - , i = 3, 4, . . . , n - 2,   bki  14     k(n-1) k k(n-1) k k(n-1) k  b11 Rn-3 + b12 Rn-2 + b13 Rn    - ,  k(n-1) b14 163 М а к л а к о в В. Н., С т е л ь м а х Я. Г. или  k1 k-1 k k1 k c11 R0 + ck1  12 R2 + c13 R3  , -   ck1 14 k (33) δfh = k ki k ki k  bki  11 Ri-2 + b12 Ri-1 + b13 Ri+1  - , i = 2, 3, . . . , n - 1, bki 14 где, напомним, первая компонента характеризует величину невязки, появление которой обусловлено вторым граничным условием в (2). Непосредственными вычислениями можно для любого k > 3 убедиться в справедливости оценок ki 3i ≈ (si )k-1 M1j , M1j i = 2, 3, . . . , n - 1, j = 1, 2, 3, 4, (34) ki - алгебраическое дополнение элемента aki матрицы (Aki ) , полученгде M1j 1j ной ПП, использующей правый шаблон. Воспользуемся в дальнейшем известными свойствами определителя [15]. Действительно, пренебрегая старшими степенями, например при j = 1, имеем ki M11 = -h h h2 2! 3 - h3! h2 2! h3 3! ... ... hk-1 (-1)k-1 (k-1)! k (-1)k hk! k k-1 hk-1 (k-1)! hk k! = h k! pi qi + pi ri pi + ri si ri + si ... ... 0 0 0 0 . . . ui . . . vi . . . wi ... ... = . . . zi ... si k-1 bm hm + (-1)k+1 m=1 (k-1)i cm hm + si M11 (k-1)i ≈ si M11 , m=1 где bm , cm - коэффициенты, не зависящие от h; ui , vi , wi , zi - некоторые функции от qi , pi , ri , si и их производных. Повторное неоднократное использование последней формулы приводит к (34) при j = 1. Аналогично доказывается справедливость оценок (34) при j = 2, 3. ki : Рассмотрим M14 ki M14 = 164 hk k! -2h 2 22 h2! 3 -23 h3! =- ... hk-1 k-1 (-2) (k-1)! -h h h2 2! 3 - h3! h2 2! h3 3! ... ... q i + pi pi + ri ri + si ... 0 hk-1 hk-1 (-1)k-1 (k-1)! (k-1)! k k hk (-2)k hk! (-1)k hk! k! 2k-3 2k-3 2k-3 -2k bm hm + cm hm -(-1)k dm hm m=3 m=3 m=3 0 . . . ui . . . vi . . . wi = ... . . . zi ... si (k-1)i +si M14 (k-1)i ≈ si M14 , Численное интегрирование матричным методом краевых задач. . . где bm , cm , dm - коэффициенты, не зависящие от h. Повторное неоднократное использование последней формулы приводит к (34) при j = 4. На основании оценок (34) и очевидных равенств bki 1j bki 14 имеем bki 1j bki 14 ≈ 3i M1j 3i M14 = ki M1j ki M14 , j = 1, 2, 3 i = 2, 3, . . . , n - 1, , j = 1, 2, 3. (35) Вывод оценок ck1 1j ck1 14 ≈ 31 M1j 31 M14 , j = 1, 2, 3, (36) 31 , j = 1, 2, 3, 4 - алгебраическое дополнение элемента a31 матрицы где M1j 1j (A31 ) , полученной ПП, использующей левый шаблон, для ведущего узла t1 оказался аналогичным выводу формул (35). Невязка (33), с учетом оценок (35), (36) примет вид  31 Rk + M 31 Rk   M 31 Rk-1 + M12  2 13 3 k , - 11 0 , δfh1  31 M14 k (37) = δfh = 3i Rk +M 3i Rk +M 3i Rk  k  M11 12 i-1 13 i+1 i-2 δfhi , i = 2, 3, . . . , n - 1   . -  3i M14 3i , i = 1, 2, . . . , n - 1, j = 1, 2, 3, 4. Пренебрегая старВычислим оценки M1j шими степенями, для i = 1 имеем 31 M11 = h 2h h2 2! h3 3! 22 h2 2! 23 h3 3! p1 h5 r1 = s1 h3 - r1 h4 + p1 ≈ s1 h3 , 3! s1 h3 7 + p1 h4 ≈ -4s1 h2 , 3 3 3 3 h 5 3 11 31 = s1 h2 + r1 - p1 h4 ≈ s1 h2 , M14 = h5 ; 2 3 12 2 6 31 M12 = -4s1 h2 + r1 31 M13 (38) (39) (40) и для i = 2, 3, . . . , n - 1: -h 3i M11 2 = h2! 3 - h3! 3i M12 = h h2 2! h3 3! pi h5 ri = -si h3 + pi ≈ -si h3 , 3! si 3si h3 + ri h4 + pi h5 ≈ 3si h3 , 3i M13 = si h3 + ri h4 + pi h5 3 ≈ si h3 , 3i M14 = -h6 . (41) (42) (43) Вычислим невязки (37). Используя оценки (38)-(40), для ведущего узла t1 найдем 165 М а к л а к о в В. Н., С т е л ь м а х Я. Г. k δfh1 =- 31 Rk-1 + M 31 Rk + M 31 Rk M11 12 2 13 3 0 ≈ 31 M14 3 2s1 h3 R0k-1 - 8s1 h2 R2k + 3s1 h2 R3k = 11h5 6s1 O(hk ) -24s1 O(hk+1 ) + 9s1 O(hk+1 ) =- - = 11h2 11h3 = O(hk-2 ) + O(hk-2 ) + O(hk-2 ) = O(hk-2 ) (44) ≈- и, используя оценки (41)-(43), для ведущих узлов ti , i = 2, 3, . . . , n - 1 найдем k δfhi =- 3i Rk 3i k 3i k M11 i-2 + M12 Ri-1 + M13 Ri+1 ≈ 3i M14 k k k -si h3 Ri-2 + 3si h3 Ri-1 + si h3 Ri+1 = -h6 -si O(hk+1 ) + 3si O(hk+1 ) + si O(hk+1 ) = = h3 = O(hk-2 ) + (hk-2 ) + O(hk-2 ) = O(hk-2 ). (45) ≈- В итоге в соответствии с (24) из оценок (44), (45) для рассматриваемой задачи имеем δfhk = O(hk-2 ), откуда следует, что порядок аппроксимации пропорционален степени используемого многочлена Тейлора и меньше нее на две единицы независимо от четности k, в отличие от ОДУ2 с граничными условиями первого рода [12], где порядок аппроксимации задачи оказался зависимым от четности k. При исследовании в [12] краевых задач для ОДУ2 с граничными условиями второго и третьего рода оказалось, что именно наличие производной в граничных условиях понизило на единицу порядок аппроксимации всей задачи, чего не оказалось в рассмотренном выше случае, как это следует из оценок (44), (45). Путем непосредственных вычислений можно убедиться в справедливости а) оценок (35) для левого шаблона и оценок (36) для правого шаблона; б) оценок (38)-(43) с точностью до постоянного сомножителя для указанных выше шаблонов. Следствием этого является независимость порядка аппроксимации разностной краевой задачи (1), (2) от выбора формы четырехточечного шаблона. 4. Численное интегрирование краевых задач для ОДУ3 при использовании пятиточечного шаблона. Выше исследована возможность численного интегрирования матричным методом [1] краевых задач для линейных неоднородных обыкновенных дифференциальных уравнений третьего порядка с переменными коэффициентами (1) с граничными условиями (2). Установлен порядок аппроксимации метода при использовании в нем четырехточечного шаблона. 166 Численное интегрирование матричным методом краевых задач. . . Поставим целью исследование возможности численного интегрирования матричным методом краевых задач для ОДУ3 при использовании пятиточечного шаблона с последующим вычислением порядка аппроксимации в зависимости от степени k многочлена Тейлора. Решение остановиться на пятиточечном шаблоне связано лишь с тем, что производные третьего порядка, вычисленные с помощью отношения центральных разностей (что соответствует пятиточечному шаблону), имеют второй порядок аппроксимации относительно шага h, тогда как вычисленные с помощью отношения конечных разностей слева и справа (что соответствует четырехточечному шаблону) - первый порядок [6]. В соответствии с матричным методом составим СЛАУ для пятиточечного шаблона ti-2 , ti-1 , ti , ti+1 , ti+2 , i = 2, 3, . . . , n - 2, в которую при фиксированном k 4 внесем: а) четыре многочлена Тейлора степени k, полученных из четырех разложений в ряд Тейлора искомого точного решения x(t) в окрестностях слева и справа от некоторого внутреннего узла ti (центрального узла шаблона), i = 2, 3, . . . , n - 2, сетки Dh ; б) уравнения qi xi + pi xi + ri xi + si xi (r) (r) = fi , r = 0, 1, . . . , k - 4, полученные дифференцированием обеих частей ОДУ3 (1) и записанные в узле ti . В итоге получим СЛАУ  3 k 2  2 h x - 23 h x + · · · + (-2)k h x(k) = x  x - 2hx + 2  i-2 , i i  2! i 3! i k! i      k  h2 h3  k h x(k) = x  x - hx + x - x + · · · + (-1)  i i-1 , i i i   2! 3! k! i      h2 h3 hk (k)   x + hx + x + x + · · · + x = xi+1 ,  i i   2! i 3! i k! i   3 k 2 2 h x + 23 h x + · · · + 2k h x(k) = x  x + 2hx + 2 i+2 , i i   2! i 3! i k! i      qi xi + pi xi + ri xi + si xi = fi ,       (IV )   qi xi + (pi + qi )xi + (ri + pi )xi + (si + r)xi + si xi = fi ,       ...      (k-4) (k-1) (k-4) qi xi + · · · + s i xi = fi . (46) В матричной форме система уравнений (46) в обозначениях 167 М а к л а к о в В. Н., С т е л ь м а х Я. Г.  Aki 2 22 h2! 3 (-2)3 h3! -2h  1   h2  1 -h 2!    h2  h 2! = 1  2  2h 22 h2!  1   qi pi ri   ... ... ... (k-3) qi ... ... hk (-2)k ... 3 - h3! ... h3 3! ... 3 23 h3! si ... ... k!   hk  k  (-1) k!   hk   , k!   k h  2k  k!   0  ...  ... ... ... ... (k) W ki = xi xi xi xi . . . xi  (47) si , (k-4) Gki = xi-2 xi-1 xi+1 xi+2 fi fi fi . . . fi . имеет вид Aki W ki = Gki . Здесь и далее второй из пары верхних индексов в наименованиях матриц и их элементов означает номер центрального узла пятиточечного шаблона. Предполагая существование обратной матрицы B ki = (Aki )-1 от локальной матрицы Aki , найдем W ki = (Aki )-1 Gki , или в координатной форме (i = 2, 3, . . . , n - 2) k+1 xi = bki 11 xi-2 + bki 12 xi-1 + bki 13 xi+1 + bki 14 xi+2 + bki 15 fi (m-5) , (48) (m-5) , (49) bki 1m fi + m=6 k+1 ki ki ki ki xi = bki 21 xi-2 + b22 xi-1 + b23 xi+1 + b24 xi+2 + b25 fi + bki 2m fi m=6 k+1 ki ki ki ki xi = bki 31 xi-2 + b32 xi-1 + b33 xi+1 + b34 xi+2 + b35 fi + (m-5) , (50) (m-5) (51) bki 3m fi m=6 ... (k) xi = bki (k+1)1 xi-2 + bki (k+1)2 xi-1 + bki (k+1)3 xi+1 + k+1 + bki (k+1)4 xi+2 + bki (k+1)5 fi bki (k+1)m fi + , m=6 ki в узле t . где bki i lm - элементы матрицы B Из равенств (48), которые являются разностными уравнениями для пятиточечного шаблона ti-2 , ti-1 , ti , ti+1 , ti+2 , составим СЛАУ и запишем ее в компактной символической форме: Lkh x = fhk , 168 (52) Численное интегрирование матричным методом краевых задач. . . где  k2 x2 bk2 bk2 b12  13 14   x + - x - x , - 1 3  k2 k2 k2 k2 4  b b b b  15 15 15 15     bki bki xi bki bki 11 12 13 14 k - x - x + - x - x , i = 3, 4, . . . , n - 3, Lh x = ki i-2 ki i-1 ki ki i+1 ki i+2 b b b b b  15 15 15 15 15    k(n-2) k(n-2)   bk(n-2) b12 xn-2 b13  11   - x - x + - x ,  k(n-2) n-4 k(n-2) n-3 k(n-2) k(n-2) n-1 b15 b15 b15 b15 (53)  k+1 k2  b1m (m-5) bk2   f2 + 11 x0 , f2 +  k2  b15 bk2  15 m=6    k+1 ki  b1m (m-5) k fh = fi + (54) f , i = 3, 4, . . . , n - 3, ki i  b  15 m=6    k(n-2) k+1 k(n-2)   b14 b1m  (m-5)  f + k(n-2) xn .  fn-2 + k(n-2) n-2 b15 m=6 b15 Задачу (52) для краткости будем далее обозначать как Lkh . Отметим, что разностная краевая задача (52) содержит (n - 3) уравнения с (n-1) неизвестными, т.е. для замыкания системы необходимо составить еще два уравнения. В [5] для оценки порядка аппроксимации обоснована целесообразность разбиения разностной краевой задачи на подзадачи (подсистемы); поэтому в отдельные подзадачи выделим систему (52) и систему, содержащую недостающие два уравнения, а затем отдельно вычислим порядки аппроксимации каждой из этих двух подзадач. Для учета граничного условия x (a) = x0 сохраним все приведенные выше выкладки для задачи (52) с тем лишь отличием, что при составлении СЛАУ вместо первого приближенного равенства в (46) при i = 2 используем следующий многочлен: x2 - 2hx2 + 22 hk-1 (k) h2 x2 + · · · + (-2)k-1 x = x0 . 2! (k - 1)! 2 (55) Для составления второго недостающего уравнения задачи Lkh сформируем фиктивное граничное условие в узле t0 следующим образом. Из точного равенства s0 x0 + r0 x0 + p0 x0 + q0 x0 = f0 найдем s0 x0 + r0 x0 = f0 - p0 x0 - q0 x0 = X0 , (56) где число X0 может быть вычислено с использованием первых двух равенств граничных условий (2). Равенство (56) будем трактовать как фиктивное граничное условие в узле t0 . 169 М а к л а к о в В. Н., С т е л ь м а х Я. Г. Имеем следующие многочлены Тейлора: h2 (IV ) hk-2 (k) x2 + · · · + (-2)k-2 x = x0 , 2! (k - 2)! 2 h2 (V ) hk-3 (k) (IV ) x2 - 2hx2 + 22 x2 + · · · + (-2)k-3 x = x0 . 2! (k - 3)! 2 x2 - 2hx2 + 22 (57) (58) Умножая обе части многочлена (57) на r0 , многочлена (58) - на s0 и складывая полученное, с учетом (56) получим фиктивное граничное условие: r0 x2 + (s0 - 2hr0 )x2 + · · · + (-1)k-1 s0 (2h)k-3 (2h)k-2 (k) - r0 x = X0 . (59) (k - 3)! (k - 2)! 2 Для учета фиктивного граничного условия (56) сохраним все приведенные выше выкладки для задачи (52) с тем лишь отличием, что в СЛАУ (46) при i = 2 вместо первого приближенного равенства используем соотношение (59). Указанные действия позволяют сформировать последнее недостающее разностное уравнение. В итоге вторую подзадачу в компактной символической форме запишем как (60) lhk x = ghk , где  k2 x2 ck2 ck2 c12  13 14  x + - x - x , -  1 3  ck2 k2 k2 k2 4 c c c 15 15 15 15 k lh x =  x2 dk2 dk2  dk2 13 14  - 12 x + - x - x , 1 k2 k2 3 k2 4 dk2 d d d 15 15 15 15  k+1 k2  c1m (m-5) ck2 11   f + k2 x0 , f +  k2 2  2 c c 15 15 m=6 ghk = k+1 k2   d1m (m-5) dk2   f + f + 11 X .  2 k2 2 k2 0 d d 15 15 m=6 (61) (62) Здесь ck2 1m - элементы обратной матрицы от матрицы СЛАУ (46) при i = 2, которая вычислена для граничного условия x (a) = x0 ; dk2 1m - элементы обратной матрицы от матрицы СЛАУ (46) при i = 2, которая вычислена для фиктивного граничного условия (59). Решение разностной краевой задачи (52), (60) дает значения искомой сеточной функции xi во внутренних узлах ti , i = 1, 2, . . . , n - 1, сетки Dh ; значения в граничных узлах t0 , tn определены граничными условиями (2); равенства (49)-(51) позволят при найденных xi вычислить производные вплоть до порядка k в узлах ti , i = 2, 3, . . . , n - 2. 5. Вычисление порядка аппроксимации разностной краевой задачи при использовании пятиточечного шаблона. Вычислим порядок аппроксимации разностной краевой задачи (52), (60). 170 Численное интегрирование матричным методом краевых задач. . . Для оценки величины невязки задачи (52) в соответствии с [5] примем норму k k k δfhk = max |δfh2 |, |δfh3 |, . . . , |δfh,n-2 | , где компонента, второй нижний индекс которой, совпадающий с номером центрального узла шаблона, характеризует меру отличий, появление которой обусловлено уравнением задачи (52) c номером на единицу меньше этого второго нижнего индекса. Для задачи (60) - норму k k δghk = max |δghc |, |δghd | , (63) где нижний индекс c означает принадлежность компоненты к первому уравнению задачи (60), нижний индекс d - ко второму уравнению. Первая компонента характеризует меру отличий, появление которой обусловлено вторым граничным условием в (2), вторая - фиктивным граничным условием (59). Норму всей разностной краевой задачи (52), (60) в соответствии с [5] запишем как (64) δf = max δfhk , δghk . Исследуем задачу (52). При фиксированном k 4 сохраним все приведенные выше выкладки для задачи (52) с тем лишь отличием, что в СЛАУ (46) вместо четырех первых приближенных равенств используем точные равенства  h3 hk (k) h2   k ,  [xi ] - 2h[xi ] + 22 [xi ] - 23 [xi ] + · · · + (-2)k [xi ] = [xi-2 ] - Ri-2   2! 3! k!     k   h2 h3  k h [x(k) ] = [x k  [x ] - h[x ] + [x ] - [x ] + · · · + (-1) i i i-1 ] - Ri-1 ,  i 2! 3! i k! i (65)   h2 h3 hk (k)  k  [xi ] + h[xi ] + [xi ] + [xi ] + · · · + [xi ] = [xi+1 ] - Ri+1 ,   2! 3! k!       h2 h3 hk  k , [xi ] + 2h[xi ] + 22 [xi ] + 23 [xi ] + · · · + 2k [x(k) ] = [xi+2 ] - Ri+2 2! 3! k! i k , Rk , Rk , Rk - дополнительные члены разложений в ряд Тейлогде Ri-2 i+2 i+1 i-1 ра в форме Лагранжа [14]: Rjk = hk+1 (k+1) x (ξj ) = O(hk+1 ), (k + 1)! ξj ∈ (ti , ti+2 ), j = i - 2, i - 1, i + 1, i + 2. В итоге вместо (52) получим задачу (23), в которой левая часть сохранит вид (53) с заменой искомой сеточной функции xi на [xi ], а вместо (54) получим  k+1 k2  b1m (m-5) bk2  k f2 + f + 11 x + δfh2 ,  k2 2 k2 0  b b  15 15  m=6   k+1 ki  b1m (m-5) k k k fh + δfh = fi + fi + δfhi , i = 3, 4, . . . , n - 3, ki  b15  m=6    k(n-2) k+1 k(n-2)   b14 b1m (m-5)  k  f + f + x + δfh(n-2) ,  n-2 k(n-2) n-2 k(n-2) n b15 m=6 b15 171 М а к л а к о в В. Н., С т е л ь м а х Я. Г. где  k2 k k2 k k2 k k2 k  - b11 R0 + b12 R1 + b13 R3 + b14 R4 ,    bk2  15  ki k  ki Rk ki k ki k  b11 R + b 12 i-1 + b13 Ri+1 + b14 Ri+2 i-2 , i = 3, 4, . . . , n - 3, δfhk = - bki  15    k(n-2) k k(n-2) k k(n-2) k k(n-2) k  b11 Rn-4 + b12 Rn-3 + b13 Rn-1 + b14 Rn    - ,  k(n-2) b15 или δfhk = - k ki k ki k ki k bki 11 Ri-2 + b12 Ri-1 + b13 Ri+1 + b14 Ri+2 , bki 15 i = 2, 3, . . . , n - 2. (66) Исследуем задачу (60). При фиксированном k 4 сохраним все приведенные выше выкладки для задачи (60) с тем лишь отличием, что в СЛАУ (46) вместо второго, третьего, четвертого приближенных равенств используем второе, третье, четвертое точные равенства (65), а вместо первого и второго приближенных равенств используем следующее: а) при построении первого уравнения задачи (60) - точное равенство с учетом второго граничного условия в (2) вида [x2 ] - h[x2 ] + hk-1 h2 (k) [x2 ] + · · · + (-1)k-1 [x ] = [x0 ] - R0k-1 ; 2! (k - 1)! 2 б) при построении второго уравнения задачи (60) - точное равенство с учетом фиктивного граничного условия (59) вида r0 [x2 ] + (s0 - 2 hr0 )[x2 ] + · · · + + (-1)k-1 s0 (2h)k-2 (k) (2h)k-3 - r0 [x ] = X0 - s0 R0k-3 - r0 R0k-2 , (k - 3)! (k - 2)! 2 при составлении которого вместо (57), (58) используем точные равенства h2 (IV ) hk-2 (k) [x2 ] + · · · + (-2)k-2 [x ] = [x0 ] - R0k-2 , 2! (k - 2)! 2 h2 (V ) hk-3 (IV ) (k) [x2 ] - 2h[x2 ] + 22 [x2 ] + · · · + (-2)k-3 [x ] = [x0 ] - R0k-3 . 2! (k - 3)! 2 [x2 ] - 2h[x2 ] + 22 В итоге вместо (60) получим задачу lhk [x] = ghk + δghk , 172 Численное интегрирование матричным методом краевых задач. . . в которой левая часть сохранит вид (61) с заменой искомой сеточной функции xi на [xi ], а вместо (62) получим  k+1 k2  c1m (m-5) ck2  k  f2 + 11 x0 + δghc , f +  2  k2 k2 c15 c15 k k m=6 gh + δgh = k+1 k2  d1m (m-5) dk2  k   f + f + 11 X0 + δghd ,  2 dk2 2 dk2 m=6 15 15 где δghk =  k2 , δghc  k2 δghd  k k2 k k2 k ck2 Rk-1 +ck2  12 R1 + c13 R3 + c14 R4  - 11 0 ,  ck2 15 = k2 (s Rk-3 - r Rk-2 ) + dk2 Rk + dk2 Rk + dk2 Rk  d11 0 0 0 0  12 1 13 3 14 4  . - k2 d15 (67) Вычислим оценки невязок (66), (67) задач (52), (60) соответственно. Оценки ki 4i M1j ≈ (si )k-1 M1j , i = 1, 2, . . . , n - 1, j = 1, 2, . . . , 5, (68) аналогичные оценкам (34), полученным при использовании четырехточечного шаблона, имеют место и для матриц (Aki ) , где Aki - локальная матрица (47), и для двух транспонированных матриц (47) при i = 2, в которых первые строки заменены коэффициентами разложений (55), (59) соответственно, что выше было выполнено при построении матриц C k2 , Dk2 . Отсутствие в настоящей работе полного вывода формулы (68) обусловлено лишь громоздкостью выкладок. При исследовании краевых задач для ОДУ2 вывод формул, аналогичных (68), дан в [12]. На основании оценок (68) и очевидных равенств bki 1j bki 15 имеем bki 1j bki 15 ≈ 4i M1j 4i M15 , = ki M1j ki M15 , j = 1, 2, 3, 4 i = 1, 2, . . . , n - 1, j = 1, 2, 3, 4. (69) 4i , i = 1, 2, . . . , n-1, j = 1, 2, . . . , 5. Используя известВычислим оценки M1j ные свойства определителя [15] и пренебрегая старшими степенями, имеем 4i M11 = -h h 2h h2 2! 3 - h3! h4 4! h2 2! h3 3! h4 4! 22 h2 2! 23 h3 3! 24 h4 4! p1 ri h7 h7 = (6si + ri h - pi h2 ) ≈ (6si + ri h) , (70) si 12 12 0 h7 h7 ≈ (-3si - 4ri h) , 3 3 7 7 h h = (3si - 4ri h - 2pi h2 ) ≈ (3si - 4ri h) , 3 3 4i M12 = (-3si - 4ri h + 2pi h2 ) 4i M13 (71) (72) 173 М а к л а к о в В. Н., С т е л ь м а х Я. Г. 4i M14 = (-6si + ri h + pi h2 ) h7 h7 ≈ (-6si + ri h) , 12 12 4i M15 = qi h10 . (73) Вычислим оценки невязок (66) задачи (52). Для рядов Тейлора xi-1 = xi - hxi + h3 h2 xi - xi + · · · = 2! 3! k (-1)m = m=0 xi+1 = xi + hxi + hm (m) x + m! i ∞ (-1)m m=k+1 hm (m) k k x = Pi-1 + Ri-1 , m! i h2 h3 xi + xi + · · · = 2! 3! k = m=0 hm (m) x + m! i ∞ m=k+1 hm (m) k k x = Pi+1 + Ri+1 , m! i k , Rk - дополнительные члены: где Ri-1 i+1 Rjk = hk+1 (k+1) x (ξj ) = O(hk+1 ), (k + 1)! ξj ∈ (ti , ti+1 ), j = i - 1, i + 1, в [12] показана справедливость следующих формул: k k Ri+1 - Ri-1 = O(hk+2 ), k k Ri+1 + Ri-1 = O(hk+1 ) (74) k k Ri+1 + Ri-1 = O(hk+2 ) (75) для нечетного k и k k Ri+1 - Ri-1 = O(hk+1 ), для четного k. С учетом (69)-(73) из (66) для i = 2, 3, . . . , n - 2 имеем δfhki = - k ki k ki k ki k bki 11 Ri-2 + b12 Ri-1 + b13 Ri+1 + b14 Ri+2 ≈ bki 15 k k h7 (6si + ri h)Ri-2 + (-6si + ri h)Ri+2 - 12qi h10 k k h7 (-3si - 4ri h)Ri-2 + (3si - 4ri h)Ri+2 - = 3qi h10 k k ) + r h(Rk k -6si (Ri+2 - Ri-2 i i+2 + Ri-2 ) =- - 12qi h3 k k ) - 4r h(Rk k 3si (Ri+1 - Ri-1 i i+1 + Ri-1 ) - . (76) 3qi h3 ≈- Подстановка оценок (74) в (76) дает 174 Численное интегрирование матричным методом краевых задач. . . δfhki ≈ - -6si O(hk+2 ) + ri hO(hk+1 ) 3si O(hk+2 ) - 4ri hO(hk+1 ) - = 12qi h3 3qi h3 = O(hk-1 ) + O(hk-1 ) + O(hk-1 ) + O(hk-1 ) = O(hk-1 ), откуда δfhk = O(hk-1 ) (77) для нечетного k. Подстановка оценок (75) в (76) дает δfhki ≈ - -6si O(hk+1 ) + ri hO(hk+2 ) 3si O(hk+1 ) - 4ri hO(hk+2 ) - = 12qi h3 3qi h3 = O(hk-2 ) + O(hk ) + O(hk-2 ) + O(hk ) = O(hk-2 ), откуда δfhk = O(hk-2 ) (78) для четного k. и L2m Из (77), (78) следует, что при m 3 задачи L2m-1 h имеют одинаh ковый порядок аппроксимации. Отметим, что, как показано в [12, 18], для ОДУ2 и систем ОДУ2 с граничными условиями первого рода ситуация ока2m+1 имели одинаковый порядок залась противоположной - задачи L2m h и Lh аппроксимации при m 1. Вычислим оценки невязок (67) задачи (60). Для матрицы C k2 на основании оценок (68) и очевидных равенств ck2 1j ck2 15 имеем ck2 1j ck2 15 = ≈ k2 M1j , j = 1, 2, 3, 4 , j = 1, 2, 3, 4. k2 M15 42 M1j 42 M15 (79) 42 , j = 1, 2, . . . , 5. Пренебрегая старшими степенями, Вычислим оценки M1j найдем h7 h7 h7 ≈ (6s2 + r2 h) ≈ s2 , 12 12 2 6 6 h h 42 M12 ≈ -(-3s2 - 88r2 h) ≈ s2 , 9 3 6 6 h 11h 42 M13 ≈ (-11s2 - 12r2 h) ≈ -s2 , 4 4 h6 7h6 340h10 42 ≈ -(-7s2 + r2 h) ≈ s2 , M15 = q2 . 6 6 3 42 M11 = (6s2 + r2 h - p2 h2 ) 42 M14 (80) (81) (82) (83) С учетом (79)-(83) вычислим оценку первой компоненты в (67): k2 δghc =- k-1 k2 Rk + ck2 Rk + ck2 Rk ck2 + c12 11 R0 1 13 3 14 4 ≈ k2 c15 175 М а к л а к о в В. Н., С т е л ь м а х Я. Г. s2 h6 (6hR0k-1 + 4R1k - 33R3k + 14R4k ) = 1360q2 h9 = O(hk-2 ) + O(hk-2 ) + O(hk-2 ) + O(hk-2 ) = O(hk-2 ). (84) ≈- Для матрицы Dk2 на основании оценок (68) и очевидных равенств dk2 1j dk2 15 имеем dk2 1j dk2 15 = ≈ k2 M1j , j = 1, 2, 3, 4 , j = 1, 2, 3, 4. k2 M15 42 M1j 42 M15 (85) 42 , j = 1, 2, . . . , 5. Пренебрегая старшими степенями, Вычислим оценки M1j найдем h7 h7 h7 ≈ (6s2 + r2 h) ≈ s2 , 12 12 2 ≈ (s0 + 3r0 h)3s2 h4 ≈ 3s0 s2 h4 , 42 M11 = (6s2 + r2 h - p2 h2 ) 42 M12 42 M13 4 4 ≈ (3s0 + 4r0 h)2s2 h ≈ 6s0 s2 h , 42 M14 ≈ -(s0 - r0 h)2s2 h4 ≈ -2s0 s2 h4 , 42 M15 = s 0 q2 (86) (87) (88) h7 3 . (89) С учетом (85)-(89) вычислим оценку второй компоненты в (67): k2 =- δghd k-3 k k2 k k2 k - r0 R0k-2 ) + dk2 dk2 12 R1 + d13 R3 + d14 R4 11 (s0 R0 ≈ dk2 15 3s2 (s0 h3 R0k-3 - r0 h3 R0k-2 + 6s0 R1k + 12s0 R3k - 4s0 R4k ) = 2s0 q2 h3 = O(hk-2 ) + O(hk-1 ) + O(hk-2 ) + O(hk-2 ) + O(hk-2 ) = O(hk-2 ). (90) ≈- Из (63), (67), (84), (90) следует оценка δghk = O(hk-2 ) (91) для задачи (60). Подстановка (77) или (78), (91) в (64) дает порядок аппроксимации задачи (52), (60) меньшим на две единицы степени используемого многочлена Тейлора k независимо от ее четности. Объяснением этого является то, что для пар индексов (1, 1), (1, 4) и (1, 2), (1, 3) в оценках (80)-(84) соответственно и в оценках (86)-(89) отсутствуют какие-либо закономерности на модули слагаемых, тогда как имеющаяся закономерность для перечисленных пар индексов (первые слагаемые различаются лишь знаками, вторые совпадают) в оценках (70)-(73) привела к зависимости порядка аппроксимации от четности используемого многочлена Тейлора k при исследовании задачи (52). Отметим, что при использовании четырехточечного шаблона отсутствие какихлибо закономерностей в оценках (38)-(40) и в оценках (41)-(43) для перечисленных пар индексов привело к порядку аппроксимации задачи, меньшему 176 Численное интегрирование матричным методом краевых задач. . . на две единицы степени используемого многочлена Тейлора k независимо от ее четности. 6. Оценка погрешностей. При выполнении численного эксперимента использованы следующие нормы: - в качестве суммарной оценки относительной погрешности Dxk = n-1 2 i=1 (xi - [xi ]) n-1 i=1 [xi ] · 100 %. (92) Оценку можно трактовать как некий аналог коэффициента вариации в статистике, который характеризует меру разброса в процентах [16]; - в качестве оценки абсолютной погрешности [5, 6] Exk = max xi - [xi ] , i = 1, 2, . . . , n - 1. (93) В качестве примера рассматривалось ОДУ3 (sin t + x)x + 3(cos t + 1)x - 3 sin t · x - cos t · x = sin t, t ∈ [7, 11] (94) с граничными условиями x(7) = 8.5211, x (7) = 0.2236, x(11) = 14.5995. (95) Общее решение ОДУ3 (94) заимствовано из [17]. В расчетах принято, что n = 20, h = 0.20. Результаты численного эксперимента для краевой задачи (94), (95) приведены в табл. 1-3. В табл. 1 приведены результаты, когда при построении СЛАУ был использован смешанный вариант ПП; в табл. 2 - результаты, когда были использованы только левые шаблоны. Как видно, результаты, представленные в табл. 1, 2, сравнимы между собой, что подтверждает теоретический вывод о том, что выбор того или иного четырехточечного шаблона не оказывает влияния на порядок аппроксимации. В табл. 1 нормы Dxk , Exk для производных x (t) характеризуют оценки относительной и абсолютной погрешностей, вычисленных по формулам (92), (93), в которых значения функций заменены на значения своих первых производных, найденных по формулам (21), (22) при l = 2 во внутренних узлах области интегрирования. Результаты численного эксперимента при использовании пятиточечного шаблона для краевой задачи (94), (95) приведены в табл. 3. В табл. 3 нормы Dxk , Exk для производных x (t) характеризуют оценки относительной и абсолютной погрешностей, вычислены по формулам (92), (93). В них значения функций заменены на значения своих первых производных, найденных по формулам (49), в которых i = 2, 3, . . . , n - 2. Анализ таблиц свидетельствует, что для рассматриваемой задачи для четырехточечного и пятиточечного шаблонов с увеличением степени k используемого многочлена Тейлора относительная и абсолютная погрешности уменьшаются, как это имело место при исследовании краевых задач для ОДУ2 [1, 12, 13] и систем ОДУ2 [18]. 177 178 k is the degree of the Taylor polynomial (P k ) used for the approximation Таблица 3 Значения погрешностей для решения краевой задачи (94), (95) при использовании пятиточечного шаблона [Error estimation in the numerical solution of the boundary value problem (94), (95) when the five-point difference scheme was used] k 4 5 6 7 8 9 10 для решения (for the solution) Dxk , % 5.02 · 10-3 6.56 · 10-3 9.44 · 10-5 2.13 · 10-5 6.34 · 10-7 1.72 · 10-7 7.16 · 10-9 k -3 -3 -5 -5 -7 -7 Ex 5.06 · 10 5.19 · 10 7.85 · 10 2.03 · 10 6.21 · 10 1.87 · 10 6.12 · 10-9 для первой производной решения (for the first derivative of the solution) k -2 Dx , % 4.43 · 10 4.34 · 10-2 7.48 · 10-4 2.51 · 10-4 6.94 · 10-6 1.76 · 10-6 5.04 · 10-8 Exk 5.98 · 10-3 4.72 · 10-3 9.70 · 10-5 3.88 · 10-5 1.07 · 10-6 2.90 · 10-7 7.43 · 10-9 Таблица 2 Значения погрешностей для решения краевой задачи (94), (95) при использовании левых шаблонов в процедурах построения СЛАУ [Error estimation in the numerical solution of the boundary value problem (94), (95) when a left version of the template for the four-point difference scheme was used] k 3 4 5 6 7 8 9 10 Dxk , % 1.10 · 10-1 4.65 · 10-3 1.05 · 10-3 1.04 · 10-4 9.51 · 10-6 6.51 · 10-7 6.03 · 10-8 4.38 ·10-9 Exk 8.28 · 10-2 4.75 · 10-3 9.15 · 10-4 8.54 · 10-5 9.07 · 10-6 6.25 · 10-7 6.75 · 10-8 4.84 · 10-9 Таблица 1 Значения погрешностей для решения краевой задачи (94), (95) при использовании смешанного варианта процедур построения СЛАУ [Error estimation in the numerical solution of the boundary value problem (94), (95) when a mixed version of the template for the four-point difference scheme was used] k 3 4 5 6 7 8 9 10 для решения (for the solution) Dxk , % 1.02 · 10-1 5.63 · 10-3 1.30 · 10-3 8.92 · 10-5 1.06 · 10-5 6.29 · 10-7 6.36 · 10-8 4.35 · 10-9 -2 -3 -3 -5 -6 -7 -8 k 7.88 · 10 5.55 · 10 1.08 · 10 7.49 · 10 9.96 · 10 6.05 · 10 6.96 · 10 4.75 · 10-9 Ex для первой производной решения (for the first derivative of the solution) Dxk , % 6.51 · 10-1 6.05 · 10-2 9.41 · 10-3 7.40 · 10-4 9.53 · 10-5 7.07 · 10-6 8.51 · 10-7 7.51 · 10-8 -2 -2 -3 -5 -5 -6 -7 k Ex 7.89 · 10 1.17 · 10 1.33 · 10 9.72 · 10 1.43 · 10 1.07 · 10 1.12 · 10 1.17 · 10-8 М а к л а к о в В. Н., С т е л ь м а х Я. Г. Численное интегрирование матричным методом краевых задач. . . Выводы. Сформулируем основные выводы по работе. 1. При использовании четырехточечного шаблона теоретически выявлены закономерности между порядком аппроксимации матричного метода и степенью используемого многочлена Тейлора при численном интегрировании краевых задач для линейных неоднородных ОДУ третьего порядка. Установлено, что порядок аппроксимации пропорционален степени используемого многочлена Тейлора и меньше нее на две единицы. 2. Теоретический вывод о том, что выбор того или иного четырехточечного шаблона не оказывает влияния на порядок аппроксимации, подтвержден численным экспериментом. 3. При использовании пятиточечного шаблона предложена процедура построения фиктивного граничного условия, позволяющая построить замкнутую систему разностных уравнений матричного метода численного интегрирования. 4. При использовании пятиточечного шаблона система разностных уравнений разбита на две подсистемы (подзадачи): в первую подсистему вошли два уравнения, первое из которых содержит заданное значение производной в граничных условиях задачи, второе - вычисленное из фиктивного граничного условия значение; во вторую подсистему вошли оставшиеся разностные уравнения построенной замкнутой системы. Вычислена невязка и дана оценка порядка аппроксимации метода в зависимости от выбранной степени многочлена Тейлора. Теоретически выявлены закономерности между порядком аппроксимации матричного метода и степенью используемого многочлена Тейлора. Установлено следующее: а) порядок аппроксимации первой подсистемы, второй подсистемы при четном значении степени используемого многочлена Тейлора и всей задачи пропорционален этой степени и меньше нее на две единицы; б) порядок аппроксимации второй подсистемы при нечетном значении степени используемого многочлена Тейлора пропорционален этой степени и меньше нее на единицу. Конкурирующие интересы. Заявляем, что в отношении авторства и публикации этой статьи конфликта интересов не имеем. Авторский вклад и ответственность. Все авторы принимали участие в разработке концепции статьи и в написании рукописи. Авторы несут полную ответственность за предоставление окончательной рукописи в печать. Окончательная версия рукописи была одобрена всеми авторами. Финансирование. Исследование выполнялось без финансирования.

About the authors

Vladimir N Maklakov

Samara State Technical University

Email: makvo63@yandex.ru
Cand. Phys. & Math. Sci.; Associate Professor; Dept. of Hight Mathematics & Applied Computer Science 244, Molodogvardeyskaya st., Samara, 443100, Russian Federation

Yanina G Stelmakh

Samara State Technical University

Email: yaninastelmah@rambler.ru
Cand. Educat. Sci.; Associate Professor; Dept. of Hight Mathematics & Applied Computer Science 244, Molodogvardeyskaya st., Samara, 443100, Russian Federation

References

Supplementary files