Справочник

module NonLinearReactionAdvectionDiffusionWithFrontData

Пакет реализует исследование в рамках гранта РФФИ #20-31-70016 "Численные методы решения обратных задач для нелинейных сингулярно возмущённых уравнений типа реакция-диффузия-адвекция с данными о положении фронта реакции".

Автор исходного кода пакета:

• Андрей Борзунов, Кафедра математики физического факультета МГУ им. Ломоносова.

Исполнители гранта:

• Мельникова Алина Александровна
• Левашова Наталия Тимуровна
• Лукьяненко Дмитрий Витальевич (Руководитель)
• Быцюра Светлана Владимировна
• Борзунов Андрей Анатольевич
• Исаев Темур Фуркатович
• Аргун Рауль Ларикович
• Горбачев Александр Викторович

Документация: https://github.com/aborzunov/NonLinearReactionAdvectionDiffusionWithFrontData.jl

J(uˢ::Matrix, Xₙ::Array, N::Int,
  Tₘ::Vector, M::Int,
  f1::Vector, f2::Vector,
  qₙˢ::Vector, α::Real = 0.0) -> Real

Вычисляет функционал $J(\mathbf{x}) = \int_0^T \left( u(f_1(t), t; q^s) - f_2(t) \right)^2 + \alpha \int_0^1 q^2(x) dx$.

Использует f2 для вычисления $u(f_1(t), t)$ и после по формуле трапеций.

J_q(uˢ::Matrix, ψˢ::Matrix,
    Xₙ::Array, N::Int,
    Tₘ::Vector, M::Int) -> Vector

Вычисляет градиент $J_q(x) = \int_0^T u^s(x, t) \psi^s(x, t) dt$.

apply_on_dynamic_mesh(ϕ::Matrix, XX::Matrix,
                  N::Int, M::Int) -> Matrix

Аппроксимирует ϕ на динамичски изменяющуюся на каждом временном шаге сетку XX.

Функция phidetermination определяют вырожденные корни на стартовой сетке. Φ проиозводит линейные операции над векторами на каждом временном шаге, опять же на стартовой сетке. Поэтому если в расчетах была использована динамическая сетка, то нужно переопределить ϕ на этой динамической сетке.

delta(x, Xₙ, x₀ = 0) -> ::Real

Вычисляет конечно разностную аппроксимацию дельта функции $\delta(x; x₀)$ на сетке Xₙ исходя из $\int\limits_a^b \delta(x; x₀) dx = 1$: если $x, x_0 \in [X_n, X_{n+1}]$, возвращаем $1/h$, где $h$ –- шаг сетки.

dparams()

Возвращает стандартный набор параметров для равномерной сетки.

Return

`a, b, t₀, T, N, M, ε, Xₙ, Tₘ, qₙ, ulₘ, urₘ, u₀ `

f1(ϕ::Matrix, u::Matrix, Xₙ::Array, N::Int, M::Int) -> Vector

Находит положение переходного слоя $f_1(t) = x_{t.p.}(t)$ путем поиска точки пересечения аргумента, при котором пересекаются $u(x,t)$ и $ϕ(x)$.

Точка пересечения находится путем интерполяции функции обратной к $u(x,t) - ϕ(x) = 0$, передавай аргументы в инвертированном виде в find_f_zeros.

See also: find_f_zeros

f2(ϕ::Matrix, u::Matrix, Xₙ::Array, N::Int, M::Int) -> Vector

Находит значение искомой функции на переходном слое $f_2(t) = u(x_{t.p.}, t)$. Находится путем интерполяции функции $u(x - f1(t), t) = 0$.

See also: find_f_zeros

generate_obs_data(u::Matrix, Xₙ::Vector, N::Int,
                  Tₘ::Vector, M::Int,
                  qₙ::Vector,
                  ulₘ::Vector, urₘ::Vector)

Фнукция-сокращение.

#Return Левый вырожденный корень, правый, их полуразность, положение переходного слоя, значение u на переходном слое. ϕl, ϕr, ϕ, f1_data, f2_data

See also: phidetermination, Φ, f1, f2.

make_gif(u::Matrix, Xₙ::Vector, Tₘ::Vector,
              ϕ_l::Matrix = missings(2), ϕ_r::Matrix = missings(2),
              f1::Vector = missings(2), f2::Vector = missings(2),
              analitical = nothing;
              frames_to_write::Int = -1,
              name = "solution.gif",
              convert2mp4 = false)

Рисует gif анимацию решения. вплоть по frames_to_write-ый кадр, сохраняет как "results/name".

Так же рисует аналитическое решение analitic(x,t), если таково передано. label::String — LaTeX строка подписи искомой функции с экранирование спец. символов.

TODO: Fix doc

make_minimzation_gif(Js::Vector, Qs::Matrix,
                     qₙ::Vector, Xₙ::Vector;
                     frames_to_write::Vector = Vector(),
                     name = "solution.gif", convert2mp4 = false,
                     β::Real = 0.0)

Сохраняет анимацию процесса минимизации.

make_plot(u::Matrix, Xₙ::Vector, Tₘ::Vector, m::Int,
               ϕ_l::Matrix = missings(2), ϕ_r::Matrix = missings(2),
               f1::Vector = missings(2), f2::Vector = missings(2),
               analitical = nothing;
               label::String = "u")

Рисует m-ый кадр решения u. Xₙ, Tₘ — сетки. ϕ_l, ϕ_r — вырожденные решения. f1, f2 — сгенерированные априорные данные. analitical — или функция или сеточные значения аналитического решения. label::String — LaTeX строка подписи искомой функции с экранирование спец. символов.

minimize(q₀::Vector, u₀::Vector,
         ulₘ::Vector, urₘ::Vector,
         Xₙ, N,
         Tₘ::Vector, M,
         ε,
         f1_data::Vector, f2_data::Vector;
         S::Int = 10,
         β::Real = 0.01,
         w::Real = 0.0001,
         create_mesh::Function = x -> [NaN, NaN])  -> Vector, Vector, Matrix

Вернет qˢ на оригинальной сетке Xₙ[:,1].

phidetermination(q::Vector, ub::Vector,
                 Xₙ::Vector, N::Int,
                 Tₘ::Vector, M::Int;
                 reverseX = false)

Решает ОДУ для нахождения вырожденного корня.

reverseX флаг обозначающий обратное направление интегрирования по оси $X$.

solve(y₀::Vector, Xₙ::Vector, N::Int,
      Tₘ::Vector, M::Int,
      ε::Real, ulₘ::Vector, urₘ::Vector,
      qₙ::Vector,
      RP::Function = directRP,
      jac::Function = ∂directRP_∂y;
      α::Complex = complex(0.5, 0.5),
      create_mesh::Function = x -> [NaN, NaN]) -> Matrix, Matrix, Vector

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

Arguments

• y₀::Vector: Сеточные значения начального условия.
• Xₙ::Vector: Пространственная сетка по $X$.
• N::Int: Число интервалов сетки по $X$.
• Tₘ::Vector: Пространственная сетка по $t$.
• M::Int: Число интервалов сетки по $t$.
• ε::Real: Малый параметр при старшей производной.
• ulₘ::Function: Сеточные значения левого ГУ.
• urₘ::Function: Сеточные значения правого ГУ.
• qₙ::Vector: Сеточные значения "неоднородности", см. постановку задачи.
• RP::Function: Функция вычисления вектора правой части — directRP.
• jac::Function: Якобиан правой части по вектору y — ∂DRP_∂y.
• α::Complex: Коэффициент схемы. При α = 0 — схема Эйлера, при α = complex(0.5, 0.5) — схема Розенброка с комплексным коэффициентом.
• create_mesh::Function: Принимает один аргумент x_tp положение переходного слоя и формирует подходяющую сетку.

Return

Тройку: - Матрицу размера $(N+1, M+1)$, содержащую значения искомой функции на сетках $X\_n, T\_m$,. - Матрицу размера $(N+1, M+1)$, содержащую значения стеки $X\_n$ на каждом моменте времени. - Вектора размера $M+1$, содержащий координаты положения переходного слоя в каждый момент времени.

u_init(x::Real; ε::Real = 0.2, x_tp = 0.15) -> Real

Начальные условия в виде $(x^2 - x -2) -6 \tanh( -3 \xi)$, где $\xi = \frac{x - x_tp}{ε}$. x_tp — положение внутрннего переходного слоя.

Example

julia> NonLinearReactionAdvectionDiffusionWithFrontData.u_init.(0:0.1:1)
11-element Array{Float64,1}:
 -7.868156688432881
 -5.900893714323723
  1.650893714323725
  3.6581566884328813
  3.753366656356917
  3.749669571706625
  3.7599835485135173
  3.7899991809276505
  3.839999959220786
  3.909999997969723
  3.999999999898918

Φ(ϕ_l::Matrix, ϕ_r::Matrix, N::Int, M::Int) -> ::Vector

Вычисляет полуразность вырожденных решений $|\phi_l^m - \phi_r^m|/2 + \phi_l^m$ на каждом шаге по времени m с помощью матриц вырожденных решений $\phi_l$ и $\phi_r$ вырожденного корня.

" ARP_y(y::Vector, m::Int, Xₙ::Vector, N::Int, Tₘ::Vector, M::Int, ε::Real, ulₘ::Vector, urₘ::Vector, qₙ::Vector, Uₙₘ::Matrix, f1::Vector, f2::Vector) -> (::Vector, ::Vector, ::Vector)

Возращает три диагонали якоибана правой части сопряженной задачи. Входные параметры полностью аналогичны adjointRP.

Return

Возвращает три вектора dl, d, dl элементов якоибана. Нижняя и верхняя диагонали длины $N-2$, главная — $N-1$.

    [ d[1]  du[1]                         ]
    [ dl[1] d[2]  du[2]                   ]
    [ 0     dl[2] d[3] du[3]              ]
    [           ...  ...  ...             ]
    [                ...  ...     du[N-2] ]
    [                     dl[N-2] d[N-1]  ]

" DRP(y::Int, m::Int, Xₙ::Vector, N::Int, ε::Real, ulₘ::Vector, urₘ::Vector, qₙ::Vector) -> (::Vector, ::Vector, ::Vector)

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

Return

Возвращает три вектора dl, d, dl элементов якоибана. Нижняя и верхняя диагонали длины $N-2$, главная — $N-1$.

    [ d[1]  du[1]                         ]
    [ dl[1] d[2]  du[2]                   ]
    [ 0     dl[2] d[3] du[3]              ]
    [           ...  ...  ...             ]
    [                ...  ...     du[N-2] ]
    [                     dl[N-2] d[N-1]  ]

Example

julia> # Данные необохдимо подготовить, например, как в examples/example_direct.jl

julia> # Подготовим массивы, выбросив граничные точки

julia> # ведь тестируемая функция — для внутреннего использования

julia> qq = qₙ[2:N];

julia> y = y₀[2:N];

julia> dl, d, du = NonLinearReactionAdvectionDiffusionWithFrontData.DRP_y(y, 1, Xₙ, N, ε, ulₘ, urₘ, qq);

julia> Tridiagonal( dl, d, du )
49×49 Tridiagonal{Float64,Array{Float64,1}}:
 -991.038   305.462      ⋅         ⋅         ⋅         ⋅         ⋅     …       ⋅         ⋅         ⋅         ⋅         ⋅        ⋅
[...]

adjointRP(y::Vector, m::Int,
          Xₙ::Vector, N::Int,
          Tₘ::Vector, M::Int,
          ε::Real,
          ulₘ::Vector, urₘ::Vector,
          qₙ::Vector,
          Uₙₘ::Matrix, f1::Vector, f2::Vector)

Вычисляет вектор правой части с помощью конечно-разностной аппроксимации пространственных производных.

Arguments

• y::Vector: Вектор размера N-1, решение в текущий момент без ГТ.
• m::Int: Номер шага в сетке по времени.
• Xₙ::Vector: Вектор размера N+1 сетки по $X$, вместе с ГТ.
• N::Int: Число интервалов в полной сетке по $X$.
• Tₘ::Vector: Вектор размера M+1 сетки по $T$.
• M::Int: Число интервалов в полной сетке по $T$.
• ε::Real: Малый параметр при старшей производной.
• ulₘ::Vector: Вектор сеточных значений левого ГУ.
• urₘ::Vector: Вектор сеточных значений правого ГУ.
• qₙ::Vector: Вектор размера N-1 сеточные значений неоднородности без ГТ.
• Uₙₘ::Matrix: Матрица размера N-1, M+1 решения прямой задачи при данном qₙ.
• f1::Vector: Эспериментальные данные — положение внутреннего слоя, размера M+1.
• f2::Vector: Эспериментальные данные — значение функции на внутреннеем слое, размера M+1.
• w::Real: Априорный параметр в аппроксимации дельта-функции (см. heterogeneity, δw).

deltaw(n::Int, x₀::Real, Xₙ::Vector, N::Int, w::Real) --> ::Real

Вычисляет $\delta(x - x₀)$, где x = Xₙ[n]. n — номер узла в сетке Xₙ(без граничных точек).

Использует более точную конечно разностную аппроксимацию дельта функции, w — априорный параметр, см. δw.

directRP(y::Vector, m::Int,
         Xₙ::Vector, N::Int,
         ε::Real,
         ulₘ::Vector, urₘ::Vector,
         qₙ::Vector)

Вычисляет вектор правой части с помощью конечно-разностной аппроксимации пространственных производных.

Arguments

• y::Vector: Вектор размера N-1, решение в текущий момент без ГТ.
• m::Int: Номер шага в сетке по времени.
• Xₙ::Vector: Вектор размера N+1 сетки по $X$, вместе с ГТ.
• N::Int: Число интервалов в полной сетке по $X$.
• ε::Real: Малый параметр при старшей производной.
• ulₘ::Vector: Вектор сеточных значений левого ГУ.
• urₘ::Vector: Вектор сеточных значений правого ГУ.
• qₙ::Vector: Вектор размера N-1 сеточные значений неоднородности без ГТ.

dparams_nonuniform()

Возвращает стандартный набор параметров для неравномерной сетки.

Return

return a, b, t₀, T, N, M, ε, Xₙ, Tₘ, qₙ, ulₘ, urₘ, u₀, meshf;

find_f_zeros(f::Vector, Xₙ::Vector)

Находит такой $x$, что $f(x) = 0$. f — сеточные значения функции на сетке Xₙ.

heterogeneity(n::Int, m::Int,
             Xₙ::Vector, N::Int,
             Tₘ::Vector, M::Int,
             Uₙₘ::Matrix,
             f1::Vector, f2::Vector,
             w::Real)

Неоднородность выражающая невязку текущего решения с искомым. $2 \delta( x - f_1(t) ) ( u^s(x,t) - f_2(t) )$.

heterogeneity_map(XX::Matrix, N::Integer,
                  Tₘ::Vector, M::Integer,
                  u::Matrix,
                  f1_data::Vector, f2_data::Vector,
                  w::Real) --> ::Matrix

Возвращает матрицу неоднородности $-\delta(x - f_1(t))(u - f_2(t))$. Следует использовать для контроля выбора параметра аппроксимации $w$ в дельта-функции.

shishkin_mesh(a::Real, b::Real,
              x_tp::Real, ε::Real,
              N::Int = 50,
              C_i::Real = 1.0, K_i::Real = 1.0,
              C_b::Real = 1.0, K_b::Real = 1.0) --> Vector

Возвращает кусочнораномерную сетку со сгущением на границах и на переходном слое.

Экскиз:

.....  .  .  .  .  .  .  .  .  .  .......  .  .  .  .  .  .  .  .  .  .  ......

Arguments

• a::Int Левая граница.
• b::Int Правая граница.
• x_tp::Real Положение переходного слоя.
• ε::Real Малый параметр при старшей производной.
• C_i::Real Масштабирующий коэффициент для ширины внутреннего сгущения.
• K_i::int Масштабирующий коэффициент количества интервалов внутреннего сгущения.
• C_b::Real Масштабирующий коэффициент для ширин пограничных сгущений.
• K_b::Int Масштабирующий коэффициент количеств интервалов пограничных сгущений.
• K::Int Кол-во интервалов вне всех сгущений.

strip_borderPoints(a::Vector, N) -> Vector

Функция для внутренного использования.

Входящий массив должен быть размера N+1. Обрезает граничные точки слева и справа. Возвращает массив размера N-1.

Example

julia> N = 10; a = collect(1:N+1);

julia> NonLinearReactionAdvectionDiffusionWithFrontData.strip_borderPoints(a, N)
9-element Array{Int64,1}:
  2
  3
  4
  5
  6
  7
  8
  9
 10

δw(x::Real, ω::Real) -> Real

Линейная аппроксимация $δ(x)$, ω — эвристический парамет. Подбирается так, чтобы невязка сопряженной задачи $2 \delta( x - f_1(t)) ( u(x,t) - f_2(t))$ при подстановке в неё решения u(x,t;q) при истинном q, обнулялась почти везде, но не всюду.

∂ARP_∂y(y::Vector, m::Int,
      Xₙ::Vector, N::Int,
      Tₘ::Vector, M::Int,
      ε::Real,
      ulₘ::Vector, urₘ::Vector,
      qₙ::Vector,
      Uₙₘ::Matrix, f1::Vector, f2::Vector) -> Tridiagonal

Обертка фнукции ARP_y, которая возвращает Tridiagonal( ARP_y(...)) трехдиагональную матрицы из векторов, которые возвращает ARP_y.

∂DRP_∂y(y::Vector, m::Int,
        Xₙ::Vector, N::Int,
        ε::Real,
        ulₘ::Vector, urₘ::Vector,
        qₙ::Vector) -> Tridiagonal

Обертка фнукции DRP_y, которая возвращает Tridiagonal( DRP_y(...)) трехдиагональную матрицу из векторов, которые возвращает DRP_y(...).

∂adjointRP_∂y(y::Vector, m::Int,
               Xₙ::Vector, N::Int,
               Tₘ::Vector, M::Int,
               ε::Real,
               ulₘ::Vector, urₘ::Vector,
               qₙ::Vector,
               Uₙₘ::Matrix, f1::Vector, f2::Vector,
               w::Real)

Функция якобиана для adjointRP. Все аргументы, размерности входных векторов такие же, как и у adjointRP.

∂directRP_∂y(y::Vector, m::Int,
              Xₙ::Vector, N::Int,
              ε::Real,
              ulₘ::Vector, urₘ::Vector,
              qₙ::Vector)

Функция якобиана для adjointRP. Размерности входных векторов такие же, как и у directRP.