Журнал вычислительной математики и математической физики, 2023, T. 63, № 9, стр. 1428-1437

Вычисление с минимальной погрешностью n-й производной по данным измерения функции

А. С. Демидов^{1, 2, 3, *}, А. С. Кочуров^1, 4, **

¹ МГУ
119234 Москва, ул. Колмогорова, 1, Россия

² МФТИ
123098 Москва, ул. Максимова, 4, Россия

³ РУДН
115419 Москва, ул. Орджоникидзе, 3, Россия

⁴ МЦФПМ-МГУ
119991 Москва, Россия

^* E-mail: demidov.alexandre@gmail.com
^** E-mail: kchrvas@yandex.ru

Поступила в редакцию 01.02.2023
После доработки 09.03.2023
Принята к публикации 29.05.2023

EDN: XMHVDC
DOI: 10.31857/S0044466923090077

Полный текст (PDF)

Аннотация

Предложено решение вопроса, возникающего во всех задачах, где по экспериментальным дискретным данным априори гладкой функции требуется приближенно вычислить ее производные. Вся проблема сводится к поиску “оптимального” шага разностной аппроксимации. Эту проблему исследовали многие математики. Оказалось, что для выбора “оптимального” шага аппроксимации производной $k$-го порядка надо знать как можно более точную оценку модуля производной порядка $k + 1$. Предложенный в статье алгоритм, дающий такую оценку, применен к задаче о концентрации тромбина, который определяет динамику свертываемости крови. Эта динамика представлена графиками и дает интересующий биофизиков ответ о концентрации тромбина. Библ. 8. Фиг. 7.

Ключевые слова: восстановление $n$-й производной, концентрация тромбина, оптимальный шаг аппроксимации, оценка модуля старших производных.

1. ВВЕДЕНИЕ

Задача, которой посвящена данная работа, инициирована проблемой, возникшей в ходе исследований, проводимых в лаборатории физической биохимии Национального медицинского исследовательского центра гематологии под руководством академика РАН Ф.И. Атауллаханова. Имеются в виду работы по нахождению концентрации C тромбина (фактора свертывания II) по данным интенсивности I, введенного в кровь фермента, активирующего ее свертываемость. Измерения интенсивности I размерности $[I]$, как функции времени $\theta \in [0,T]$ и глубины $\xi \in [0,X]$ проникновения в кровь, проводились 110 раз с шагом 31.2 sec за время $T = 109 \times 31.2[T] \approx $ ≈ 56.7 min, где $[T] = {\text{sec}}$ – размерность $T$. А именно, через каждые 31.2 sec фиксировалась фотография интенсивности проникновения фермента в кровь в 1485 точках, расположенных с шагом 4.3$[X]$ вплоть до глубины $X = 1484 \times 4.3[X] \approx 6.4$ mm, где $[X] = {\text{mkm}}$ – размерность $X$.

Следует отметить, что многие трудные и пока не решенные био-физиологические проблемы свертывания крови были освещены в недавнем обзоре [1] по современным методам компьютерного моделирования тромбоза. Но математическая проблема, с которой столкнулись в лаборатории Ф.И. Атауллаханова, в этом обзоре не обозначена. Как отмечается в статье [2], “в современной концепции свертывания крови является растущее согласие с тем, что для понимания регуляции свертывания крови in vivo и нарушений ее функции необходимо принимать во внимание ее пространственную неоднородность, диффузию и поток.” Это утверждение базируется, в частности, на одном из наиболее значимых, и часто цитируемом исследовании [3] о свертывании крови, где проведен детальный анализ соответствующей системы дифференциальных уравнений с частными производными. В простейшей модельной ситуации измеряемая интенсивность I введенного в кровь фермента и искомая концентрация C тромбина связаны уравнением д-иффузии

$\frac{{\partial I}}{{\partial \theta }} - D\frac{{{{\partial }^{2}}I}}{{\partial {{\xi }^{2}}}} = C(\theta ,\xi ),\quad \theta \in [0,T],\quad \xi \in [0,X],\quad D = 350\;\frac{{{{{{\text{(mkm)}}}}^{{\text{2}}}}}}{{{\text{sec}}}}.$

При переходе к безразмерным параметрам

$t = \frac{\theta }{T}{\kern 1pt} ,\quad x = \frac{\xi }{X},\quad a = D\frac{T}{{{{X}^{2}}}}$

получаем для нормированной интенсивности $Y = \frac{I}{T}$ размерности $\frac{{[I]}}{{[T]}}$ следующее уравнение:

(1)

${{Y}_{t}}(t,x) - a{{Y}_{{xx}}}(t,x) = f(t,x),\quad {\text{где}}\quad f(t,x) = C(tT,xX),\quad {\text{а}}\quad (t,x) \in [0,1] \times [0,1].$

Чтобы найти искомую концентрацию C тромбина, представленную функцией $f$, остается найти с приемлемой точностью производные ${{Y}_{t}}(t,x)$ и ${{Y}_{{xx}}}(t,x)$ функции $Y = \frac{I}{T}$. При этом

(2)

$Y({{t}_{i}},{{x}_{k}}) = \widetilde Y({{t}_{i}},{{x}_{k}}) + \delta ({{t}_{i}},{{x}_{k}}),\quad \left| {\delta ({{t}_{i}},{{x}_{k}})} \right| \leqslant \sigma \left| {\widetilde Y({{t}_{i}},{{x}_{k}})} \right|,$

где $\widetilde Y({{t}_{i}},{{x}_{k}})$ – экспериментальные значения в точках $({{t}_{i}},{{x}_{k}}) = ({{\theta }_{i}}{\text{/}}T,{{\xi }_{k}}{\text{/}}T) \in [0,1] \times [0,1]$ нормированной интенсивности введенного в кровь фермента, а $\sigma = 0.05$ – относительная точность измерения прибора.

Казалось бы, все просто: искомая концентрация тромбина задана простой формулой $f = {{\widetilde Y}_{t}} - a{{\widetilde Y}_{{xx}}}$, а дискретные данные функции $\widetilde Y$ известны.

Но как найти с приемлемой точностью производные ${{\widetilde Y}_{t}}({{t}_{i}},{{x}_{k}})$ и ${{\widetilde Y}_{{xx}}}({{t}_{i}},{{x}_{k}})$?

С этим вопросом к авторам статьи обратился Ф.И. Атауллаханов. И вопрос этот вовсе не праздный. Его трудность иллюстрируется на фиг. 1a и особенно на фиг. 1б “невразумительными” графиками разностных аппроксимаций производных ${{\left. {{{{\widetilde Y}}_{t}}(t,x)} \right|}_{{x = 0.202}}}$ и ${{\left. {{{{\widetilde Y}}_{{xx}}}(t,x)} \right|}_{{t = 0.367}}}$, полученных непосредственно по “вполне вразумительным” данным измерений $\widetilde Y(t,0.202)$ и $\widetilde Y(0.367,x)$, представленным на фиг. 2.

Фиг. 1.

(a) график ${{\widetilde Y}_{t}}(t,{{x}_{0}})$; (б) график ${{\widetilde Y}_{{xx}}}({{t}_{0}},x)$.

Фиг. 2.

(a) график $\widetilde Y(t,{{x}_{0}})$; (б) график $\widetilde Y({{t}_{0}},x)$.

Данные измерений $\widetilde Y$ гладкой функции $Y$ необходимо дискретны, а приближение дискретнозначной функции $\widetilde Y$ той или иной гладкой функцией $P$ не влечет близости ${{P}_{t}}(t,x)$ и ${{P}_{{xx}}}(t,x)$ к искомым производным ${{Y}_{t}}(t,x)$ и ${{Y}_{{xx}}}(t,x)$. Этот общеизвестный факт проиллюстрирован на фиг. 3, 4 . А именно, на фиг. 3a даны графики разности между $\widetilde Y(t,0.202)$ и аппроксимирующими полиномами степени n = 4 и n = 8 (с относительной точностью приближения порядка 0.037 и 0.014 соответственно), а на фиг. 3б графики разности между $\widetilde Y(0.367,x)$ и соответствующими аппроксимирующими полиномами степени n = 7 и n = 8 (с относительной точностью приближения порядка 0.082 и 0.029 соответственно).

Фиг. 3.

Фиг. 4.

(a) P¹ = P_t – для n = 4 и n = 8; (б) P² = P_xx – для n = 7 и n = 8 .

И несмотря на такую относительно приемлемую аппроксимацию функций $\widetilde Y(t,0.202)$ и $\widetilde Y(0.367,x)$ аппроксимирующими полиномами, различие производных этих полиномов, графики которых представлены на фиг. 4a и фиг. 4б, достигает $100\% $, что не позволяет сказать что-либо определенное об искомых производных ${{Y}_{t}}(t,x)$ и ${{Y}_{{xx}}}(t,x)$.

Обозначенная проблема построения производных функции по данным ее экспериментальных измерений относится к задаче о наилучшем приближении неограниченного линейного оператора линейными операторами на классе элементов банахова пространства. Эта задача разрабатывается со второй половины 20-го века многими видными математиками. Она и родственные ей экстремальные задачи подробно освещены в недавнем обзоре [4]. Но собственно обозначенная проблема, конкретизация которой представлена в разд. 1, там не исследовалась. Частично ответ на нее дан в разд. 2, где используется приведенный в разд. 3 алгоритм.

2. ВОПРОС АТАУЛЛАХАНОВА КАК ПРОБЛЕМА ВЫБОРА ШАГА АППРОКСИМАЦИИ

Экспериментальные измерения функции I проводились в лаборатории физической биохимии Гематологического Научного Центра РАМН посредством регистрации фото- или кино-камерой интенсивности люминесценции, введенного в кровь фермента, активирующего свертываемость крови. Эти, как и любые экспериментальные измерения, необходимо дискретны по времени и пространственной переменной. В данном случае по экспериментальным измерениям $\widetilde I({{\theta }_{i}}{\text{/}}T,{{\xi }_{k}}{\text{/}}X)$ были получены значения $\widetilde Y({{t}_{i}},{{x}_{k}})$ в дискретных точках ${{t}_{i}} = i\tau $ ($i = 0, \ldots ,109$) и ${{x}_{k}} = kh$ ($k = 0, \ldots ,1484$), где

$\tau = 31.2{\text{/}}3400.8 = 1{\text{/}}109\quad {\text{и}}\quad h = 4.3{\text{/}}6381.2 = 1{\text{/}}1484{\kern 1pt} ,$

соответствующие шагу измерения 31.2 sec по времени и шагу 4.3 mkm по глубине проникновения фермента I.

По значениям $\widetilde Y$ в точках $({{t}_{i}},{{x}_{k}}) \in [0,1] \times [0,1]$

$({{t}_{i}},{{x}_{k}}) = (i\tau ,kh),\quad \left( {{{t}_{i}} \pm j(\tau ,0),{\kern 1pt} {{x}_{k}} \pm l(0,h)} \right)({{t}_{i}} \pm j\tau ,{{x}_{k}} \pm lh),$

где $j$ и $l$ – некоторые натуральные числа, вычислялись приближения $\widetilde Y_{{j(\tau ,0)}}^{1}(t,x)$ и $\widetilde Y_{{l(0,h)}}^{2}(t,x)$ частных производных

${{Y}_{t}}(t,x) \approx \widetilde Y_{{j(\tau ,0)}}^{1}(t,x)\;\mathop = \limits^{{\text{def}}} \;\frac{{\widetilde Y(t + j\tau ,x) - \widetilde Y(t - j\tau ,x)}}{{2j\tau }},$

${{Y}_{{xx}}}(t,x) \approx \widetilde Y_{{l(0,h)}}^{2}(t,x)\;\mathop = \limits^{{\text{def}}} \;\frac{{\widetilde Y(t,x + lh) - 2\widetilde Y(t,x) + \widetilde Y(t,x - lh)}}{{{{{(lh)}}^{2}}}}.$

В том случае, когда $j = l = 1$, т.е. когда разностные аппроксимации производных вычисляются с шагом проведенных экспериментальных измерений, результаты вычислений $\widetilde Y_{{j(\tau ,0)}}^{1}(t,{{x}_{0}})$ и $\widetilde Y_{{l(0,h)}}^{2}({{t}_{0}},x)$ были выше представлены на фиг. 1 (для $k = 300$, $i = 40$, т.е. для ${{x}_{0}} = kh = \frac{{300}}{{1484}} \approx 0.202$, ${{t}_{0}} = i\tau = \frac{{40}}{{109}} \approx 0.367$). Бросается в глаза особенная по сравнению с графиком ${{\widetilde Y}^{1}} = \widetilde Y_{{(\tau ,0)}}^{1}(t,{{x}_{0}})$ “невразумительность” графика ${{\widetilde Y}^{2}} = \widetilde Y_{{(0,h)}}^{2}({{t}_{0}},x)$. Это связано не столько с тем, что это графики производным различных порядков (а именно, 1-го по $t$ и 2-го по $x$), сколько с существенным различием величины шагов $\tau = 1{\text{/}}109$ и $h = 1{\text{/}}1484$. Действительно, чем меньше шаг, тем больше отклонение от истинного значения производной за счет слагаемого в виде разностного отношения для погрешности измерения, в знаменателе которого стоит шаг аппроксимации. С другой стороны, с ростом шага аппроксимации растет остаточный член в аппроксимационной формуле Тейлора, что увеличивает погрешность аппроксимации рассматриваемой производной.

Проблема выбора шага аппроксимации обычно решается последовательным его увеличением (например, в два раза) и этот процесс завершается тогда, когда графическое изображение разностной аппроксимации визуально будет мало отличаться от графика гладкой функции. Но при этом есть немалый риск сгладить существенные особенности искомой производной. А там, где производная мала, на первый план выступают погрешности измерений в виде “шерстки”, с которыми можно мириться, поскольку это не затрагивает принципиальные особенности производной. Кроме того, конечно, желательно иметь более обоснованный выбор шага аппроксимации и избавиться от необходимости визуального наблюдения. Следует при этом сказать, что вопрос выбора оптимального шага аппроксимации ${{s}_{{{\text{opt}}}}}$ при реконструкции производной $n$-го порядка функции одной переменной $y:s \mapsto y(s)$, равной, аналогично (2), сумме измеренной и погрешности δ(s), был рассмотрен в работе [5] для $n = 1$ и $n = 2$, в [6] для $n = 3$ и $n = 4$, а для произвольного $n$ – в [7]. Так, С.Б. Стечкин доказал, что шаг аппроксимации для функции ${{s}_{{{\text{opt}}}}}$ в следующих оценках:

(3)

$\left| {\dot {y}(s) - \frac{{\tilde {y}(s + {{s}_{{{\text{opt}}}}}) - \tilde {y}(s - {{s}_{{{\text{opt}}}}})}}{{2{{s}_{{{\text{opt}}}}}}}} \right| \leqslant \frac{{\left\| \delta \right\|}}{{{{s}_{{{\text{opt}}}}}}} + \frac{{{{s}_{{{\text{opt}}}}}}}{2}\left\| {\ddot {y}} \right\| = \sqrt {2\left\| \delta \right\| \cdot \left\| {\ddot {y}} \right\|} ,\quad {\text{где}}\quad {{s}_{{{\text{opt}}}}} = \sqrt {\frac{{2\left\| \delta \right\|}}{{\left\| {\ddot {y}} \right\|}}} ,$

(4)

$\left| {\ddot {y}(s) - \frac{{\tilde {y}(s + {{s}_{{{\text{opt}}}}}) - 2\tilde {y}(s) + \tilde {y}(s - {{s}_{{{\text{opt}}}}})}}{{s_{{{\text{opt}}}}^{2}}}} \right| \leqslant \frac{{4\left\| \delta \right\|}}{{s_{{{\text{opt}}}}^{2}}} + \frac{{{{s}_{{{\text{opt}}}}}}}{3}\left\| {{{y}^{{(3)}}}} \right\|,\quad {\text{где}}\quad {{s}_{{{\text{opt}}}}} = {{\left( {\frac{{24\left\| \delta \right\|}}{{\left\| {{{y}^{{(3)}}}} \right\|}}} \right)}^{{1{\text{/}}3}}}$

для функции $y \in {{C}^{{n + 1}}}$ (соответственно, $n = 1,\;2$), будет оптимальным, в том смысле, что неравенства (3), (4), где $\left\| \delta \right\| = \mathop {\sup }\limits_s \left| {\delta (s)} \right| \approx \sigma \left\| {\tilde {y}} \right\|$, выполняются для всех функций $y \in {{C}^{{n + 1}}}$ и обращаются в равенства для некоторых из них.

Поскольку функция $\tilde {y}$ измеряется в дискретном наборе точек ${{s}_{i}}$, параметр ${{s}_{{{\text{opt}}}}}$ необходимо заменяется на некоторое ближайшее к нему. В том случае, когда ${{s}_{i}} = i\eta $, параметр ${{s}_{{{\text{opt}}}}}$ заменяется на

(5)

${{\Delta }_{{{\text{opt}}}}} = \mu \eta \in \left[ {(\mu - 1{\text{/}}2)\eta ,(\mu + 1{\text{/}}2)\eta } \right),$

где натуральное число $\mu $ соответствует оптимальному шагу $\mu \eta $ аппроксимации рассматриваемой производной.

В реальных задачах, кроме принадлежности функции $y:\mathbb{R} \ni s \mapsto y(s)$ к классу ${{C}^{n}}$, практически всегда имеется дополнительная информация, например, бывает известно, где функция сильно (мало) меняется, где она принимает большие (малые) значения. Таким образом, возникает чисто математическая проблема: улучшить оценки (3) и (4) при наличии той или иной дополнительной информации. Впрочем и формулы (3), (4) дают ориентир для выбора подходящего шага приближения производной ${{y}^{{(n)}}}$. Но для этого требуется как можно более точно оценить норму производной ${{y}^{{(n + 1)}}}$. Алгоритм, позволяющий получить такую оценку для большинства реальных задач, дан в разд. 3.

3. ОТВЕТ НА ВОПРОС АТАУЛЛАХАНОВА

В случае упомянутой во введении задачи о концентрации тромбина, приближенные значения $\widetilde Y_{t}^{'}(t,x)$ и $\widetilde Y_{x}^{{''}}(t,x)$ производных $Y_{t}^{'}(t,x)$ и $Y_{x}^{{''}}(t,x)$ в левой части уравнения (1) задавались, соответственно, в виде

$\widetilde Y_{t}^{'}\left( {t,x_{*}^{{}}} \right) = \dot {y}(t)\quad {\text{для}}\quad x_{*}^{{}} = kh,\quad h = 1{\text{/}}1484,\quad k = 0, \ldots ,1484,$

$\widetilde Y_{t}^{{''}}\left( {t_{*}^{{}},x} \right) = \ddot {y}(t)\quad {\text{для}}\quad t_{*}^{{}} = i\tau ,\quad \tau = 1{\text{/}}109,\quad i = 0, \ldots ,109{\kern 1pt} ,$

где производные $\dot {y}(s)$ и $\ddot {y}(s)$ функции $y:[0,1] \ni s \mapsto y(s)$ вычислялись по формулам (3) и (4) с заменой $s$ на ${{s}_{i}} = i\eta $ и ${{s}_{{{\text{opt}}}}}$ на ${{\Delta }_{{{\text{opt}}}}}\mathop = \limits^{(5)} \mu \eta $. При этом нормы $\left\| {\ddot {y}} \right\|$ и $\left\| {{{y}^{{(3)}}}} \right\|$ оценивались, следуя алгоритму разд. 4.

Графики $\widetilde Y_{t}^{'}\left( {t,x_{*}^{{}}} \right)$ и $\widetilde Y_{x}^{{''}}({{t}_{*}},x)$ изображены на фиг. 5. Сравнение этих графиков с представленными на фиг. 1 демонстрируют важность выбора оптимального шага аппроксимаций производных.

Фиг. 5.

На фиг. 6 представлены графики

${{\widetilde Y}_{t}}:[0,1] \times [0,1] \ni (t,x) \mapsto {{Y}_{t}}(t,x)\quad {\text{и}}\quad {{\widetilde Y}_{{xx}}}(t,x):[0,1] \times [0,1] \ni (t,x) \mapsto {{Y}_{{xx}}}(t,x).$

Фиг. 6.

А на фиг. 7 представлен график искомого приближения правой части $f(t,x) = C(tT,xX)$ формулы (1). При этом шаг для вычисления производной по каждой переменной определялся в два этапа. Сначала при каждом фиксированном значении $x$ находился шаг по $t$ в соответствии с представленным в п. 3 алгоритмом. А затем в качестве итогового шага по $t$ был взят шаг, равный среднему арифметическому шагов, вычисленных при каждом фиксированном значении $x$. Он оказался равным $10\tau $, в то время как минимальный (соответственно, максимальный) составлял $7\tau $ (соответственно, $18\tau )$. Шаг по переменной $x$ выбирался переменой местами $t$ и $x$. Этот средний шаг оказался равным $24h$, в то время как минимальный (соответственно, максимальный) составил $7h$ (соответственно, $35h)$.

Фиг. 7.

Отметим также, что на прямоугольнике $[{{\Delta }_{{{\text{op}}{{{\text{t}}}_{t}}}}},1 - {{\Delta }_{{{\text{op}}{{{\text{t}}}_{t}}}}}] \times [{{\Delta }_{{{\text{op}}{{{\text{t}}}_{x}}}}},1 - {{\Delta }_{{{\text{op}}{{{\text{t}}}_{x}}}}}]$ производные были подсчитаны по симметричным разностям, а вне этого прямоугольника по соответствующим левым или правым разностям.

4. ОЦЕНКА НОРМЫ ПРОИЗВОДНЫХ ФУНКЦИИ ОДНОЙ ПЕРЕМЕННОЙ ПО ЕЕ ДИСКРЕТНЫМ ЗНАЧЕНИЯМ

Обозначим через $N$ число точек ${{s}_{i}} = i\eta \in [0,1]$, $i = 0, \ldots ,N - 1$, в которых известны приближенные значения $\tilde {y}:s \mapsto y(s) - \delta (s)$ гладкой функции $y$. Для выбора оптимального шага аппроксимации производной ${{y}^{{(n)}}}$ порядка $n$ требуется, как было отмечено выше, знание весьма точной оценки нормы

$\left\| {{{y}^{{(n + 1)}}}} \right\| = \mathop {\max }\limits_{s \in [0,1]} \left| {{{y}^{{(n + 1)}}}(s)} \right|.$

С этой целью введем двумерный массив отрезков $[{{s}_{i}},{{s}_{j}}]$ по $i = 0, \ldots ,N - 1$ и $j = i + (n + 1)\ell \leqslant N - 1$, где $\ell $ – натуральное число, такое, что ${{s}_{j}} = \eta j \leqslant 1$. Заметим, что

$\left\| {{{y}^{{(n + 1)}}}} \right\| = \mathop {\max }\limits_{i,j} \mathop {\max }\limits_{s \in [{{s}_{i}},{{s}_{j}}]} \left| {{{y}^{{(n + 1)}}}(s)} \right|.$

Для функции $\tilde {y} = s \mapsto y(s) - \delta (s)$ рассмотрим так называемую (cм., например, [8] стр. 231) разделенную разность $(n + 1)$-го порядка в точке ${{s}_{i}}$ с шагом $\ell \eta $, относящуюся к отрезку $[{{s}_{i}},{{s}_{{i + (n + 1)\ell }}}]$. Она задается формулой

(6)

$\tilde {y}_{{i,j}}^{{(n + 1)}}\;\mathop = \limits^{{\text{def}}} \;\frac{{\sum\limits_{k = 0}^{n + 1} \,{{{( - 1)}}^{{n + 1 - k}}}C_{{n + 1}}^{k}\tilde {y}({{s}_{{i + k\ell }}})}}{{{{{(\ell \eta )}}^{{n + 1}}}}} = \frac{{\sum\limits_{k = 0}^{n + 1} \,{{{( - 1)}}^{{n + 1 - k}}}C_{{n + 1}}^{k}y({{s}_{{i + k\ell }}})}}{{{{{(\ell \eta )}}^{{n + 1}}}}} - \frac{{\sum\limits_{k = 0}^{n + 1} \,{{{( - 1)}}^{{n + 1 - k}}}C_{{n + 1}}^{k}\delta ({{s}_{{i + k\ell }}})}}{{{{{(\ell \eta )}}^{{n + 1}}}}}$

и удовлетворяет условию

(7)

${{y}^{{(n + 1)}}}(\zeta ) = \tilde {y}_{{i,j}}^{{(n + 1)}} + O(\left\| \delta \right\|{\text{/}}{{(\ell \eta )}^{{n + 1}}}){\kern 1pt} ,$

где $\zeta $ – некоторая точка на $[{{s}_{i}},{{s}_{j}}] = [{{s}_{i}},{{s}_{{i + (n + 1)\ell }}}]$, а $\left\| \delta \right\| = \mathop {\sup }\limits_s \left| {\delta (s)} \right|$ есть норма погрешности измерения. Учитывая (7), будем принимать за искомое значение нормы $\left\| {{{y}^{{(n + 1)}}}} \right\|$ оценку такой величины

$\left\| {{{{\tilde {y}}}^{{(n + 1)}}}} \right\|\;\mathop = \limits^{{\text{def}}} \;\mathop {\max }\limits_{i,j} \left| {\tilde {y}_{{i,j}}^{{(n + 1)}}} \right|{\kern 1pt} {\text{,}}\quad j = i + (n + 1)\ell {\kern 1pt} .$

Как следует из формулы (6), величины $\left| {y_{{i,j}}^{{(n + 1)}}} \right|$, а тем более, $\left| {\tilde {y}_{{i,j}}^{{(n + 1)}}} \right|$ очень сильно зависят от положения и длины отрезка $[{{s}_{i}},{{s}_{j}}]$. Это затрудняет поиск устойчивого максимума этих величин, который можно было бы принять за приемлемую оценку $\left\| {{{{\tilde {y}}}^{{(n + 1)}}}} \right\|$. Ниже излагается алгоритм из трех последовательных действий (шагов), позволяющий преодолеть эту трудность. На первых двух шагах элементы $\left| {y_{{i,j}}^{{(n + 1)}}} \right|$ выстраиваются в неубывающую последовательность, в “хвосте” которой концентрируются разностные отношения, вычисленные с большой погрешностью. Это построение осуществляется с помощью функций $\mathcal{Y}$ и $\mathcal{M}$, формальное определение которых иллюстрируется простым примером. “Хвост” последовательности отделяется на 3-м шаге алгоритма параметром $\hat {\varkappa }$, а за приемлемую оценку $\left\| {{{{\tilde {y}}}^{{(n + 1)}}}} \right\|$ принимается $\mathcal{M}(\hat {\varkappa })$.

Шаг 1. Из двумерного массива элементов $\left| {\tilde {y}_{{i,j}}^{{(n + 1)}}} \right|$ формируется набор

${{\mathcal{Y}}^{{(n + 1)}}} = \left\{ {\mathcal{Y}_{0}^{{(n + 1)}},\mathcal{Y}_{1}^{{(n + 1)}}, \ldots ,\mathcal{Y}_{{N - n - 2}}^{{(n + 1)}}} \right\}$

из семейств

$\mathcal{Y}_{i}^{{(n + 1)}} = \left\{ {\left| {\tilde {y}_{{i,i + 1(n + 1)}}^{{(n + 1)}}} \right|,\;\left| {\tilde {y}_{{i,i + 2(n + 1)}}^{{(n + 1)}}} \right|,\;\left| {\tilde {y}_{{i,i + 3(n + 1)}}^{{(n + 1)}}} \right|,\; \ldots ,\;\left| {\tilde {y}_{{i,i + {{k}_{i}}(n + 1)}}^{{(n + 1)}}} \right|} \right\},\quad i = 0, \ldots ,N - n - 2.$

В семействе $\mathcal{Y}_{i}^{{(n + 1)}}$ есть ровно ${{k}_{i}} = \left[ {(N - 1 - i){\text{/}}(n + 1)} \right]$ элементов. В частности, в последнем семействе лишь один элемент $\left| {\tilde {y}_{{N - n - 2,N - 1}}^{{(n + 1)}}} \right|$, а общее число элементов в наборе ${{\mathcal{Y}}^{{(n + 1)}}}$ равно

$\mathcal{K} = {{k}_{0}} + {{k}_{1}} + \ldots + {{k}_{{_{{N - n - 2}}}}}{\kern 1pt} .$

Пусть $\varkappa \in [1,\mathcal{K}]$ – порядковый номер элемента $\left| {\tilde {y}_{{i,j}}^{{(n + 1)}}} \right|$ в последовательном одномерном (линейном) наборе ${{\mathcal{Y}}^{{(n + 1)}}}$. Будем ассоциировать набор ${{\mathcal{Y}}^{{(n + 1)}}}$ с функцией

${{\mathcal{Y}}^{{(n + 1)}}}:[1,\mathcal{K}] \ni \varkappa \mapsto {{\mathcal{Y}}^{{(n + 1)}}}(\varkappa ){\kern 1pt} ,$

где ${{\mathcal{Y}}^{{(n + 1)}}}(\varkappa )$ есть элемент $\left| {\tilde {y}_{{i,j}}^{{(n + 1)}}} \right|$ набора ${{\mathcal{Y}}^{{(n + 1)}}}$ с номером $\varkappa $.

Шаг 2. Положим $\bar {\omega } = ({{\omega }_{1}},{{\omega }_{2}}, \ldots ,{{\omega }_{\mathcal{K}}})$, так, что $\mathcal{Y}({{\omega }_{\varkappa }}) \leqslant \mathcal{Y}({{\omega }_{m}})$, если $\varkappa < m$. Введем функцию

$\mathcal{M} = \mathcal{M}_{{\tilde {y}}}^{{(n + 1)}}:[1,\mathcal{K}] \ni \varkappa \mapsto \mathcal{M}_{{\tilde {y}}}^{{(n + 1)}}(\varkappa ) = {{\mathcal{Y}}^{{(n + 1)}}}({{\omega }_{\varkappa }}){\kern 1pt} .$

Пример с $\left\| \delta \right\| = 0$ при $n = 1$, $N = 10$ (тем самым, $\mathcal{K}$), $y:[0,1] \ni s \mapsto y(s) = 3{{s}^{3}}{\text{/}}2$. Тогда

$\left| {\tilde {y}_{{i,i + 2k}}^{{(n + 1)}}} \right| = \left| {y_{{i,i + 2k}}^{{(n + 1)}}} \right| = 9\frac{{i + k}}{{N - 1}} = i + k$ и на шаге 1 получаем

${{\mathcal{Y}}^{{(2)}}} = \left\{ {{{{(1,2,3,4)}}_{0}},\;{{{(2,3,4,5)}}_{1}},\;{{{(3,4,5,)}}_{2}},\;{{{(4,5,6)}}_{3}}{\kern 1pt} ,\;{{{(5,6)}}_{4}}{\kern 1pt} ,\;{{{(6,7)}}_{5}},\;{{{(7)}}_{6}},\;{{{(8)}}_{7}}} \right\}.$

Здесь нижний индекс $i$ у скобок соответствует семейству $\mathcal{Y}_{i}^{{(2)}}$. Шаг 2 в этом примере дает

$\begin{gathered} \bar {\omega } = (1,2,5,3,6,9,4,7,10,12,8,11,13,15,14,16,17,18,19,20)\quad {\text{и}} \\ \mathcal{M}_{{\tilde {y}}}^{{(2)}}(\varkappa ) \in \{ 1, 2,2, 3,3,3, 4,4,4,4, 5,5,5,5, 6,6,6, 7,7, 8\} . \\ \end{gathered} $

Шаг 3. Приблизим график функции $\mathcal{M}_{{\widetilde y}}^{{(n + 1)}}$ графиком кусочно-линейной функции

$k \mapsto {{V}_{{{{a}_{1}},{{b}_{1}},{{a}_{2}},{{b}_{2}},\varkappa }}}(k) = \left\{ \begin{gathered} {{a}_{1}}k + {{b}_{1}}\quad {\text{при}}\quad 1 \leqslant k \leqslant \varkappa , \hfill \\ {{a}_{2}}k + {{b}_{2}}\quad {\text{при}}\quad \varkappa < k \leqslant \mathcal{K}, \hfill \\ \end{gathered} \right.$

параметризованной числами ${{a}_{1}},\;{{b}_{1}},\;{{a}_{2}},\;{{b}_{2}},\;\varkappa $. Приближение достигается минимизацией

(8)

$S({{a}_{1}},{{b}_{1}},{{a}_{2}},{{b}_{2}},\varkappa ) = \sum\limits_{k = 1}^\varkappa {{\left( {\tilde {y}(k) - {{a}_{1}}k - {{b}_{1}}} \right)}^{2}} + \sum\limits_{k = \varkappa + 1}^\mathcal{K} {{\left( {k - \frac{{\tilde {y}(k) - {{b}_{2}}}}{{{{a}_{2}}}}} \right)}^{2}}\left( {\tilde {y}(k) - \tilde {y}(k - 1)} \right)$

по параметрам ${{a}_{1}}$, ${{b}_{1}}$, ${{a}_{2}}$, ${{b}_{2}}$, $\varkappa $.

Предложение 1. Для нахождения минимума функции (8) требуется линейно зависящее от $\mathcal{K}$ число операций сложения и умножения.

В самом деле, достаточно решить для каждого $\varkappa $ две задачи минимизации:

$\mathop {\min }\limits_{{{a}_{1}},{{b}_{1}}} \sum\limits_{k = 1}^\varkappa {{\left( {\tilde {y}(k) - {{a}_{1}}k - {{b}_{1}}} \right)}^{2}},\quad \mathop {\min }\limits_{{{a}_{2}},{{b}_{2}}} \sum\limits_{k = \varkappa + 1}^\mathcal{K} {{\left( {k - \frac{{\tilde {y}(k) - {{b}_{2}}}}{{{{a}_{2}}}}} \right)}^{2}}\left( {\tilde {y}(k) - \tilde {y}(k - 1)} \right).$

Если $\varkappa = \mathcal{K}$, то остается лишь первая из этих задач минимизации и вторая, если $\varkappa = 0$.

Рассмотрим 1-ю из этих задач для какого-либо $\varkappa = 1, \ldots ,\mathcal{K}$. Раскрывая выражение, стоящее под знаком суммы, и производя суммирование по $k$ получившихся слагаемых, получим

$ - 2{{a}_{1}}\sum\limits_{k = 1}^\varkappa \,k\tilde {y}(k) - 2{{b}_{1}}\sum\limits_{k = 1}^\varkappa \,\tilde {y}(k) + \sum\limits_{k = 1}^\varkappa \,\tilde {y}{{(k)}^{2}} + \frac{1}{6}a_{1}^{2}\varkappa (\varkappa + 1)(2\varkappa + 1) + {{a}_{1}}{{b}_{1}}\varkappa (\varkappa + 1) + b_{1}^{2}\varkappa .$

Поэтому минимум по ${{a}_{1}} = {{a}_{1}}(\varkappa )$, ${{b}_{1}} = {{b}_{1}}(\varkappa )$ в 1-й задаче минимизации полностью определяется значениями

а чтобы решить первые задачи нахождения минимума для всех $\varkappa $, нужно знать все $\alpha (\varkappa )$, $\beta (\varkappa )$, $\gamma (\varkappa )$, $\varkappa = 1, \ldots ,\mathcal{K}$. Так как

$\alpha (\varkappa ) = \alpha (\varkappa - 1) + \varkappa \tilde {y}(\varkappa ),\quad \beta (\varkappa ) = \beta (\varkappa - 1) + \tilde {y}(\varkappa ),\quad \gamma (\varkappa ) = \gamma (\varkappa - 1) + \tilde {y}{{(\varkappa )}^{2}},$

то для нахождения всех $\alpha (\varkappa )$, $\beta (\varkappa )$, $\gamma (\varkappa )$, $\varkappa = 1, \ldots ,\mathcal{K}$, нужно выполнить операций сложения и умножения в количестве, линейно зависящем от $\mathcal{K}$. Точно также и для вторых задач минимизации (по ${{a}_{2}},{{b}_{2}}$) находятся функции $\varkappa \mapsto {{a}_{2}}(\varkappa )$ и $\varkappa \mapsto {{b}_{2}}(\varkappa )$. Для их нахождения также требуется число сложений и умножений, линейно зависящее от $\mathcal{K}$.

Положив $S{\kern 1pt} *{\kern 1pt} (\varkappa )\;\mathop = \limits^{(8)} \;S({{a}_{1}}(\varkappa ),{{b}_{1}}(\varkappa ),{{a}_{2}}(\varkappa ),{{b}_{2}}(\varkappa ),\varkappa )$ и определив равенством

$S{\kern 1pt} *{\kern 1pt} \left( {\varkappa {\kern 1pt} *} \right)\;\mathop = \limits^{{\text{def}}} \;\mathop {\min }\limits_{\varkappa \in [1,\mathcal{K}]} S{\kern 1pt} *{\kern 1pt} (\varkappa ){\kern 1pt} ,$

число $\varkappa {\kern 1pt} *$, завершаем доказательство предложения 1.

По найденным параметрам $\varkappa {\kern 1pt} *$, ${{a}_{1}}\left( {\varkappa {\kern 1pt} *} \right)$, ${{b}_{1}}\left( {\varkappa {\kern 1pt} *} \right)$ задачи минимизации (8) определяется следующим образом искомое число $\hat {\varkappa }$, отделяющее “хвост” последовательности $\mathcal{M}(\varkappa )$. А именно,

$\hat {\varkappa } = \max \left\{ {\varkappa < \varkappa {\kern 1pt} *:\mathcal{M}(\varkappa ) < {{a}_{1}}\left( {\varkappa {\kern 1pt} *} \right)\varkappa + {{b}_{1}}\left( {\varkappa {\kern 1pt} *} \right)} \right\}.$

Итак, за приемлемую оценку $\left\| {{{{\tilde {y}}}^{{(n + 1)}}}} \right\|$ принимается $\mathcal{M}(\hat {\varkappa })$.

5. ЗАКЛЮЧЕНИЕ

Рассмотрена трудная и важная прикладная задача о вычислении производных априори гладкой функции по экспериментальным данным ее значений в дискретных точках. На практике этот вопрос обычно решается последовательным увеличением шага аппроксимации (например, в два раза) и этот процесс завершается тогда, когда графическое изображение разностной аппроксимации визуально будет мало отличаться от графика гладкой функции. Однако при этом есть немалый риск сгладить существенные особенности искомой производной. Кроме того, конечно, желательно иметь более обоснованный выбор шага аппроксимации и избавиться от необходимости визуального наблюдения. Обе эти негативные особенности отсутствуют в алгоритме, предложенном в данной работе. Это позволило дать ответ на вопрос Ф.И. Атауллаханова, касающийся проблемы определения концентрации тромбина по данным интенсивности введенного в кровь фермента, активирующего ее свертываемость. Суть проблемы в том, что искомая концентрации тромбина задается правой частью уравнения диффузии, в левой части которого стоит разность производной первого порядка по времени и производной второго порядка по пространственной переменной интенсивности фермента, которая определяется экспериментально.

Список литературы

Dunster J.L., Gibbins J.M., Panteleev M.A., Volpert V. Modeling thrombosis in silico: Frontiers, challenges, unresolved problems and milestones // Physics of Life Reviews. 2018. Vol. 26–27. P. 57–95. https://doi.org/10.1016/j.plrev.2018.02.005
Panteleev M.A., Dashkevich N.M., Ataullakhanov F.I. Hemostasis and thrombosis beyond biochemistry: roles of geometry, flow and diffusion // Thrombosis Research. 2015. Vol. 136. No 4. P. 699–711. Epub 2015 Jul 29. Review. PubMed PMID: 26278966.https://doi.org/10.1016/j.thromres.2015.07.025
Атауллаханов Ф.И., Лобанова Е.С., Морозова О.Л., Шноль Э.Э., Ермакова Е.А., Бутылин А.А., Заикин А.Н. Сложные режимы распространения возбуждения и самоорганизация в модели свертывания крови // Успехи физ. наук. 2007. Т. 177. № 1. С. 87–104.
Арестов В.В., Акопян Р.Р. Задача Стечкина о наилучшем приближении неограниченного оператора ограниченными и родственные ей задачи // Тр. Ин-та матем. и механ. УрО РАН. 2020. Т. 26. № 4. С. 7‒31.
Стечкин С.Б. Неравенства между нормами производных произвольной функции // Acta scient. math. 1965. Vol. 26. № 3–4. P. 225–230.
Арестов В.В. О наилучшем приближении операторов дифференцирования // Матем. заметки. 1967. Т. 1. № 2. С. 149–154.
Буслаев А.П. О приближении оператора дифференцирования // Матем. заметки. 1981. Т. 29. № 5. С. 372–378.
Бабенко К.И. Основы численного анализа. Москва-Ижевск: НИЦ “Регулярная и хаотическая динамика”, 2002. С. 3–847.

Дополнительные материалы отсутствуют.

Инструменты

следующая статья выпуска предыдущая статья выпуска содержание выпуска

Журнал вычислительной математики и математической физики

Архивы выпусков Информация о журнале Отправить рукопись в журнал