Numerical integration of the boundary value problems for the second order nonlinear ordinary differential equations of an arbitrary structure using an iterative procedure

Full Text

Численное интегрирование краевых задач для нелинейных обыкновенных дифференциальных. . . На основании предложенного в работе [1] метода интегрирования краевых задач для линейных обыкновенных дифференциальных уравнений второго порядка (ОДУ2) в [2] рассмотрена итерационная процедура нахождения численного решения краевых задач для нелинейных ОДУ2 специального вида r(x, u, u , u )u + p(x, u, u , u )u + q(x, u, u , u )u = f (x, u, u , u ), (1) где r(x, u, u , u ), p(x, u, u , u ), (x, u, u , u ), f (x, u, u , u ) - заданные функции, дифференцируемые по своим аргументам нужное число раз. Поставим целью построить итерационную процедуру численного интегрирования краевых задач для нелинейных ОДУ2 произвольной структуры. 1. Построение итерационной процедуры. Исследуем краевую дифференциальную задачу для ОДУ2 произвольной структуры Ψ(x, u, u , u ) = 0, x ∈ [a, b], u(a) = u0 , u(b) = un , (2) где отрезок [a, b] - область интегрирования D, Ψ(x, u, u , u ) - заданная функция, дифференцируемая по своим аргументам нужное число раз, u0 , un - заданные числа. На области D зададим сетку Dh равноотстоящими узлами, определяемыми значениями xi = x0 + ih, i = 1, 2, . . . , n - 1, x0 = a, xn = b, h = (b - a)/n. Очевидно, ОДУ2 в краевой задаче (2) и следующее ОДУ2 u + pu + qu = u + pu + qu + λΨ(x, u, u , u ), где p, q, λ - постоянные величины, требующие дальнейшего уточнения, эквивалентны в смысле совпадения решений. Поэтому далее вместо (2) будем рассматривать краевую задачу вида u + pu + qu = F (x, u, u , u ), u(a) = u0 , u(b) = un , (3) где F (x, u, u , u , λ) = u + pu + qu + λΨ(x, u, u , u ) - известная функция при заданных значениях p, q, λ. В дальнейшем для краткости примем для любой функции обозначение ϕ(xi ) = ϕi , где xi - узел сетки Dh . При некотором фиксированном натуральном числе k 2 для узла сетки ti , i = 1, 2, . . . , n - 1, составим систему уравнений, состоящую из (k + 1)-го уравнения. В качестве первого уравнения системы примем первые k + 1 члена разложения в ряд Тейлора функции u в узле xi-1 , в качестве второго - в узле xi+1 . Дополнительно в систему внесём равенства, полученные дифференцированием по аргументу x обеих частей уравнения в задаче (3), так что (r) (qui + pui + ui )(r) = Fi , где r = 0, 1, . . . , k - 2, которые запишем в узле xi . Получим  k h2 h3  k h u(k) = u   u - hu + - + · · · + (-1) u u i-1 , i  i  2! i 3! i k! i   2 3 k   u + hu + h u + h u + · · · + h u(k) = u ,   i+1 i i  2! i 3! i k! i (4) qui + pui + ui = Fi ,    qui + pui + ui = Fi ,     ..   .     (k-2) (k-1) (k) (k-2) qui + pui + ui = Fi . 355 М а к л а к о в В. Н. В матричной форме система уравнений (4) имеет вид Ak Vik = Gki в обозначениях  1 -h   1 h  k A = q p  0 q .  .. ... 0 0 h3 h2 - 2! 3! h2 h3 2! 3! 1 0 p 1 .. .. . . 0 0 ... ... ... ... .. . ...  hk k (-1) hk k! 0 0 .. . 1     ui ui-1 k!     ui   ui+1   u   F  i      i  , Vik =  u  , Gki =  F  . i    i     ..   ..    .   .   (k-2) (k)  Fi ui Элементы матрицы Ak могут быть вычислены, следовательно, (4) является системой линейных алгебраических уравнений. Предполагая существование обратной матрицы (Ak )-1 от матрицы Ak , найдём Vik = (Ak )-1 Gki , или в координатной форме: (k-2) , (5) (k-2) , (6) ui = ak11 ui-1 + ak12 ui+1 + ak13 Fi + ak14 Fi + · · · + ak1, k+1 Fi ui = ak21 ui-1 + ak22 ui+1 + ak23 Fi + ak24 Fi + · · · + ak2, k+1 Fi (k-2) ui = ak31 ui-1 + ak32 ui+1 + ak33 Fi + ak34 Fi + · · · + ak3, k+1 Fi (7) , ... (k) ui (k-2) = akk+1, 1 ui-1 + akk+1, 2 ui+1 + akk+1, 3 Fi + · · · + akk+1, k+1 Fi , (8) где под aklm , l, m = 1, 2, . . . , k + 1, понимаем соответствующие элементы матрицы (Ak )-1 . Составим систему уравнений, в которую внесём равенства (5) для всех i = 1, 2, . . . , n - 1:  k-2   (m) k   -u + a u = - ak1, m+3 F1 - ak11 u0 , 1 12 2    m=0    k-2   (m)  k k  a11 u1 - u2 + a12 u3 = - ak1, m+3 F2 ,     m=0 k-2 (9) (m) k k   a u - u + a u = - ak1, m+3 F3 , 2 3 4 11 12    m=0     . . .    k-2    (m) k  ak1, m+3 Fn-1 - ak12 un .  a11 un-2 - un-1 = - m=0 Отметим, что (9) есть система разностных уравнений с постоянными коэффициентами для трёхточечного шаблона xi-1 , xi , xi+1 , i = 1, 2, . . . , n - 1, 356 Численное интегрирование краевых задач для нелинейных обыкновенных дифференциальных. . . при построении которой не использовались конечные разности для аппроксимации производных, как это было сделано при численном интегрировании краевых задач для обыкновенных дифференциальных уравнений в [3-6] и для ряда дифференциальных уравнений в частных производных в [5-8]. В матричной форме система (9) имеет вид B k U = Gk в обозначениях -1 ak12 0 0 ak11 -1 ak12 0  0 ak11 -1 ak12 Bk =   . .. .. ..  . . . . . 0 0 0 0 ... ... ... .. . ...    0 0  0 , ..   . 1   u1  u2    u3  , U =  .   ..  un-1  k-2 (m) ak1, m+3 F1 - - ak11 u0      m=0   k-2   (m) k   - a1, m+3 F2     m=0   k-2 k k (1) (k)  . G = G (U, U , . . . , U ) =  (m) k  - a F   1, m+3 3   m=0   .   ..    k-2    (m) - k k a1, m+3 Fn-1 - a12 un  m=0 Bk В предположении существования обратной матрицы (B k )-1 от матрицы найдём U = (B k )-1 Gk , или в развёрнутом виде U = (B k )-1 Gk (U, U (1) , . . . , U (k) ), (r) (r) (10) (r) где U (r) = u1 u2 · · · un-1 - обозначение вектор-столбца производных степени r функции u во внутренних узлах сетки Dh . Матричное уравнение (10) наряду с (6)-(8) будем использовать при построении итерационной процедуры для нахождения приближённого решения задачи (3). В качестве начального (нулевого) приближения U0 примем расположенные на прямой u = cx + d, проведённой через точки u(a) = u0 , u(b) = un , значения функции u во внутренних узла сетки. В этом случае, очевидно, значения производных нулевого приближения определим как (1) U0 = c c ··· c , (r) U0 = 0 0 ··· 0 , r = 2, 3, . . . , k. (r) Подставим начальные приближения U0 , U0 , r = 1, 2, . . . , k, в правые части матричного равенства (10), их вычислим и найденные значения примем в качестве первого приближения: (1) (k) U1 = (B k )-1 Gk (U0 , U0 , . . . , U0 ). 357 М а к л а к о в В. Н. (r) Компоненты векторов первого приближения U1 , r = 1, 2, . . . , k, вычислим по формулам (6)-(8) соответственно, в правые части которых подставим найденные значения вектора первого приближения U1 и векторов начального (r) приближения U0 , r = 1, 2, . . . , k. (r) Найденные векторы первого приближения U1 , U1 , r = 1, 2, . . . , k, вновь подставим в правую часть равенства (10), их вычислим и полученные значения примем в качестве второго приближения U2 для искомой функции u. (r) Для вычисления U2 , r = 1, 2, . . . , k, воспользуемся формулами (6)-(8). Описанная итерационная процедура может быть выполнена произвольное число раз при условии её сходимости по следующим формулам (1) (k) Ut = (B k )-1 Gk (Ut-1 , Ut-1 , . . . , Ut-1 ), (11) (r) uit = ak(r+1) 1 u(i-1) t + ak(r+1) 2 u(i+1) t + k-2 (m) ak(r+1) (m+3) Fi + i = 1, 2, . . . , n - 1, r = 1, 2, . . . , k, (12) , m=0 где t - номер итерации, uit - компоненты вектора Ut . 2. Выбор постоянных. Постоянные величины p, q выберем из требования обеспечения устойчивости задачи на каждой итерации. Разностную краевую задачу с переменными коэффициентами ai ui-1 + bi ui + ci ui+1 = fi , u(a) = u0 , u(b) = un , i = 1, 2, . . . , n - 1, (13) называют устойчивой [5], если она имеет единственное решение при произвольных u0 , un , f1 , f2 , . . . , fn-1 и имеет место оценка для норм u (14) C f , где C - не зависящее от h число, а в качестве нормы некоторого вектора w = w0 w1 · · · wn принято (15) . w = max w0 , w1 , . . . , wn В соответствии с (15) в неравенстве (14) положено u = u0 u1 · · · un , f = u0 un f1 f2 · · · fn-1 . Правомерность выбора нормы в виде (15) для разностной краевой задачи (13) обоснована в [5]. Разностную краевую задачу с постоянными коэффициентами aui-1 + bui + cui+1 = fi , u(a) = u0 , i = 1, 2, . . . , n - 1 u(b) = un , (16) называют хорошо обусловленной [5], если она при всех достаточно больших n имеет единственное решение при произвольных u0 , un , f1 , f2 , . . . , fn-1 и если значения u0 , u1 , . . . , un , образующие решение, удовлетворяют неравенству ui 358 P max u0 , un , f1 , f2 , . . . , fn-1 , i = 0, 1, . . . , n, (17) Численное интегрирование краевых задач для нелинейных обыкновенных дифференциальных. . . где P - не зависящее от h число. Выполнение неравенства (17) означает, что погрешности, допущенные при задании u0 , un , f1 , f2 , . . . , fn-1 , не накапливаются и не приводят к возрастающим с увеличением n ошибкам [5]. Неравенство (14) и совокупность неравенств (17) эквивалентны в смысле совпадения оценок ввиду выбора нормы в форме (15) и в силу совпадения их правых частей; следовательно, хорошая обусловленность влечёт устойчивость на каждой итерации для рассматриваемой задачи. Достаточным условием (критерием) хорошей обусловленности задачи (16) является выполнение неравенства [5] b - a+c a + b + c или Θ>0 b - a + c > 0. (18) Критерий (18) для задачи (9) имеет вид 1 - ak11 + ak12 > 0 или k k det Ak - M11 + M12 > 0, (19) k , M k - алгебраические дополгде det Ak - определитель матрицы Ak ; M11 12 нения элементов, расположенных на пересечениях первой строки с первым и вторым столбцами, соответственно, матрицы (Ak ) . Критерий (19) для задачи (9) при k = 2 будет записан как 2h - qh3 -2h > 0; откуда следуют решения q < 0 и q > 4/h2 . Отбросим второе решение как зависящее от величины шага h. Критерий (19) при k = 3 примет вид 2h + откуда (p2 - 4q)h3 q 2 h5 (p2 - q)h3 + - 2h + > 0, 3 6 3 6 + (p2 - q)h2 > 0, q(qh2 - 6) > 0, 6 + (p2 - q)h2 < 0, q(qh2 - 6) < 0. Решением первой системы является q < 0 для любого значения p, вторая система несовместна. Нахождение аналитического решения неравенства (19) при k > 3 оказалось довольно трудоёмким. В силу того, что элементы обратной матрицы (Ak )-1 зависят только от параметров p и q, с использованием персонального компьютера (ПК) была смоделирована процедура вычисления разности k + M k , k = 4, 5, . . . , 8, для различных p и q. Для нескольких det Ak - M11 12 произвольных значений p неравенство (19) обращалось в истину для ряда отрицательных значений q, в частности и для q = -1, которое и было, наряду с p = 1, использовано при численном эксперименте, что гарантированно обеспечивало устойчивость на каждой итерации. 359 М а к л а к о в В. Н. 3. Численный эксперимент. При выполнении численного эксперимента были использованы ОДУ2 специального вида (1) (20) uu - u u = 1, имеющее аналитическое решение ua = 1 ch(C1 x + C2 ), C1 c краевыми условиями u(-3) = 5.10, (21) u(3) = 7.71, и ОДУ2 произвольной структуры (22) eu - (u )2 - x = 0, имеющее аналитическое решение в параметрической форме ua (s) = 41 e2s (2s - 3) + 13 es (3C1 - 2s3 + 6s - 6) + x(s) = es - s2 , 1 5 15 4s - C1 s2 + C2 , с краевыми условиями u(0) = 1.11, u(6) = 7.43. (23) Общее решение ОДУ2 (20) и метод нахождения общего решения ОДУ2 (22) заимствованы из [9]. При выполнении численного эксперимента был использован метод перебора [6], или метод сканирования [10], согласно которому между некоторыми выбранными начальным λstr и конечным λfin значениями с периодом дискретизации λfin - λstr ∆λ = r вычислялись промежуточные значения λm = λstr + m∆λ, m = 0, 1, . . . , r и с каждым из них была реализована итерационная процедура (11), (12) для краевых задач (20), (21) и (22), (23). Для каждого λm при всех k = 2, 3, . . . , 8 на каждой итерации c номером t вычислялись следующие оценки погрешности на всём отрезке интегрирования: 1) D(k, λ, t) = n-1 i=1 (uit - n-1 a i=1 ui uai )2 · 100% - суммарная оценка относительной погрешности на итерации с номером t, которую можно трактовать как некий аналог коэффициента вариации в статистике, характеризующий меру разброса в процентах [11]; величина D отличается от коэффициента вариации тем, что стандартное отклонение заменено корнем квадратным из остаточной дисперсии [11]; 360 Численное интегрирование краевых задач для нелинейных обыкновенных дифференциальных. . . a 2) E(k, λ, t) = n-1 i=1 uit - ui - суммарная оценка абсолютной погрешности; 3) H(k, λ, t) = max uit - uai , i = 1, 2, . . . , n - 1 - максимальная оценка абсолютной погрешности; 4) Q(k, λ, t) = n-1 2 i=1 (uit - ui (t-1) ) n-1 i=1 ui (t-1) · 100% - суммарная оценка относительной погрешности между приближёнными решениями на текущей и предыдущей итерациями с номерами t и (t - 1) соответственно; 5) R(k, λ, t) = n-1 i=1 uit - ui (t-1) - суммарная оценка абсолютной погрешности между текущей и предыдущей итерациями; 6) S(k, λ, t) = max uit - ui (t-1) , i = 1, 2, . . . , n - 1 - максимальная оценка абсолютной погрешности между текущей и предыдущей итерациями. Приведём некоторые результаты, выявленные сканированием по параметру λ в задачах (20), (21) и (22), (23). Расчётным путём для рассматриваемых задач установлено, что при фиксированных k и λ оценки погрешностей D, E, H между точным и приближённым решениями с увеличением номера итерации t изменялись практически синхронно (росли или убывали одновременно). Аналогичная ситуация складывается с оценками погрешностей Q, R, S между приближёнными решениями на текущей и предыдущей итерациями. Для каждого набора чисел k и λ при отсутствии переполнения ПК реализация итерационной процедуры прекращалась на итерации с номером tust при выполнения условия D(k, λ, tust ) - D(k, λ, tust-1 ) 10-15 . (24) В расчётах для обеих задач принято n = 20, h = 0.3. В итоге предварительного грубого сканирования области λ ∈ [-10; 10] с периодом дискретизации ∆λ = 5 · 10-2 для обеих задач были найдены интервалы изменения параметра λ, в которых выполнено условие (24), а далее было положено ∆λ = = 4.5·10-3 в задаче (20), (21) и ∆λ = 5.5·10-3 в задаче (22), (23) для итогового сканирования. Для дискретно заданных функций D(t), Q(t) численным экспериментом установлено, что в интервале изменения параметра λ для всех k = 2, 3, . . . , 8: а) график функции D(t) имеет один локальный минимум в точке tmin и одну точку перегиба; б) график функции D(t) имеет горизонтальную асимптоту D = Dust = = D(tust ); в) точки tmin , tust смещаются вправо при уменьшении λ ; г) при фиксированном k погрешность Dmin = D(tmin ) практически не зависит от λ; д) функция Q(t) оказалась невозрастающей, её график имеет горизонтальную асимптоту Q = 0; е) при фиксированном k модуль производной Q(t) уменьшается при уменьшении λ . Результаты расчётов для задачи (20), (21) приведены в табл. 1, для задачи (22), (23) - в табл. 2. В табл. 1, 2 принято Qmin = Q(tmin ), Qust = Q(tust ); в качестве tmin указано наименьшее значение из всех вычисленных tmin (λ) 361 М а к л а к о в В. Н. Таблица 1 Значения погрешностей краевой задачи (20), (21) [The values of the errors for the boundary value problem (20), (21)] k λ tmin Dmin , % Qmin , % tust Dust , % Qust , % 2 3 4 5 6 7 8 -0.3240 -0.2970 -0.2790 -0.2745 -0.2700 -0.2700 -0.2700 34 39 51 53 49 49 49 9.74 · 10-4 3.42 · 10-4 6.94 · 10-7 6.79 · 10-7 6.74 · 10-7 6.74 · 10-7 6.74 · 10-7 2.56 · 10-4 1.33 · 10-4 1.69 · 10-7 1.93 · 10-7 7.47 · 10-8 7.49 · 10-8 7.46 · 10-8 37 42 76 83 87 87 88 1.19 · 10-3 3.73 · 10-4 7.19 · 10-7 6.82 · 10-7 6.94 · 10-7 6.94 · 10-7 6.94 · 10-7 2 3 4 5 6 7 8 -0.3960 -0.3740 -0.3630 -0.3575 -0.3465 -0.3465 -0.3465 29 30 67 71 69 69 69 6.13 · 10-2 9.78 · 10-3 9.69 · 10-5 9.71 · 10-5 8.18 · 10-6 8.18 · 10-6 8.17 · 10-6 1.71 · 10-2 1.92 · 10-2 6.78 · 10-6 3.64 · 10-6 2.89 · 10-6 3.53 · 10-6 4.20 · 10-6 66 79 98 102 105 111 113 6.22 · 10-2 6.22 · 10-2 9.78 · 10-5 9.78 · 10-5 1.11 · 10-5 1.11 · 10-5 1.15 · 10-5 2.49 · 10-5 7.21 · 10-6 9.94 · 10-9 7.14 · 10-10 7.24 · 10-11 4.48 · 10-11 2.32 · 10-11 Таблица 2 Значения погрешностей краевой задачи (22), (23) [The values of the errors for the boundary value problem (22), (23)] k λ tmin Dmin , % Qmin , % tust Dust , % Qust , % 4.19 · 10-6 3.71 · 10-7 1.12 · 10-9 3.45 · 10-10 2.54 · 10-10 1.61 · 10-10 6.97 · 10-11 при фиксированном k, что соответствует левой границе интервала изменения λ. Правая граница упомянутого интервала оказалась меньше нуля на 3-4 периода дискретизации ∆λ для обеих задач. Данные табл. 1, 2 свидетельствуют о том, что, как и в работах [1, 2, 12], имеет место уменьшение погрешностей с увеличением числа k. Анализ итогов численного эксперимента позволяет сделать вывод, что параметр λ, как оказалось, являющийся «ответственным» за сходимость итерационной процедуры, можно трактовать как некую характеристику «противовеса», уравнивающую «вес» влияния слагаемых u + pu + qu и λΨ(x, u, u , u ) друг на друга в правой части задачи (3). Теоретическое обоснование выбора значения параметра λ требует дополнительного исследования. Также требуют дальнейшей теоретической проработки вопросы, связанные с определением момента достижения заданной точности. 4. Выводы. Предложена итерационная процедура численного интегрирования краевых задач для обыкновенных дифференциальных уравнений второго порядка произвольной структуры. Опытным путём выявлены некоторые закономерности влияния значения параметра λ на сходимость. Теоретически установлены значения коэффициентов, обеспечивающие устойчивость на каждой итерации независимо от вида исходного уравнения произвольной структуры. Численным экспериментом показано, что использование многочленов Тейлора произвольных степеней вместо конечных разностей при аппроксимации производных приводит к уменьшению погрешностей между точным и найденным численным решениями.

About the authors

Vladimir N Maklakov

Samara State Technical University

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

References

Supplementary files