Консервативная монотонизация варианта CIP схемы для решения уравнения переноса

Аристова Е.Н.; Караваева Н.И.

doi:https://doi.org/10.20948/prepr-2020-121

№ 121 за 2020 г.

https://doi.org/10.20948/prepr-2020-121

Консервативная монотонизация варианта CIP схемы для решения уравнения переноса

Для интерполяционно-характеристического метода решения уравнения переноса, основанного на эрмитовой интерполяции, исследованы условия нарушения монотонности. Эрмитова интерполяция строится на предыдущем временном слое с использованием сеточных значений самой функции и ее пространственных производных. Замыкание алгоритма производится с использованием интегральных средних и формулы Эйлера–Маклорена. В случае, когда производные на концах отрезка интерполяции имеют одинаковый знак, предложен консервативный вариант монотонизации, сохраняющий интегральные средние и значения функции в узлах. В случае, когда производные имеют разные знаки, сохранить узловые значения функции не удается. Здесь приходится использовать кусочно-линейную аппроксимацию.

уравнение переноса, сеточно-характеристический метод, эрмитова интерполяция, CIP метод, консервативная монотонизация.

Сonservative limiting method for modification of the CIP method for solving the transport equation

Elena Nikolaevna Aristova
Nataliia Igorevna Karavaeva

The monotonicity breaking conditions are investigated for numerical solving of transport equation by grid-characteristic method based on the Hermitian interpolation. Hermitian interpolation is built on the previous time layer using the grid values of the function itself and its spatial derivatives. The algorithm is closed using integral averaged values and the Euler–Maclaurin formula. In the case of the same signs of the derivatives at the ends of the interpolation segment, a conservative modification of limiting method is proposed, which preserves the integral averaged values and the nodal values of the function. In the case of different signs of the derivatives, it is not possible to preserve the nodal values of the function. In this case piecewise linear approximation is used.

transport equation, grid-characteristic method, Hermitian interpolation, CIP method, conservative limiting method.

Работа выполнена при поддержке Российского фонда фундаментальных исследований, проект 18-01-00857-а.

1. Введение

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

Интерполяционно-характеристические методы решения уравнения переноса основаны на интерполяции сеточной функции и интегрировании уравнения переноса вдоль обратной характеристики от точки ее пересечения с границами ячейки до заданной точки нового временного слоя. Использование компактного двухточечного шаблона в простейшем варианте позволяет построить только линейную аппроксимацию, и, соответственно, получаются явные уголки. Одним из способов повышения порядка аппроксимации с сохранением компактного шаблона является использование эрмитовой интерполяции. Эрмитова интерполяция строится на основании узловых значений функции и ее первых пространственных производных. Метод CIP (Cubic Interpolation Polinom) третьего порядка аппроксимации по обеим переменным был предложен около тридцати лет назад и хорошо разработан [1-7]. В нем для нахождения производных на новом временном слое используется продолженное уравнение переноса, т.е. уравнение, записанное для пространственной производной. Однако такой способ замыкания не слишком удобен при решении неоднородного уравнения переноса с переменным коэффициентом поглощения. Во-первых, решения уравнения переноса могут быть разрывными, и, во-вторых, при переменном коэффициенте поглощения в продолженное уравнение войдет и исходная функция, что усложняет решение данного уравнения. Поэтому в работе [8] была предложена модификация метода CIP. Она основана на использовании интегральных средних и формулы Эйлера–Маклорена, связывающей интегральные средние в ячейке, узловые значения и узловые значения производных. Использование интегральных средних и формулы Эйлера–Маклорена делает данную схему родственной бикомпактным схемам Рогова [9-15], однако в бикомпактных схемах построения делаются методом прямых, что позволяет строить независимые аппроксимации по времени и пространству. Для сеточно-характеристических схем порядки аппроксимации по времени и пространству жестко связаны. Восстановление производных из интегральных средних приводит к консервативности предлагаемой модификации метода, основанного на эрмитовой интерполяции.

Исследование диссипативно-дисперсионных свойств ряда компактных схем показало, что схема CIP обладает относительно небольшой диссипаций и экстрамалой дисперсией [16]. Однако схемы высокого порядка немонотонны и, как следствие, неположительны. Это свойство является очень нежелательным при численном решении уравнения переноса. В данной работе рассматривается консервативная монотонизация предложенного в [8] метода, основанного на эрмитовой интерполяции.

2. Постановка задачи

Рассмотрим задачу нахождения численного решения нестационарного уравнения переноса незаряженного излучения или нейтральных частиц в плоской одномерной геометрии вида

L u (x, t) \equiv \frac{\partial u (x, t)}{\partial t} + μ \frac{\partial u (x, t)}{\partial x} + κ (x, t) u (x, t) = f (x, t)

.

(1)

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

u (x, 0) = u_{0} (x), 0 \leq x \leq X

.

(2)

Рассмотрим случай положительной скорости переноса $μ$ . Тогда поставим на левой границе классические граничные условия вида

u (0, t) = ϕ (t), 0 \leq t \leq T

.

(3)

Для начала приведем подход к построению схемы для однородного уравнения переноса

L u (x, t) \equiv \frac{\partial u (x, t)}{\partial t} + μ \frac{\partial u (x, t)}{\partial x} = 0

.

(4)

Интерполяционно-характеристическая схема будет построена в рамках одной расчетной ячейки, что упрощает ситуацию вблизи границ и контактных разрывов. Будем использовать аппроксимацию уравнения в рамках квадратной ячейки с вершинами $(x_{m}, t^{n}), (x_{m + 1}, t^{n}), (x_{m}, t^{n + 1}), (x_{m + 1}, t^{n + 1})$ . Для положительной скорости переноса $μ$ неизвестным сеточным значением искомой функции будет $y_{m + 1}^{n + 1}$ . Выпустим из точки $(x_{m + 1}, t^{n + 1})$ характеристику $x - μ t = c o n s t$ назад до пересечения либо с нижней, либо с боковой гранью ячейки (Рис. 1). В случае a) число Куранта $σ = \frac{μ τ}{h} \leq 1$ , в б) $σ \geq 1$ . Обозначим координату пересечения обратной характеристики с границами ячейки или $x *$ ( $σ \leq 1$ ), или $t *$ ( $σ \geq 1$ ). Для точного решения (1) значение из этой точки переносится без изменений в точку $(x_{m + 1}, t^{n + 1})$ . Таким образом, точность метода определяется точностью восстановления неузлового значения $y *$ .

Кубический интерполянт Эрмита для случая $σ \leq 1$ в барицентрических пространственных координатах $p$ и $q$ имеет вид

$P_{3} (p, q) = y_{m + 1}^{n} H^{R} (p, q) + y_{m}^{n} H^{L} (p, q) + d_{m + 1}^{n} h G^{R} (p, q) + d_{m}^{n} h G^{L} (p, q)$ ,

$p = \frac{x * - x_{m}}{h}, q = \frac{x_{m + 1} - x *}{h}$ , $p + q = 1$ , $h = x_{m + 1} - x_{m},$

$H^{R} = p (p + 2 p q), H^{L} = q (q + 2 p q)$ ,

$G^{R} = - p \cdot q p, G^{L} = q \cdot q p,$

(5)

где $y_{m}^{n}$ – сеточные значения искомой функции, а $d_{m}^{n}$ – сеточные значения пространственных производных искомой функции.

Случай $σ \geq 1$ рассматривается аналогично, интерполяция производится на левой границе ячейки с учетом производных по времени для искомой функции.

В предложенном варианте CIP метода не используются дифференциальные продолжения уравнения (1). Предположим, что в начальный момент времени мы знаем не только узловые значения функции $u_{0} (x_{m})$ , но и значения ее производной:

d_{m}^{0} = u_{0}^{'} (x_{m}) .

(6)

Для замыкания алгоритма нам необходима процедура получения сеточных значений $d_{m + 1}^{n + 1}$ . Для граничного узла значение производной в случае краевых условий (3) может быть рассчитано из (4) как $d_{0}^{n + 1} = - \dot{ϕ} (t^{n + 1}) / μ$ . Для однородного уравнения переноса производная $d_{m + 1}^{n + 1}$ может быть вычислена как производная интерполянта в точке $x *$ , однако такой вариант не обобщается на неоднородное уравнение переноса (1). Поэтому в [8] предлагается вычислить интегральное среднее на нижнем отрезке $[x *, x_{m + 1}]$ либо от интерполянта $P_{3}$ , либо по формуле Симпсона. В силу характеристических свойств однородного уравнения (4) это интегральное среднее будет совпадать с интегральным средним по отрезку $[t^{n}, t^{n + 1}]$ в точке $x_{m + 1}$ справа. Значение пространственной производной $d_{m + 1}^{n}$ может быть из уравнения (4) пересчитано в производную по времени: $g_{m + 1}^{n} = - μ d_{m + 1}^{n}$ . Тогда на правой границе нам известны два узловых значения $y_{m + 1}^{n}$ , $y_{m + 1}^{n + 1}$ , интегральное среднее ${\bar{y}}_{m + 1}$ и значение временной производной искомой функции $g_{m + 1}^{n}$ . Эти данные позволяют получить значение $g_{m + 1}^{n + 1}$ из формулы Эйлера–Маклорена:

{\bar{y}}_{m + 1} \equiv \frac{1}{τ} \int_{t^{n}}^{t^{n + 1}} u (x_{m + 1}, t) d t = \frac{1}{2} (y_{m + 1}^{n} + y_{m + 1}^{n + 1}) - \frac{τ}{12} (g_{m + 1}^{n + 1} - g_{m + 1}^{n})

.

(7)

В свою очередь, величина $g_{m + 1}^{n + 1}$ может быть пересчитана в пространственную производную $d_{m + 1}^{n + 1} = - g_{m + 1}^{n + 1} / μ$ из уравнения (4).

Интерполяционно-характеристический подход может быть распространен и на случай неоднородного уравнения (1). Заменой переменных вида

{\begin{cases} ξ = x - μ t \\ η = x + μ t \end{cases}

(8)

уравнение (1) может быть приведено к обыкновенному дифференциальному уравнению

2 μ \frac{\partial u}{\partial η} + κ u = f (\frac{ξ + η}{2}, \frac{η - ξ}{2 μ})

(9)

вдоль характеристики $ξ \equiv x - μ t = c o n s t$ . Решение этого неоднородного уравнения может быть найдено методом вариации постоянной, что позволяет получить решение в виде

u (η) = u * \exp (- \frac{κ (η - η *)}{2 μ}) + \frac{1}{2 μ} \int_{η *}^{η} \exp (- \frac{κ (η - η^{'})}{2 μ}) f (\frac{ξ + η^{'}}{2}, \frac{η^{'} - ξ}{2 μ}) d η^{'}

,

(10)

вычисляемое вдоль характеристики $ξ = c o n s t$ . Численное значение $y *$ может быть получено по формулам интерполяции Эрмита (5),(6), как и в однородном случае, только теперь для получения решения $y_{m + 1}^{n + 1}$ необходимо вычислить интеграл в (10) вдоль отрезка характеристики. Это может быть сделано, например, с помощью формулы Симпсона. Использование формулы Симпсона при умеренных оптических толщинах ячейки не ухудшает точность метода. Описание метода будет завершено, если мы укажем способ вычисления пространственной производной в точке $(x_{m + 1}, t^{n + 1})$ . В этом пункте численного метода, в отличие от однородного уравнения переноса, использование интегральных средних становится существенным. Интегральное среднее вдоль правого ребра ячейки от $(x_{m + 1}, t^{n})$ до $(x_{m + 1}, t^{n + 1})$ получается интегрированием интегрального среднего от $(x *, t^{n})$ до $(x_{m + 1}, t^{n})$ с преобразованием каждого значения подынтегральной функции типа (10) вдоль линейно меняющегося по длине отрезка характеристики. Для вычисления интегрального среднего по правому ребру ячейки будем использовать формулу Симпсона, а для вычисления $y_{m + 1}^{n + 1 / 2}$ , которое входит в формулу Симпсона, используем аналог формулы (10) по отрезку характеристики вдвое меньшей длины (Рис. 2).

Рис. 2. Конфигурация используемых характеристик для вычисления интегрального среднего на правом ребре расчетной ячейки.

Интегральное среднее на отрезке от $(x_{m + 1}, t^{n})$ до $(x_{m + 1}, t^{n + 1})$ позволяет по формуле (7) получить производную по времени $g_{m + 1}^{n + 1}$ , которая может быть пересчитана в пространственную с использованием уравнения (1). Тем самым, алгоритм получения решения с использованием эрмитовой интерполяции для неоднородного уравнения переноса завершен.

Монотонность метода определяется монотонностью построенного интерполянта.

3. Исследование монотонности эрмитовой интерполяции

Исследуем экстремумы функции (5). Индекс по времени для удобства будем опускать. И пусть для наглядности $d_{m}^{n} = f_{m}^{'}$ . Уравнение для нахождения экстремумов (5) будет определяться требованием $P_{3}^{'} (x) = 0$ и иметь вид

\frac{6 p q}{h} (y_{m + 1} - y_{m}) = 3 p q (f_{m + 1}^{'} + f_{m}^{'}) - p f_{m + 1}^{'} - q f_{m}^{'}

.

(11)

Уравнение (11) можно привести к виду

\frac{f_{m}^{'}}{p} + \frac{f_{m + 1}^{'}}{1 - p} = 6 (\frac{f_{m + 1}^{'} + f_{m}^{'}}{2} - \frac{y_{m + 1} - y_{m}}{h})

.

(12)

Изучим, когда уравнение (12) будет иметь корни на [0,1], что соответствует немонотонному профилю решения. При проведении такого исследования нам будет удобно использовать следующие обозначения

a = f_{m}^{'}, b = f_{m + 1}^{'}, c = 6 (\frac{f_{m + 1}^{'} + f_{m}^{'}}{2} - \frac{y_{m + 1} - y_{m}}{h})

.

(13)

Тогда уравнение принимает вид

\frac{a}{p} + \frac{b}{1 - p} = c

.

(14)

Исследование решений уравнения (14) представим как нахождение количества точек пересечения графиков функции $y (p) = \frac{a}{p} + \frac{b}{1 - p}$ и $\tilde{y} = c$ . Рассмотрим всевозможные варианты знаков величин $a$ и $b$ .

Вариант 1. Пусть $a b < 0$ , тогда на отрезке [0,1] функция $y (p)$ принимает значения от $- \infty$ до $+ \infty$ и корень у уравнения (14) есть при любом $c$ . Интерполянт является немонотонной функцией. В этом случае монотонизация проводится линейным профилем с сохранением $\bar{y}$ .

Вариант 2. Пусть $a b > 0$ , рассмотрим отдельно случаи положительных и отрицательных коэффициентов $a$ и $b$ .

Вариант 2.1. $a > 0$ и $b > 0$

Найдем минимальное значение, которое функция $y (p)$ принимает на [0,1]. Точки экстремума найдем из уравнения $y^{'} (p) = 0$ , которое примет вид

- \frac{a}{p^{2}} + \frac{b}{{(1 - p)}^{2}} = 0

.

(15)

Легко видеть, что ${(\frac{p - 1}{p})}^{2} = \frac{b}{a}$ , $1 - \frac{1}{p} = \pm \sqrt{b / a}$ . Тогда $\frac{1}{p} = 1 \mp \sqrt{b / a}$ и из двух корней мы выберем лежащий на [0,1], то есть $p_{\min} = \frac{\sqrt{a}}{\sqrt{a} + \sqrt{b}}$ . Найдем значение функции $y (p)$ в этой точке, получится $y (p_{\min}) = y_{\min} = {(\sqrt{a} + \sqrt{b})}^{2}$ .

Для получения условия монотонности остается сравнить $y_{\min}$ и $c$ . При $y_{\min} > c$ графики функций $y (p)$ и $\tilde{y} = c$ не пересекаются, решений у уравнения (14) нет, интерполянт монотонен. В обратной ситуации, когда $y_{\min} > c$ , имеем

6 (\frac{f_{m + 1}^{'} + f_{m}^{'}}{2} - \frac{y_{m + 1} - y_{m}}{h}) > {(\sqrt{f_{m}^{'}} + \sqrt{f_{m + 1}^{'}})}^{2}

.

(16)

У интерполянта существует точка минимума и максимума. Окончательно запишем условие немонотонности (16) в виде

\frac{y_{m + 1} - y_{m}}{h} < \frac{1}{3} (f_{m}^{'} + f_{m + 1}^{'} - \sqrt{f_{m}^{'} f_{m + 1}^{'}})

.

(17)

Это условие выполняется при небольшой разнице значений функции и больших значениях производных.

Вариант 2.2. $a < 0$ и $b < 0$

Рассматривается аналогично Варианту 2.1. Ищем точки экстремума функции $y (p)$ , но на этот раз нас интересует $p_{\max}$ и $y_{\max} = y (p_{\max})$ . Проведя аналогичные вычисления, получим, что $y_{\max} = - {(\sqrt{- a} + \sqrt{- b})}^{2}$ . И условие существования точек экстремума, то есть немонотонности, примет вид

\frac{y_{m + 1} - y_{m}}{h} > \frac{1}{3} (f_{m}^{'} + f_{m + 1}^{'} + \sqrt{f_{m}^{'} f_{m + 1}^{'}})

.

(18)

Подводя итог, скажем, что эрмитова интерполяция является монотонной при условиях

$f_{m}^{'} > 0, f_{m + 1}^{'} > 0$ и $\frac{y_{m + 1} - y_{m}}{h} \geq \frac{1}{3} (f_{m}^{'} + f_{m + 1}^{'} - \sqrt{f_{m}^{'} f_{m + 1}^{'}})$ или

$f_{m}^{'} < 0, f_{m + 1}^{'} < 0$ и $\frac{y_{m + 1} - y_{m}}{h} \leq \frac{1}{3} (f_{m}^{'} + f_{m + 1}^{'} + \sqrt{f_{m}^{'} f_{m + 1}^{'}})$ .

(19)

4. Монотонизация схемы при одинаковых знаках
узловых значений производной

Предложим способ осуществить монотонизацию разностной схемы, сохранив узловые значения $y_{m}, y_{m + 1}$ и величину интегрального среднего по ячейке $\bar{y}$ . Таким образом, остается возможность изменить значения производных, это соответствует изменению потоков, перетекающих из ячейки в ячейку, но не меняет общее количество «вещества» $y$ , которое является сохраняющейся величиной для уравнения адвекции.

Построим на $[x_{m}, x_{m + 1}]$ функцию $y = α {(x - x *)}^{3} + ω$ , которая является монотонной, выбрав параметры $α, ω, x *$ из требований (20-22)

y_{m} = α {(x_{m} - x *)}^{3} + ω

,

(20)

y_{m + 1} = α {(x_{m + 1} - x *)}^{3} + ω

,

(21)

\bar{y} = \frac{1}{h} \int_{x_{m}}^{x_{m + 1}} (α {(x - x *)}^{3} + ω) d x .

(22)

Для дальнейшего нам понадобится преобразовать (22), а именно

\begin{array}{l} \bar{y} = \frac{1}{h} \int_{x_{m}}^{x_{m + 1}} (α {(x - x *)}^{3} + ω) d x = \frac{1}{4 h} (α {(x_{m + 1} - x *)}^{4} - α {(x_{m} - x *)}^{4}) + ω = \\ = \frac{1}{4 h} ((y_{m + 1} - ω) (x_{m + 1} - x *) - (y_{m} - ω) (x_{m} - x *)) + ω = \\ = \frac{1}{4 h} (y_{m + 1} (x_{m + 1} - x *) - y_{m} (x_{m} - x *)) - \frac{ω}{4 h} (x_{m + 1} - x * - x_{m} + x *) + ω = \\ = y_{m + 1} \frac{x_{m + 1} - x *}{4 h} - y_{m} \frac{x_{m} - x *}{4 h} + \frac{3 ω}{4} . \end{array}

(23)

Мы получили форму уравнения (22), в которую входят только неизвестные $x *$ и $ω$ , причем линейно.

Теперь вернемся к уравнениям (20),(21). Найдем производную от каждого из них и вычтем друг из друга полученные выражения, чтобы получить связь $x *$ и $α$

y_{m + 1}^{'} - y_{m}^{'} = 3 α {(x_{m + 1} - x *)}^{2} - 3 α {(x_{m} - x *)}^{2} = 3 α (x_{m + 1} + x_{m} - 2 x *) h

.

(24)

Выразим $x *$ из (24)

x * = \frac{1}{2} (x_{m + 1} + x_{m}) - \frac{y_{m + 1}^{'} - y_{m}^{'}}{6 α h}

.

(25)

Отметим, что в силу формулы Эйлера–Маклорена $\bar{y} = \frac{1}{2} (y_{m + 1} + y_{m}) - \frac{h}{12} (y_{m + 1}^{'} - y_{m}^{'})$ и того, что в предлагаемом варианте монотонизации не терпят изменений величины $\bar{y}, y_{m + 1}, y_{m}$ , разность производных $y_{m + 1}^{'} - y_{m}^{'}$ сохраняется с точностью до $O (h^{4})$ .

Теперь с помощью (25) получим выражения для разностей $x_{m + 1} - x *$ и $x_{m} - x *$ , а именно

$x_{m + 1} - x * = \frac{h}{2} + \frac{y_{m + 1}^{'} - y_{m}^{'}}{6 α h}$ ,

$x_{m} - x * = - \frac{h}{2} + \frac{y_{m + 1}^{'} - y_{m}^{'}}{6 α h}$ .

(26)

Воспользуемся (26) и преобразуем (23)

\begin{array}{l} \bar{y} = y_{m + 1} \frac{1}{4 h} (\frac{h}{2} + \frac{y_{m + 1}^{'} - y_{m}^{'}}{6 α h}) - y_{m} \frac{1}{4 h} (- \frac{h}{2} + \frac{y_{m + 1}^{'} - y_{m}^{'}}{6 α h}) + \frac{3 ω}{4} = \\ = \frac{1}{8} (y_{m + 1} + y_{m}) + \frac{y_{m + 1}^{'} - y_{m}^{'}}{24 α h^{2}} (y_{m + 1} - y_{m}) + \frac{3 ω}{4} . \end{array}

(27)

Откуда

ω = \frac{4}{3} \bar{y} - \frac{1}{6} (y_{m + 1} + y_{m}) - \frac{y_{m + 1}^{'} - y_{m}^{'}}{18 α h^{2}} (y_{m + 1} - y_{m}) .

(28)

Теперь получим уравнение на $α$ . Вычтем из (21) (20), получится

y_{m + 1} - y_{m} = α {(\frac{h}{2} + \frac{y_{m + 1}^{'} - y_{m}^{'}}{6 α h})}^{3} - α {(- \frac{h}{2} + \frac{y_{m + 1}^{'} - y_{m}^{'}}{6 α h})}^{3}

.

(29)

После преобразований (29) примет вид

y_{m + 1} - y_{m} = \frac{1}{4} α h^{3} + \frac{1}{12 α h} {(y_{m + 1}^{'} - y_{m}^{'})}^{2}

.

(30)

Мы получили квадратное уравнение для нахождения $α$ . Придадим ему вид

α^{2} - \frac{4}{h^{3}} (y_{m + 1} - y_{m}) α + \frac{1}{3 h^{4}} {(y_{m + 1}^{'} - y_{m}^{'})}^{2} = 0

(31)

и найдем корни. Их количество будет определяться значением дискриминанта

D_{1} = \frac{4 {(y_{m + 1} - y_{m})}^{2}}{h^{6}} - \frac{{(y_{m + 1}^{'} - y_{m}^{'})}^{2}}{3 h^{4}}

.

(32)

При $D_{1} < 0$ предлагаемая монотонизация непригодна, в этом случае можно использовать линейную интерполяцию. При $D_{1} = 0$ получим 1 корень

α = \frac{2}{h^{3}} (y_{m + 1} - y_{m})

.

(33)

При $D_{1} > 0$ уравнение (31) имеет 2 корня

α = \frac{2}{h^{3}} (y_{m + 1} - y_{m}) \pm \sqrt{\frac{4 {(y_{m + 1} - y_{m})}^{2}}{h^{6}} - \frac{{(y_{m + 1}^{'} - y_{m}^{'})}^{2}}{3 h^{4}}}

.

(34)

Далее по (28) можно найти $ω$ и по (25) $x *$ . Таким образом находится весь набор неизвестных, через который определятся вид решения на $[x_{m}, x_{m + 1}]$ .

Наличие двух корней у уравнения (31) и, соответственно двух решений у исходной системы, привлекает внимание. На простом примере проиллюстрируем разницу в получаемых решениях.

Пусть $x_{m} = 0, x_{m + 1} = 1, h = 1, y_{m} = 0, y_{m + 1} = 1 / 3, y_{m}^{'} = 1, y_{m + 1}^{'} = 2.$ Посчитаем $\bar{y}$ по формуле Эйлера-Маклорена, запишем уравнение на $α$ при данном наборе параметров. Получим $\bar{y} = 1 / 12$ , $α^{2} - 4 α / 3 + 1 / 3 = 0$ . Отсюда $α_{1} = 1, α_{2} = 1 / 3$ . Далее по (28) получаем $ω_{1} = 1 / 27, ω_{2} = 0$ . Осталось посчитать $x *$ по (25). Получим $x_{1} * = 1 / 3, x_{2} * = 0$ . В итоге

\begin{array}{l} y_{1} = {(x - 1 / 3)}^{3} + 1 / 27, \\ y_{2} = x^{3} / 3. \end{array}

(35)

Вид функций $y_{1}, y_{2}$ иллюстрирует Рис. 3.

Graph1 — Рис. 3. График зависимости $y_{1} (x), y_{2} (x)$ , $P_{3} (x)$ .

Из двух вариантов мы предпочтем тот, в котором значения производных меняются меньше. В данном примере выбирается $y_{1}$ , так как $y_{1}^{'} (0) = \frac{1}{3}, y_{1}^{'} (1) = \frac{4}{3}$ , а $y_{2}^{'} (0) = 0, y_{2}^{'} (1) = 1$ .

Заключение

В работе кратко изложена консервативная модификация CIP метода решения уравнения переноса, основанная на интерполяции Эрмита для сеточно-характеристического метода. Проанализированы случаи, когда профиль интерполянта немонотонен, предложен вариант монотонизации решения. Изучены условия возможности проведения такой монотонизации. Это позволило получить монотонизированный профиль решения при совпадающих знаках производных в узлах ячейки, а также при выполнении ряда условий.