Интерполяционные методы решения дифференциальных уравнений

Задачи с начальными условиями для систем обыкновенных дифференциальных уравнений

Рассмотрим задачу Коши для системы обыкновенных дифференциальных уравнений $$ \begin \tag <1>\frac &= f_i (t, u_1, u_2, \ldots, u_n), \quad t > 0\\ \tag <2>u_i(0) &= u_i^0, \quad i = 1, 2, \ldots, m. \end $$

Используя векторные обозначения, задачу (1), (2) можно записать как задачу Коши $$ \begin \tag <3>\frac> &= \pmb(t, \pmb), \quad t > 0, \\ \tag <4>\pmb(0) &= \pmb_0 \end $$ В задаче Коши необходимо по известному решению в точке $ t = 0 $ необходимо найти из уравнения (3) решение при других $ t $.

Численные методы решения задачи Коши

Существует большое количество методов численного решения задачи (3), (4). Вначале рассмотрим простейший явный метод Эйлера и его программную реализацию. Затем будут представлены методы Рунге—Кутта и многошаговые методы.

При построении численных алгоритмов будем считать, что решение этой дифференциальной задачи существует, оно единственно и обладает необходимыми свойствами гладкости.

Идея численных методов решения задачи (3), (4) состоит из четырех частей:

1. Вводится расчетная сетка по переменной $ t $ (время) из $ N_t + 1 $ точки $ t_0 $, $ t_1 $, $ \ldots $, $ t_ $. Нужно найти значения неизвестной функции $ \pmb $ в узлах сетки $ t_n $. Обозначим через $ \pmb^n $ приближенное значение $ \pmb(t_n) $.

2. Предполагаем, что дифференциальное уравнение выполнено в узлах сетки.

3. Аппроксимируем производные конечными разностями.

4. Формулируем алгоритм, который вычисляет новые значения $ \pmb^ $ на основе предыдущих вычисленных значений $ \pmb^k $, $ k 0 $ при $ \tau\to 0 $.

Явный метод Эйлера

Проиллюстрируем указанные шаги. Для начала введем расчетную сетку. Очень часто сетка является равномерной, т.е. имеет одинаковое расстояние между узлами $ t_n $ и $ t_ $: $$ \omega_\tau = \< t_n = n \tau, n = 0, 1, \ldots, N_t \>. $$

Затем, предполагаем, что уравнение выполнено в узлах сетки, т.е.: $$ \pmb^\prime (t_n) = \pmb(t_n, u(t_n)), \quad t_n \in \omega_\tau. $$

Заменяем производные конечными разностями. С этой целью, нам нужно знать конкретные формулы, как производные могут быть аппроксимированы конечными разностями. Простейший подход заключается в использовании определения производной: $$ \pmb^\prime(t) = \lim_ <\tau \to 0>\frac<\pmb(t+\tau) — \pmb(t)><\tau>. $$

В произвольном узле сетки $ t_n $ это определение можно переписать в виде: $$ \begin \pmb^\prime(t_n) = \lim_ <\tau \to 0>\frac<\pmb(t_n+\tau) — \pmb(t_n)><\tau>. \end $$ Вместо того, чтобы устремлять шаг сетки к нулю, мы можем использовать малый шаг $ \tau $, который даст численное приближение $ u^\prime(t_n) $: $$ \begin \pmb^\prime(t_n) \approx \frac<\pmb^ — \pmb^><\tau>. \end $$ Такая аппроксимация известна как разностная производная вперед и имеет первый порядок по $ \tau $, т.е. $ O(\tau) $. Теперь можно использовать аппроксимацию производной. Таким образом получим явный метод Эйлера: $$ \begin \tag <5>\frac<\pmb^ — \pmb^n> <\tau>= \pmb(t_n, \pmb^). \end $$

Четвертый шаг заключается в получении численного алгоритма. Из (5) следует, что мы должны знать значение $ y^n $ для того, чтобы решить уравнение (5) относительно $ y^ $ и получить формулу для нахождения приближенного значения искомой функции на следующем временном слое $ t_ $: $$ \begin \tag <6>\pmb^ = \pmb^n + \tau \pmb(t_n, \pmb^) \end $$

При условии, что у нас известно начальное значение $ \pmb^0 = \pmb_0 $, мы можем использовать (6) для нахождения решений на последующих временных слоях.

Программная реализация явного метода Эйлера

Выражение (6) может быть как скалярным так и векторным уравнением. И в скалярном и в векторном случае на языке Python его можно реализовать следующим образом

При решении системы (векторный случай), u[n] — одномерный массив numpy длины $ m+1 $ ($ m $ — размерность задачи), а функция F должна возвращать numpy -массив размерности $ m+1 $, t[n] — значение в момент времени $ t_n $.

Таким образом численное решение на отрезке $ [0, T] $ должно быть представлено двумерным массивом, инициализируемым нулями u = np.zeros((N_t+1, m+1)) . Первый индекс соответствует временному слою, а второй компоненте вектора решения на соответствующем временном слое. Использование только одного индекса, u[n] или, что то же самое, u[n, :] , соответствует всем компонентам вектора решения.

Функция euler решения системы уравнений реализована в файле euler.py:

Строка F_ = lambda . требует пояснений. Для пользователя, решающего систему ОДУ, удобно задавать функцию правой части в виде списка компонент. Можно, конечно, требовать чтобы пользователь возвращал из функции массив numpy , но очень легко осуществлять преобразование в самой функции решателе. Чтобы быть уверенным, что результат F будет нужным массивом, который можно использовать в векторных вычислениях, мы вводим новую функцию F_ , которая вызывает пользовательскую функцию F «прогоняет» результат через функцию assaray модуля numpy .

Неявный метод Эйлера

При построении неявного метода Эйлера значение функции $ F $ берется на новом временном слое, т.е. для решении задачи (5) используется следующий метод: $$ \begin \tag <7>\frac<\pmb^ — \pmb^n> <\tau>= \pmb(t_, \pmb^). \end $$

Таким образом для нахождения приближенного значения искомой функции на новом временном слое $ t_ $ нужно решить нелинейное уравнение относительно $ \pmb^ $: $$ \begin \tag <8>\pmb^ — \tau \pmb(t_, \pmb^) — y^n = 0. \end $$

Для решения уравнения (8) можно использовать, например, метод Ньютона.

Программная реализация неявного метода Эйлера

Функция backward_euler решения системы уравнений реализована в файле euler.py:

Отметим, что для нахождения значения u[n+1] используется функция fsolve модуля optimize библиотеки scipy . В качестве начального приближения для решения нелинейного уравнения используется значение искомой функции с предыдущего слоя u[n] .

Методы Рунге—Кутта

Одношаговый метод Рунге—Кутта в общем виде записывается следующим образом: $$ \begin \tag <9>\frac<\pmb^ — \pmb^n> <\tau>= \sum_^s b_i \pmb_i, \end $$ где $$ \begin \tag <10>\pmb_i = \pmb\left( t_n + c_i\tau, \pmb^n + \tau \sum_^s a_\pmb_j \right), \quad i = 1, 2, \ldots, s. \end $$ Формула (9) основана на $ s $ вычислениях функции $ \pmb $ и называется $ s $-стадийной. Если $ a_ = 0 $ при $ j \geq i $ имеем явный метод Рунге—Кутта. Если $ a_ = 0 $ при $ j > i $ и $ a_ \ne 0 $, то $ \pmb_i $ определяется неявно из уравнения $$ \begin \tag <11>\pmb_i = \pmb\left( t_n + c_i\tau, \pmb^n + \tau \sum_^ a_\pmb_j + \tau a_ \pmb_i \right), \quad i = 1, 2, \ldots, s. \end $$ О таком методе Рунге—Кутта говорят как о диагонально-неявном.

Одним из наиболее распространенных является явный метод Рунге-Кутта четвертого порядка: $$ \begin \tag <12>\pmb_1 & = \pmb(t_n, \pmb^n), &\quad \pmb_2 &= \pmb\left( t_n + \frac<\tau><2>, \pmb^n + \tau \frac<\pmb_1> <2>\right),\\ \pmb_3 &= \pmb\left( t_n + \frac<\tau><2>, \pmb^n + \tau \frac<\pmb_2> <2>\right), &\quad \pmb_4 &= \pmb\left( t_n + \tau, \pmb^n + \tau \pmb_3 \right),\\ \frac<\pmb^ -\pmb^n> <\tau>&= \frac<1> <6>(\pmb_1 + 2\pmb_2 + 2\pmb_3 + \pmb_4) & & \end $$

Многошаговые методы

В методах Рунге—Кутта в вычислениях участвуют значения приближенного решения только в двух соседних узлах $ \pmb^n $ и $ \pmb^ $ — один шаг по переменной $ t $. Линейный $ m $-шаговый разностный метод записывается в виде $$ \begin \tag <13>\frac<1> <\tau>\sum_^m a_i \pmb^ = \sum_^ b_i \pmb(t_, \pmb^), \quad n = m-1, m, \ldots \end $$ Вариант численного метода определяется заданием коэффициентов $ a_i $, $ b_i $, $ i = 0, 1, \ldots, m $, причем $ a_0 \ne 0 $. Для начала расчетов по рекуррентной формуле (13) необходимо задать $ m $ начальных значений $ \pmb^0 $, $ \pmb^1 $, $ \dots $, $ \pmb^ $ (например, можно использовать для их вычисления метод Эйлера).

Различные варианты многошаговых методов (методы Адамса) решения задачи с начальными условиями для систем обыкновенных дифференциальных уравнений могут быть получены на основе использования квадратурных формул для правой части равенства $$ \begin \tag <14>\pmb(t_) — \pmb(t_n) = \int_^> \pmb(t, \pmb) dt \end $$

Для получения неявного многошагового метода используем для подынтегральной функции интерполяционную формулу по значениям функции $ \pmb^ = \pmb(t_, \pmb^) $, $ \pmb^n $, $ \dots $, $ \pmb^ $, т.е. $$ \begin \tag <15>\frac<\pmb^ — \pmb^n> <\tau>= \sum_^ b_i \pmb(t_, \pmb^) \end $$

Для интерполяционного метода Адамса (15) наивысший порядок аппроксимации равен $ m+1 $.

Для построения явных многошаговых методов можно использовать процедуру экстраполяции подынтегральной функции в правой части (14). В этом случае приближение осуществляется по значениям $ \pmb^n $, $ \pmb^ $, $ \dots $, $ \pmb^ $ и поэтому $$ \begin \tag <16>\frac<\pmb^ — \pmb^n> <\tau>= \sum_^ b_i \pmb(t_, \pmb^) \end $$

Для экстраполяционного метода Адамса (16) погрешность аппроксимации имеет $ m $-ый порядок.

На основе методов Адамса строятся и схемы предиктор–корректор. На этапе предиктор используется явный метод Адамса, на этапе корректора — аналог неявного метода Адамса. Например, при использовании методов третьего порядка аппроксимации в соответствии с (18) для предсказания решения положим $$ \frac<\pmb^ — \pmb^n> <\tau>= \frac<1> <12>(23 \pmb^ -16\pmb^ + 5\pmb^). $$ Для уточнеия решения (см. (17)) используется схема $$ \frac<\pmb^ — \pmb^n> <\tau>= \frac<1> <24>(9\pmb^ + 19\pmb^ — 5\pmb^ + \pmb^). $$ Аналогично строятся и другие классы многошаговых методов.

Жесткие системы ОДУ

При численном решении задачи Коши для систем обыкновенных дифференциальных уравнений (3), (4) могут возникнуть дополнительные трудности, порожденные жесткостью системы. Локальные особенности поведения решения в точке $ u = w $ передаются линейной системой $$ \begin \frac

= \sum_^ \frac<\partial f_i> <\partial u_j>(t, w) v + \bar(t), \quad t > 0. \end $$
Пусть $ \lambda_i(t) $, $ i = 1, 2, \ldots, m $ — собственные числа матрицы $$ \begin A(t) = \< a_(t) \>, \quad a_(t) = \frac<\partial f_i><\partial u_j>(t, w). \end $$ Система уравнений (3) является жесткой, если число $$ \begin S(t) = \frac <\max_<1 \leq i \leq m>|Re \lambda_i(t)|> <\min_<1 \leq i \leq m>|Re \lambda_i(t)|> \end $$ велико. Это означает, что в решении присутствуют составляющие с сильно различающимися масштабами изменения по переменной $ t $.
Для численное решения жестких задач используются вычислительные алгоритмы, которые имеют повышенный запас устойчивости. Необходимо ориентироваться на использование $ A $-устойчивых или $ A(\alpha) $-устойчивых методов.
Метод называется $ A $-устойчивым, если при решении задачи Коши для системы (3) область его устойчивости содержит угол $$ \begin |\arg(-\mu)| —>
Методы интегрально-дифференциальной интерполяции
Постановка задачи
Использование различных сочетаний функциональных, дифференциальных и интегральных условий согласования (4.2), (4.3), (4.5) позволяет конструировать различные типы интерполяционных или интерполяционно-сглаживающих многочленов. В данном разделе рассматриваются задачи построения двух интегрально-дифференциальных интерполяционных параболических многочленов [math]S_2^<(I)>(x)[/math] и [math]S_2^<(II)>(x)[/math] .
Первый определяется по сеточной функции, заданной не только своими значениями [math]f_i
(i=\overline<0,n-1>)[/math] , но и значениями интегралов [math]I_^
(i=\overline<0,n-1>)[/math] . Так как для нахождения многочлена [math]S_2^<(I)>(x)[/math] используются совместно функциональное условие (4.3) и интегральное условие (4.5), то многочлен называется интегрально-функциональным. Поскольку функциональное условие (4.3) является частным случаем дифференциального условия (4.2), соответствующий метод относится к интегрально-дифференциальным.
Второй многочлен определяется по сеточной функции, заданной не только значениями производных [math]f_^
(i=\overline<0,n>)[/math] , но и значениями интегралов [math]I_^
(i=\overline<0,n-1>)[/math] . Так как для нахождения многочлена [math]S_2^<(II)>(x)[/math] совместно используются дифференциальное условие (4.2) при [math]p=1[/math] и интегральное условие (4.5), то многочлен (и соответствующий метод) называется интегрально-дифференциальным.
Интерполяционный параболический интегрально-функциональный многочлен
Пусть некоторая сеточная функция, определенная в общем случае на неравномерной сетке [math]\Omega_n= \[/math] задана своими значениями [math]f_0,f_1,\ldots, f_n[/math] и интегралами [math]I_0^1,I_1^2,\ldots,I_^n[/math] [math]\textstyle<\Big(I_^= \int\limits_^> f(x)dx\Big)>[/math] , которые могут быть предварительно рассчитаны по одной из квадратурных формул. Будем считать, что точность заданных значений [math]f_i[/math] и вычисленных интегралов [math]I_^[/math] для всех узлов [math]x_i
(i=\overline<0,n>)[/math] и отрезков [math][x_i,x_]
(i=\overline<0,n-1>)[/math] не ниже [math]O(h_^3)[/math] и [math]O(h_^4)[/math] соответственно. Это требование следует из утверждения В.2, устанавливающего принцип соответствия порядков аппроксимации (интерполяции).
Применим кусочный способ интерполяции. Для некоторого отрезка [math][x_i,x_]
(0 \leqslant i \leqslant n-1)[/math] на двухточечном шаблоне требуется определить многочлен
где [math]a_<0,i>,\, a_<1,i>,\, a_<2,i>[/math] — неизвестные коэффициенты, удовлетворяющие функциональным (4.3) и интегральным условиям согласования (4.5):
Для определения тех коэффициентов параболического интегрально-функционального многочлена можно получить фи линейных уравнения, являющиеся следствием трех условий согласования:
Вводя фазу интерполяции [math]u=\frac>[/math] , последнее условие после интегрирования можно записать в виде
Итак, для оставшихся двух коэффициентов [math]a_<1,i>,\, a_<2,i>[/math] имеется система линейных алгебраических уравнений:
\Delta f_i= f_-f_i[/math] . Очевидно, определитель этой системы [math]\begin3h_^2& 2h_^3\\ h_& h_^2 \end\ne0[/math] , что свидетельствует о существовании и единственности интегрально-функционального многочлена [math]S_<2,i>^<(I)>(x)[/math] на отрезке [math][x_i,x_][/math] . В силу произвольности отрезка [math][x_i,x_][/math] данный многочлен [math]S_<2,i>^<(I)>(x)[/math] существует и является единственным на каждом частичном отрезке [math][x_i,x_],
i=\overline<0,n-1>[/math] . Поэтому многочлен [math]S_<2,i>^<(I)>(x)[/math] может быть положен в основу кусочной многозвенной интерполяции, в том числе и сплайн-интерполяции.
Определяя [math]a_<1,i>,\, a_<2,i>[/math] из приведенной системы, записываем параболический интегрально-функциональный многочлен [math]S_<2,i>^<(I)>(x)[/math] в полиномиальной форме:
Если в (4.37) подставить вместо [math](x-x_i)[/math] произведение [math]uh_[/math] и выделить слагаемые при [math]I_^,\, f_i,\, f_[/math] , то получим параболический интегрально-функциональный многочлен в лагранжевой форме:
P_<2,i+1>(u)= u(3u-2)[/math] — коэффициенты Лагранжа, удовлетворяющие функциональным и интегральному условиям:
Так же как и классические (функциональные) многочлены Лагранжа, ин-тефально-функциональный многочлен [math]S_<2,i>^<(I)>(x)[/math] линейно зависит от параметров [math]I_^,\, f_i,\, f_[/math] . При использовании кусочной интерполяции для определения значения функции [math]f(x_<\ast>)[/math] в точке [math]x_<\ast>[/math] с помощью многочлена [math]S_<2,i>^<(I)>(x)[/math] выполняются действия, аналогичные применяемым при функциональной интерполяции.
Методика решения задачи интегрально-функциональной интерполяции
1. Выбрать «окно» интерполяции [math][x_i,x_][/math] так, чтобы [math]x_i \leqslant x_<\ast>\leqslant x_[/math] .
2. Определить фазу интерполяции [math]u=\frac-x_i>>[/math] и коэффициенты [math]P_<2,k>(u)
3. По формуле (4.38) вычислить значение [math]f(x_<\ast>)\approx S_<2,i>^<(I)>(u)[/math] .
1. Правильность вычисления коэффициентов Лагранжа проверяется по условию [math]P_<2,I>+ P_<2,i>+ P_<2,i+1>=1[/math] получающемуся путем подстановки в (4.38) функции
2. При аппроксимации значения [math]f(x_<\ast>)[/math] значением многочлена [math]S_<2,i>^<(I)>(x_<\ast>)[/math] возникает погрешность [math]R_<2>^<(I)>(x_<\ast>)\colon[/math] [math]f(x_)= S_<2,i>^<(I)>(x_<\ast>)+ R_<2>^<(I)>(x_<\ast>)[/math] . При условии, что [math]f(x)\in C_3[a,b][/math] , можно получить оценку погрешности на отрезке [math][x_i,x_]\colon[/math]
Если класс гладкости функции [math]f(x)[/math] понижен и [math]f(x)\in C_2[a,b][/math] , то порядок аппроксимации интегрально-функциональным многочленом [math]S_<2,i>^<(I)>(x)[/math] также понижается на единицу:
3. Сравним условия применения функционального (классического) многочлена Лагранжа [math]L_2(x)[/math] , определенного на трехточечном шаблоне [math](x_i,x_,x_)[/math] , и интегрально-функционального многочлена [math]S_<2,i>^<(I)>(x)\colon[/math]
а) аппроксимируемая функция для применения многочлена Лагранжа задается только своими значениями [math]f_i
(i=\overline<0,n>)[/math] , т.е. функциональными параметрами, а для использования [math]S_<2>^<(I)>(x)[/math] функция определяется еще и интегралами [math]I_^
(i=\overline<0,n-1>)[/math] . Если в постановке задачи интегралы неизвестны, то, как отмечено выше, они вычисляются предварительно по значениям [math]f_i
(i= \overline<0,n>)[/math] по тем или иным квадратурным формулам. Интегрально-функциональная интерполяция предоставляет дополнительные возможности для более полного учета локальных свойств интерполируемых функций, которые в отдельных точках и областях могут иметь особенности: разрывы производных, локальные экстремумы и др. При этом локальные свойства учитываются путем соответствующего выбора трехточечного шаблона, на котором аппроксимируется интеграл [math]I_^[/math] ;
б) многочлен [math]L_2(x)[/math] записывается на трехточечном шаблоне [math](x_i,x_, x_)[/math] , а интегрально-функциональный многочлен [math]S_<2,i>^<(I)>(x)[/math] — на двухточечном [math](x_i,x_)[/math] ;
в) функциональный многочлен [math]L_2(x)[/math] при его записи на одном шаблоне определяется только функциональными параметрами [math]f_i,f_,f_[/math] , a интегрально-функциональный многочлен [math]S_<2,i>^<(I)>(x)[/math] — двумя функциональными и интегральным [math]I_^[/math] , поэтому при его построении можно учесть интегральные свойства аппроксимируемой функции;
г) по порядку аппроксимации оба многочлена обеспечивают одинаковый третий порядок аппроксимации, который при [math]f(x)\in C_3[a,b][/math] является максимальным, однако константа в мажоранте для [math]S_<2,i>^<(I)>(x)[/math] в восемь раз меньше соответствующей константы в мажоранте для [math]L_2(x)\colon[/math]
Пример 4.7. Пусть сеточная функция на отрезке [math][-1;4][/math] задана функциональными [math]f_i
(i=\overline<0,4>)[/math] и интегральными параметрами [math]I_^,
i=\overline<0,3>[/math] , приведенными в табл. 4.12. Данные параметры получены точно путем «обработки» сеточного представления функции [math]f(x)=x^3[/math] . интегралы вычислены точно по формуле Ньютона-Лейбница [math]I_^= F_-F_i[/math] , где [math]F_i[/math] — первообразная. Значения интегралов в табл. 4.12 помещены посередине между значениями функции. Требуется найти значение [math]f(x_<\ast>)[/math] при [math]x_<\ast>=2[/math] на основе интегрально-функционального многочлена [math]S_<2,i>^<(I)>(x)[/math] .
1. По значению [math]x_<\ast>=2[/math] выбираем «окно» интерполяции [math][x_2,x_3]= [1;3][/math] , такое, что [math]x_2 .
2. Вычислим [math]u=\frac-x_2>= \frac<2-1><2>= \frac<1><2>[/math] и коэффициенты Лагранжа
Так как [math]\textstyle<\sum\limits_P_<2,k>=1>[/math] коэффициенты вычислены правильно.
3. Найдем значение /(2) по формуле (4.38) при [math]i=2\colon[/math]
Получилось значение [math]S_<2,2>^<(I)>(2)=8[/math] , совпадающее с точным значением исходной формульной функции.
Пример 4.8. Дана сеточная функция, характеризующаяся функциональными и интегральными параметрами, заданными в табл. 4.12. Найти полиномиальный вид интегрально-функционального многочлена на отрезке [math][x_2,x_3][/math] , внутри которого находится точка [math]x_<\ast>=2[/math] , и вычислить [math]S_<2>^<(I)>(2)[/math] . Результаты сравнить с решением, полученным с помощью многочлена Лагранжа в примере 4.3.
Для решения задачи в соответствии с исходными данными используем интегрально-функциональный многочлен [math]S_<2>^<(I)>(x)[/math] , записанный в форме (4.37).
1. В поставленной задаче «окно» интерполяции, как в примере 4.7, задано отрезком [math][x_2,x_3]= [1;3][/math] , поэтому [math]i=2[/math] .
2. Вычислим значения [math]\nabla I_2^3,\, \Delta f_2,\, h_3[/math] и коэффициенты многочлена [math]S_<2,2>^<(I)>(x)\colon[/math]
Полученные коэффициенты подставим в формулу (4.37):
Из сопоставления полученных результатов с результатами решения примера 4.3 следует, что значение [math]S_<2,2>^<(I)>(2)[/math] , полученное интегрально-функциональным методом, лучше соответствует точному значению, чем [math]L_2(2)=6[/math] . Это обусловлено тем, что кроме значений функции [math]f_i=f(x_i)[/math] в узлах задан еще и интеграл, значение которого является точным. В данном случае значение [math]S_<2,2>^<(I)>(2)[/math] совпало с точным значением, так как [math]x_<\ast>=2[/math] находится в середине отрезка [math][1;3][/math] .
Интерполяционный параболический интегрально-дифференциальный многочлен
Пусть некоторая сеточная функция на сетке [math]\Omega_n[/math] задана значениями производных [math]f’_i=f'(x_i)
(i=\overline<0,n>)[/math] и интегралами [math]I_^
(i=\overline<0,n>)[/math] , которые либо известны, либо могут быть предварительно рассчитаны по одной из квадратурных формул. Так же, как и для многочлена [math]S_<2,i>^<(I)>(x)[/math] рассмотренного ранее, будем считать, что [math]f’_i[/math] и [math]I_^[/math] заданы или вычислены с точностью не ниже [math]O(h_^2)[/math] и [math]O(h_^4)[/math] соответственно.
Применим кусочный способ интерполяции. Для некоторого отрезка [math][x_i,x_]
(0 \leqslant i \leqslant n-1)[/math] на двухточечном шаблоне требуется определить многочлен
где [math]a_<0,i>,\, a_<1,i>,\, a_<2,i>[/math] — неизвестные коэффициенты, удовлетворяющие дифференциальным ((4.2) при [math]p=1[/math] ) и интегральным (4.5) условиям согласования:
Для нахождения трех коэффициентов многочлена можно получить три линейных уравнения, являющиеся следствием трех условий согласования, так же, как это сделано ранее. В результате имеем выражение для параболического интегрально-дифференциального многочлена [math]S_<2,i>^<(II)>(x)[/math] в полиномиальной форме (где [math]\Delta f’_i= f’_-f’_i[/math] ):
Если подставить вместо [math](x-x_i)[/math] произведение [math]u\,h_[/math] ( [math][/math] — фаза интерполяции) и выделить слагаемые при [math]I_^,\, f’_i,\, f_[/math] , то получим параболический интегрально-дифференциальный многочлен в лагранжевой форме (где [math]u=\frac>,
0 \leqslant u \leqslant 1[/math] ):
1. При аппроксимации значения [math]f(x_<\ast>)[/math] значением многочлена [math]S_<2,i>^<(II)>(x_<\ast>)[/math] возникает погрешность [math]R_<2>^<(II)>(x_<\ast>)\colon[/math] [math]f(x_<\ast>)= S_<2,i>^<(II)>(x_<\ast>)+ R_<2>^<(II)>(x_<\ast>)[/math] . При условии, что [math]f(x)\in C_3[a,b][/math] , можно получить оценку погрешности на отрезке [math][x_,x_]\colon[/math]
2. Алгоритм аппроксимации сеточных функций на основе (4.41), (4,42) аналогичен изложенному ранее для аппроксимации многочленом [math]S_<2,i>^<(I)>(x_<\ast>).[/math]
6.6. Интерполяционные методы Адамса
При S = 1 формула (6.16) примет вид
(6.22)
Если Q = 2, получим следующее вычислительное правило:
(6.23)
Обычно на практике используют экстраполяционную формулу (6.18), а затем корректируют полученное значение по формуле (6.23). И если результат уточненного значения не превышает допустимую погрешность расчета, то шаг H считается допустимым .
Для расчетов на компьютере формулы (6.18) и (6.23) в конечно-разностном виде неудобны. С учетом (6.21) их можно представить в виде
(6.24)
Приведенные формулы имеют достаточно большую точность. Они дают погрешность порядка
О( H4 ), но сами формулы оценки погрешности достаточно сложны. Приближенно погрешность можно оценить по правилу Рунге.
Пример 6.2. Решить дифференциальное уравнение на отрезке [0, 1] c начальным условием Y(X=0) = 1. Найти решение методом Адамса (с коррекцией) в точке X4, решение в трех первых точках найти методом Рунге- Кутта, приняв шаг .
Решение. Значения функции в четырех первых точках возьмем из табл. 6.1 (см. пример в предыдущем разделе). Теперь стало понятно, зачем мы сохраняли значения первой производной в этих точках (см. формулы (6.24)).
Для того чтобы скорректировать полученный результат, необходимо вычислить значение производной в этой точке:
Теперь уточним значение по интерполяционной формуле (а можно этого и не делать, тогда погрешность метода будет больше):
Так как в качестве нового значения функции принято скорректированное, то Обязательно Следует пересчитать значение производной. В нашем случае модуль разности экстраполяционной и интерполяционной формул меньше ε, Что позволяет продолжить вычисления с тем же шагом.
Вопросы для самопроверки
· Сформулируйте задачу Коши для обыкновенных дифференциальных уравнений первого порядка.
· Что является решением дифференциального уравнения: а) в высшей математике, б) в прикладной математике?
· Какие методы решения дифференциальных уравнений называются одношаговыми, многошаговыми? Приведите примеры.
· Сравните решения, полученные на первом, втором шаге методами Эйлера, Рунге-Кутта и разложением в ряд Тейлора (трудоемкость, погрешность…).
· Как оценить погрешность применяемого метода? Как ее уменьшить?
· Сравните одношаговые и многошаговые методы решения дифференциальных уравнений, указав достоинства и недостатки первых и вторых.
· Что такое экстраполяционные и интерполяционные методы (формулы) Адамса?
· Можно ли применять: а) только экстраполяционные методы Адамса,
б) только интерполяционные?
· Можно ли использовать: а) многошаговые методы без одношаговых;
б) одношаговые методы без многошаговых?
· При решении дифференциального уравнения методом Адамса на 27-м шаге необходимо сменить шаг. Как это сделать?

источники:
http://mathhelpplanet.com/static.php?p=metody-integralno-differentsialnoy-interpolyatsii
http://matica.org.ua/metodichki-i-knigi-po-matematike/chislennye-metody-iu-ia-katcman/6-6-interpoliatcionnye-metody-adamsa