Scientific journal
Fundamental research
ISSN 1812-7339
"Перечень" ВАК
ИФ РИНЦ = 1,749

PIECEWISE INTERPOLATION SOLUTION OF THE TWO-POINT CAUCHY PROBLEM FOR ORDINARY DIFFERENTIAL EQUATIONS WITH ITERATIVE REFINEMENT

Romm Ya.E. 1 Dzhanunts G.A. 1
1 Taganrog Institute named after A.P. Chekhov (branch) Rostov State University of Economics
The method of the variable piecewise interpolation solution of the one-point and two-point Cauchy problem for ordinary differential equations is based on the system of subintervals in the adjacent borders of which are taken equal to the interpolation node values. On the basis of the primitive Newton polynomial interpolation, built to the right of the function and presented in the form of an algebraic polynomial with numerical coefficients, the iterative refinement is performed according to the type of successive approximations Picard. The resulting approximate solution is continuous and continuously differentiable on the entire interval approached. The method converges uniformly to the exact solution with increasing number of subintervals simultaneously performed evenly approximation of the derivative. Program realization enables the accuracy of the approximation, significantly surpassing the accuracy of difference methods. Approximate solution of the two-point Cauchy problem has the same properties and the correct value at the end point.
ordinary differential equations
one-point and two-point Cauchy problem
piecewise interpolation solution
interative refinement
1. Afanasev A.P., Dzjuba S.M., Kirichenko M.A., Rubanov N.A. Priblizhennoe analiticheskoe reshenie sistem obyknovennyh differencialnyh uravnenij s polinomialnoj pravoj chastju // Zhurnal vychislitelnoj matematiki i matematicheskoj fiziki, 2013, vol. 53, no. 2, pp. 321–328.
2. Afanasev A.P., Dzjuba S.M. Metod postroenija priblizhennyh analiticheskih reshenij differencialnyh uravnenij s polinomialnoj pravoj chastju // Zhurnal vychislitelnoj matematiki i matematicheskoj fiziki, 2015, vol. 55, no. 10, pp. 1694–1702.
3. Berezin I.S., Zhidkov N.P. Metody vychislenij, vol. 1, Moscow, Nauka Pupl., 1966. 632 p.
4. Dzjadyk V.K. O primenenii linejnyh metodov k priblizheniju polinomami reshenij obyknovennyh differencialnyh uravnenij i integralnyh uravnenij Gammershtejna // Izvestija AN SSSR. Ser. matem, 1970, vol. 34 (4), pp. 827 – 848.
5. Kasti Dzh., Kalaba R. Metody pogruzhenija v prikladnoj matematike. Monograph. Moscow, Mir Publ., 1976. 224 p.
6. Rihtmajer R., Morton K. Raznostnye metody reshenija kraevyh zadach. Moscow, Mir Publ., 1972. 420 p.
7. Romm Ya.E. Beskonfliktnye i ustojchivye metody determinirovannoj parallelnoj obrabotki: Avtoreferat dissertacii na soiskanie uchenoj stepeni doktora tehnicheskih nauk. Taganrog: TRTU, 1998. 42 p.
8. Romm Ya.E. Lokalizacija i ustojchivoe vychislenie nulej mnogochlena na osnove sortirovki. II // Kibernetika i sistemnyj analiz. 2007, no. 2, pp. 161–174.
9. Romm Ya.E., Dzhanunts G.A. The computer method of variable piecewise polynomial approximation of functions and solutions of ordinary differential equations (2013) Cybernetics and Systems Analysis, 49 (3), pp. 409–423.
10. Romm Ya.E., Dzhanunts G.A. Kompjuternoe kusochno-interpoljacionnoe reshenie od-notochechnoj i dvuhtochechnoj zadachi Koshi dlja obyknovennyh differencialnyh uravnenij / Dep. manuscript No. no. 47-В2016 05.04.2016.
11. Awoyemi D.O., Kayode S.J., Adoghe L.O. A Five-Step P-Stable Method for the Numerical Integration of Third Order Ordinary Differential Equations // American Journal of Computational Mathematics, no. 4, 2014, pp. 119–126.
12. Fatimah B.O., Senapon W.A., Adebowale A.M. Solving Ordinary Differential Equations with Evolutionary Algorithms // Open Journal of Optimization, no. 4, 2015, pp. 69–73, available at: http://dx.doi.org/10.4236/ojop.2015.43009.

Численное решение нелинейной двухточечной граничной задачи актуально для разделов оптимального управления, вариационных задач, астрофизики, биологии [5]. К основным методам решения относится сведение двухточечных задач к одноточечным задачам Коши с помощью конечно-разностных схем [6]. Излагаемый ниже метод отличается от известных по построению и тем, что дает непрерывное и непрерывно дифференцируемое компьютерное приближение решения двухточечной задачи. Метод основан на варьируемом кусочно-интерполяционном приближении решения задачи Коши для обыкновенных дифференциальных уравнений (ОДУ) с итерационным уточнением [9, 10]. Теоретические аналоги содержатся в [1, 2, 4], сходные эксперименты описаны в [11, 12]. В отличие от аналогов, предлагаемое итерационное уточнение строится на каждом подынтервале с помощью первообразной от интерполяционного полинома Ньютона, аппроксимирующего функцию правой части. Интерполяционный полином на подынтервале преобразуется к виду алгебраического полинома с числовыми коэффициентами, что обеспечивает табличную первообразную и позволяет программировать последовательность интегральных приближений Пикара. Для решения двухточечной задачи данный метод строится одновременно от начальной точки по направлению конечной и от конечной – к начальной. В окрестности точки пересечения встречные приближения интерполируются с применением итерационного уточнения. В результате в один проход получается непрерывное и непрерывно дифференцируемое кусочно-интерполяционное решение двухточечной задачи с точными значениями в начале и в конце промежутка.

Пусть вначале рассматривается задача Коши для ОДУ,

y′= f(x, y); y(x0) = y0, (1)

предполагается, что в области

romm01.wmf

функция f(x, y) определена, непрерывна, непрерывно дифференцируема (в точках x0 – справа, xfin – слева), удовлетворяет условию Липшица: romm02.wmf romm03.wmf. Предполагается также, что решение задачи (1) существует и единственно во всей R.

Варьируемое кусочно-интерполяционное приближение

Пусть F(x) – действительная функция одной действительной переменной, x ∈ [α, β], ε > 0 – априори заданная граница абсолютной погрешности. Строится объединение подынтервалов равной длины:

romm04.wmf P = 2k, k = 0, 1, ... (2)

На каждом подынтервале строится интерполирующий F(x) полином Ньютона Ψin(t) степени n с отстоящими на romm05.wmf узлами, n выбирается одинаковым для всех подынтервалов и минимальным при условии

romm06.wmf x ∈ [xri, xri+1], romm07.wmf romm08.wmf (3)

Интерполяционный полином

romm09.wmf romm10.wmf (4)

где xij = xi + jh – узлы интерполяции, romm11.wmf; ΔjFi0 – конечная разность j-го порядка, преобразуется к виду romm12.wmf romm13.wmf. Для этого вычисляются romm14.wmf romm15.wmf romm16.wmf, полагается bij = ΔjFi0, romm17.wmf. Полином romm18.wmf разложен по корням k = 0, 1, ..., j – 1, находятся его коэффициенты:

romm19.wmf (5)

Пусть zl = l, romm20.wmf. Для любых корней zl коэффициенты (5) дает формула [7, 8]

romm21.wmf

где k-й шаг умножения матриц справа налево записывается в виде

romm22.wmf (6)

Здесь

romm23.wmf,

romm24.wmf, romm25.wmf.

При k = j левые части (6) совпадут с коэффициентами (5). После преобразований интерполяционный полином примет вид

romm26.wmf (7)

где ai0 = F(xi0), romm27.wmf Проверка точности приближения (3) полиномом (7) выполняется на каждом подынтервале в точках, равноотстоящих на h/γ, γ ≥ 3 – параметр. Для x ∈ [xi, x i+1) адрес выборки соответственных коэффициентов (7) даст i = [(x – α)/ρ], ρ = x i+1 – xi, romm28.wmf. Построение полинома начинается от n = 1 при k = 0 для каждого i из (2), P = 2k, с проверкой (3) во всех проверочных точках. При нарушении (3) хоть в одной из них значение k увеличивается на единицу, так продолжается до априори заданной границы k ≤ k0. Если в результате искомая точность приближения не достигнута, то полагается k = 0, при этом n увеличивается на единицу, и так – до априори заданной границы n ≤ n0. Фиксируется наименьшее n, при котором (3) верно во всех проверочных точках всех подынтервалов, в соответствии с этим фиксируется k [7, 9].

Подобным способом аппроксимируется f(x, φ(x)) из (1), где y ≈ φ(x). Вначале в правую часть вместо y подставляется y0. Функция f(x, y0) приближается полиномами вида (7) по изложенной схеме так, как если бы F(x) равнялась f(x, y0). При фиксированных на выходе алгоритма n и k на подынтервале [xi–1, xi], i = 1, аналогично, при i = 2, 3, ..., выполняется итерационное уточнение. Первообразная

romm29.wmf

в виде romm30.wmf принимается за приближение решения на данном подынтервале: romm31.wmf. Полагается romm32.wmf, при той же степени n на том же подынтервале строится аппроксимирующий полином (4) в форме (7):

romm33.wmf.

Снова берется первообразная с тем же значением константы

romm34.wmf

в виде romm35.wmf, подставляется в правую часть, romm36.wmf, которая аппроксимируется аналогично:

romm37.wmf

Итерации

romm38.wmf

romm39.wmf l = 1, 2, ...,

где romm40.wmf, продолжаются до достижения заданной малости невязки, при реализации l ≤ q = const, t из (3).

Весь отрезок решения задачи (1) разбивается на отрезки равной длины:

romm41.wmf

α r+1 = βr, romm42.wmf, (8)

каждый отрезок [αr, βr] обрабатывается как [α, β] из (2), аналогично (2) рассматривается

romm43.wmf

r = 0, 1, ..., M – 1, P = 2k, k ∈ {0, 1, ...}. (9)

За начальное значение на левой границе подынтервала [xri, xri+1] всегда принимается конечное значение на правой границе предыдущего подынтервала [xri–1, xri]: romm44.wmf, i = 1, 2, ..., P – 1, аналогично, при переходе от [αr, βr] к [α r+1, βr+1].

Пусть произвольно выбрано n, 1 ≤ n ≤ n0, n0 = const. Предполагается, что f(x, y(x)) имеет n + 1 непрерывную производную в R. Обозначим romm45.wmf, где romm46.wmf – приближение y(x) с помощью полинома romm47.wmf. При P = 1, i = 0, [xr0, xr1] = [αr, βr], погрешность приближения полиномом (4), обозначаемым Ψr0n(t), оценивается из неравенства [3]:

romm48.wmf romm49.wmf (10)

где romm50.wmf romm51.wmf Оценка соответствует k = 0 в (9), что интерпретируется как 0-й этап разбиения αr, βr]. Из (10) следует [9, 10]:

P = 1, [αr, βr] = [xr0, xr1], romm52.wmf romm53.wmf (11)

где 1 ≤ n ≤ n0,

romm54.wmf c = const. (12)

Если P = 2, и на каждом из подынтервалов [xr0, xr1], [xr1, xr2] разбиения (9) выполняется интерполяция полиномом (4) с тем же значением n, то с учетом равенства числа узлов, отстоящих друг от друга на каждом подынтервале на h/2, получится аналог соотношения (11) с вдвое меньшим h romm55.wmf, i = 0, 1. В общем случае [9] –

P = 2k; romm56.wmf romm57.wmf

romm58.wmf

h – шаг интерполирования полинома Ψr0n(t) на [αr, βr] при k = 0. В предположении h < 1,

romm59.wmf romm60.wmf

Пусть romm61.wmf x ∈ [xri, x ri+1]. Для приближения y(x) показано [10], что

romm62.wmf

romm63.wmf (13)

Равномерная по n оценка получится, если учесть, что 1 ≤ n, и выбрать h < 1:

romm64.wmf

romm65.wmf (14)

В рассматриваемых условиях имеет место теорема 1.

Теорема 1 [9, 10]. При любом n, 1 ≤ n ≤ n0, n0 = const, последовательность полиномов Prin+1(x), где i из (9), r из (8), равномерно на [x0, xfin] сходится к решению y(x) задачи (1) при k → ∞. Скорость сходимости оценивается из (13), где c из (12), h – расстояние между узлами полинома Ψr0n(t), интерполирующего правую часть (1) на [αr, βr]. Значение h не зависит от r и не меняется с ростом k при разбиении данного отрезка на 2k подынтервалов. При h < 1 скорость сходимости оценивается как из (13), так и из (14).

В частности, теорема верна в случае двукратной непрерывной дифференцируемости правой части (1).

Следствие 1 [10]. В тех же условиях последовательность полиномов Ψrin(t) равномерно на [x0, xfin] приближает производную y′(x) при k → ∞ с оценками:

romm66.wmf

romm67.wmf

и если h < 1, то, в частности,

romm68.wmf

romm69.wmf

На практике вариации n, 1 ≤ n ≤ n0, и k, 1 ≤ k ≤ k0, реализуются по минимуму невязки так, чтобы дробь 2–k(n+1) оказалась наименьшей за счет произведения k(n+1).

Итерационное уточнение

Вначале на отрезке [αr, βr] выполняются вариации, на основе которых выбираются и фиксируются степень полиномов n, а также количество подынтервалов 2k. Затем на каждом подынтервале [xri, xri+1], i = 0, 1, ..., 2k – 1, выполняется итерационное уточнение:

romm70.wmf

где romm71.wmf x ∈ [xri, x ri+1]; romm72.wmf

romm73.wmf romm74.wmf (15)

romm75.wmf romm76.wmf romm77.wmf

где q – конечный номер итераций на подынтервале

[xri–1, xri], romm78.wmf romm79.wmf romm80.wmf

При l ≥ 2:

romm81.wmf romm82.wmf romm83.wmf (16)

Обозначение romm84.wmf вводится для отличия от romm85.wmf, относительно которого итерационное уточнение не предполагалось. Итерации (15), (16) продолжаются до некоторой фиксированной границы l = q. Абстрактно допускается l → ∞. Показано [10], что если

romm86.wmf (17)

то

romm87.wmf romm88.wmf (18)

где romm89.wmf Если же

romm90.wmf romm91.wmf (19)

то на этом подынтервале

romm92.wmf

romm93.wmf

В условиях (19) итерационное уточнение снижает погрешность кусочно-интерполяционного приближения к решению и его производной пропорционально romm94.wmf при l → ∞, где L – константа Липшица. При условии (17) уменьшение погрешности ckr в (18) имеет коэффициент меньший единицы, если romm95.wmf.

Структура невязки

На отрезке [xj, x] фиксированной длины, x – xj << βr – αr, невязку можно представить в виде

romm96.wmf

где romm97.wmf и i – переменные значения индекса подынтервала, включающего x. Определяется разность аналогичных приближений для текущего и предшествующего количества подынтервалов:

romm98.wmf (20)

где M из (8); i1, i2 – индексы соответственно конечного и начального подынтервалов, которым принадлежат x и xj в случае 2k подынтервалов в (9). По определению, romm99.wmf для 2k подынтервалов, аналогично для – i2, romm100.wmf для 2 k–1 подынтервалов, аналогично – для i4; для случая 2k–1 подынтервалов i3, i4 – номера конечного и начального подынтервалов, соответственно включающих x и xj на отрезке [αr, βr]. Отрезки [xj, x] покрывают [αr, βr]. Для начальной пары фиксированных n и k на объединении (9) находится и запоминается наибольшее ΔM(r, k, n, l), обозначаемое ∇M(r, k, n, l). Поиск ∇M(r, k, n, l) возобновляется при том же n и возрастании k на единицу в пределах romm101.wmf, затем при увеличении n на единицу от начального k, romm102.wmf, и т.д. Среди таких максимумов определяется минимальный по всем k и n:

romm103.wmf

romm104.wmf (21)

В соответствии со значением ∇M(r, k, n, l; min) фиксируется 2k и n в границах [αr, βr]. В программной реализации минимум (21) определяется с итерационным уточнением полиномов в (20) при некотором l = const. Для соответственно зафиксированных k и n заново выполняется окончательное итерационное уточнение. Наибольшая точность и устойчивость метода относительно вида задачи достигается со следующей невязкой [10]:

romm105.wmf.

Здесь ω – постоянное целое, выбранное экспериментально: ω = 5; ωτ – фиксированная часть длины [αr, βr], в случае 1 << βr – αr экспериментально выбранная как romm106.wmf. В случае жестких задач βr – αr < 1, и romm107.wmf p ≥ 2. Знаменатель нормирует невязку, ε > 0 – наименьшее число для исключения деления на ноль. С romm108.wmf выполняется то же, что с ΔM(r, k, n, l), для выбора romm109.wmf, определяемого аналогично ∇M(r, k, n, l; min).

Случай системы ОДУ

Рассматривается задача

Y′ = F(x, Y); Y(x0) = Y0, (22)

где F(x, Y) = (f1(x, Y), f2(x, Y), ..., fN(x, Y)); Y = (y1(x), y2(x), ..., yN(x)); Y0 = (y01, y02, ..., y0N). Здесь и ниже F обозначает вектор-функцию N + 1 переменных в отличие от обозначения функции одной переменной в (3), (4). Используется romm110.wmf Предполагается, что в области romm111.wmf выполнены условия существования и единственности, функция F(x, Y) определена, непрерывна и непрерывно дифференцируема, удовлетворяет условию Липшица:

romm112.wmf romm113.wmf romm114.wmf romm115.wmf

Кусочно-интерполяционное приближение и итерационное уточнение строятся для каждого уравнения системы (22) в полной аналогии случаю одной переменной. Здесь и ниже приближения всех компонентов на шаге рассматриваются при одних и тех же значениях l, k, n одновременно во всех уравнениях системы. Если правая часть (22) (n + 1)-кратно непрерывно дифференцируема в R0, то при каждом romm116.wmf выполнено [10]:

romm117.wmf x ∈ [xr0, xr1],

где Ψjr0n(t) – интерполяционный полином Ньютона в форме с числовыми коэффициентами, построенный для аппроксимации romm118.wmf romm119.wmf romm120.wmf Отсюда следует

romm121.wmf x ∈ [xr0, xr1], (23)

где romm122.wmf C = const,

romm123.wmf

romm124.wmf romm125.wmf (24)

h – шаг интерполирования полинома Ψjr0n(t) на [αr, βr] при k = 0, одинаковый для всех k = 1, 2, ... и для всех romm126.wmf. В предположении h < 1, –

romm127.wmf romm128.wmf (25)

На основе (23)–(25) получаются следующие оценки [10]:

romm129.wmf romm130.wmf

где P jrin+1(x) – первообразная от Ψjr0n(t), приближающая j-ю компоненту решения, при ограничении n ≤ n0, n0 = const,

romm131.wmf

romm132.wmf

и, в предположении h < 1, –

romm133.wmf romm134.wmf

romm135.wmf

romm136.wmf

Отсюда следует равномерная сходимость romm137.wmf к Y(x) на [x0, xfin]. С данными оценками формулируется аналог теоремы 1 для системы (22). Доказывается [10], что romm138.wmf равномерно на [x0, xfin] приближает производную Y′(x) при k →∞ с оценками:

romm139.wmf

romm140.wmf

и, при h < 1,

romm141.wmf

romm142.wmf

Итерационное уточнение

На [αr, βr] выполняются вариации, на основе которых выбираются и фиксируются степень полиномов n и количество подынтервалов 2k, затем на каждом подынтервале [xri, x ri+1], i = 0, 1, ..., 2k – 1, выполняется итерационное уточнение:

romm143.wmf romm144.wmf romm145.wmf

romm146.wmf romm147.wmf romm148.wmf

romm149.wmf romm150.wmf

где q – конечный номер итераций на предыдущем подынтервале [xri–1, xri]. Далее,

romm151.wmf

romm152.wmf

romm153.wmf

romm154.wmf

Если

romm155.wmf

то romm156.wmf

romm157.wmf romm158.wmf (26)

если же

romm159.wmf

l = 0, 1, ..., romm160.wmf

то [10]

romm161.wmf (27)

Согласно (26), (27) итерационное уточнение уменьшает погрешность кусочно-интерполяционных приближений решения задачи (22) аналогично одномерному случаю.

Варьируемое кусочно-интерполяционное приближение решения двухточечной задачи Коши для системы ОДУ

Двухточечная задача Коши для системы (22) на отрезке [x0, xfin] ставится в виде

Y′ = F(x, Y); Y(x0) = Y0; Y(xfin) = Yfin,

где Y, F(x, Y) – те же, что в (22), Y0 = (y01, y02, ..., y0N); Yfin = (yfin1, yfin2, ..., yfinN). Для ее решения отрезок [x0, xfin] делится на две части, примерно, возле середины:

romm162.wmf.

На отрезке [x0, xpr] задача решается в точности описанным выше способом (движение вдоль этого отрезка обозначается romm163.wmf). На отрезке [xpr, xfin] решение задачи выполняется в полной аналогии с предыдущим, однако за вектор начальных значений принимается (Y(xfin) = Yfin, и от точки xfin движение решения идет в обратном направлении – к xpr (движение в обратном направлении обозначается romm164.wmf). Обратное направление реализуется отрицательным значением шага –h во всех конечных разностях, построенных от xfin как от начальной точки, в каждом интерполяционном полиноме Ньютона на каждом подынтервале. Итерационное уточнение также строится в обратном направлении и имеет вид

romm165.wmf

romm166.wmf romm167.wmf

где romm168.wmf – правый конец подынтервала и romm169.wmf. Этих изменений достаточно для получения кусочно-интерполяционного приближения на romm170.wmf с теми же свойствами и с той же точностью, как в прямом направлении [10]. Приближения решения и его производной склеиваются на romm171.wmf и на romm172.wmf в точке xpr с помощью интерполяционного полинома Ньютона той степени, которая получена на правой границе xpr по направлению вдоль romm173.wmf. При этом часть узловых значений заимствуется слева от xpr, часть – справа от этой точки на romm174.wmf. Полученный полином подвергается итерационному уточнению. Метод реализует программа [10], результат работы которой ниже приводится для задачи

romm175.wmf romm176.wmf y1(1) = 2; y2(1) = 4, .... y1(512) = 262656; y2(512) = 263169:

1.00, 6.12, … , 26.60, 31.72, 36,84, … , 509.44, 5.12

0.00, 0.00, , … , 0.00, 1.11×10–16, 0.00, …, 0.00, 0.00

В первой строке – значения независимой переменной, во второй – соответственная норма абсолютной погрешности приближенного решения. Точность приближения сравнима с точностью изложенного метода для соответствующих одноточечных задач Коши. Программы, эксперименты, доказательства приведенных утверждений приведены в [10].

Заключение

Компьютерный метод приближенного решения одноточечной и двухточечной задач Коши для ОДУ дает аналитическое приближение решения и его производной на основе кусочной интерполяции путем соединения узлов на смежных границах подынтервалов и задания в них равных узловых значений. С ростом числа подынтервалов метод равномерно сходится к решению, одновременно равномерно приближается производная. Скорость сходимости обратно пропорциональна числу подынтервалов. На основе преобразования интерполирующих решение и производную полиномов к алгебраическому виду с числовыми коэффициентами выполняется итерационное уточнение по типу последовательных приближений Пикара. Достигается сравнительно высокая точность приближения решения и производной для жестких и нежестких задач [10].