Проверка корректности решения сопряженной задачи

Зададим следующую пробную функцию: $g(x,t) = (1-2t) \sin(\pi x)$ и найдем её производные:

$\frac{\partial g}{\partial x} = \pi (1-2t) \cos(\pi x)$
$\frac{\partial^2 g}{\partial x^2} = - \pi^2 (1 - 2t) \sin(\pi x)$
$\frac{\partial g}{\partial t} = -2 \sin(\pi x)$.

Подставим найденные производные в исходное уравнение (см.Сопряженная задача) и перенеся все в правую часть, определим $g_d$.

\[g_d(x,t) = 2 \sin(\pi x) - ( - \varepsilon \pi^2 (1 - 2t) \sin(\pi x) )+ u(x,t) \pi (1 - 2t) \cos(\pi x) + q(x) (1 - 2t) \sin(\pi x) - 2 \delta( x - f_1(t) ) (u(x,t) - f_2(t))\]

Таким образом, пробная функция $g$ будет являться решением уравнения

\[\left\{ \begin{aligned} & \frac{\partial \psi}{\partial t} = - \varepsilon\frac{\partial^2 \psi}{\partial x^2} + u \frac{\partial \psi}{\partial x} + q(x)\,\psi -\\ & \qquad - 2\delta(x - f_1(t))(u(x,t) - f_2(t)) - g_d, \quad x \in (0,1), \quad t \in (0,T), \\ &\psi(0,t) = 0, \quad \psi(1,t) = 0, \quad t \in (0,T), \\ &\psi(x,T) = 0, \qquad x \in [0,1]. \end{aligned} \right.\]

Такая проверка корректности решения применяется в юнит тесте "tests/adjoint_check.jl". Файл содержит один @testset, внутри него реализовано решение вышеописанной системы, проверка его корректности через @test. А так же, @testset возвращает ψ, ψ_model, Xₙ, Tₘ, что соответствует решению, аналитическому решению, сетке по X, T.

Непосредственная реализация проверки сопряженной задачи

using Test;
using LaTeXStrings;
using Plots;

# Юнит тест проверяет корректность решения прямой задачи.
# Алгоритм описан в /docs/src/direct/adjoint_check.md
# Возвращает решение, аналитическое решение, сетку по X, сетку по T.

using NonLinearReactionAdvectionDiffusionWithFrontData;
using NonLinearReactionAdvectionDiffusionWithFrontData: phidetermination, Φ;
using NonLinearReactionAdvectionDiffusionWithFrontData: f1, f2;
using NonLinearReactionAdvectionDiffusionWithFrontData: heterogeneity, adjointRP, ∂ARP_∂y;

using ForwardDiff;

# Сначала, сгенирируем экспериментальные данные, на увеличенном числе узлов.
u_l(t) = -8
u_r(t) =  4
qf(x) = 4*sin(3 * π * x);       # Коэффициент линейного усиления, который в обратной
                                # задаче необходимо определить, но при генерации априорной
                                # информации мы задаем некоторый коэффициент, который,
                                # собственно, после имея априорную информацию и будем определять.
ε = 0.2;                        # Малый параметр при старшей производной
a, b = 0, 1;                    # Область по X
t₀, T = 0.0, 0.28;              # Область по T
N, M = 250, 500;                # Увеличенное Кол-во разбиений по X, T
h = (b-a)/N;                    # шаг по X
τ = (T-t₀)/M;                   # шаг по T
Xₙ = [a  + n*h for n in 0:N];   # Сетка по Х
Tₘ = [t₀ + m*τ for m in 0:M];   # Сетка по Т
qₙ =      qf.(Xₙ);              # Сеточные значения коэффициента линейного усиления
ulₘ=    u_l.(Tₘ);               # Сеточные значения левого  ГУ
urₘ=    u_r.(Tₘ);               # Сеточные значения правого ГУ
u₀ = u_init.(Xₙ);               # Начальные условия
#
u, XX, TP = solve(u₀, Xₙ, N, Tₘ, M, ε, ulₘ, urₘ, qₙ);

# Генерация априорной информации
ϕl      = phidetermination(qₙ, ulₘ, Xₙ, N, Tₘ, M);                  # Левый вырожденный корень
ϕr      = phidetermination(qₙ, urₘ, Xₙ, N, Tₘ, M, reverseX = true); # Правый вырожденный корень
ϕ       = Φ(ϕl, ϕr, N, M);                                          # Серединный корень
f1_data = f1(ϕ, u, Xₙ, N, M);                                       # Положение переходного слоя
f2_data = f2(f1_data, u, Xₙ, N, M);                                 # Значение функции на переходном слое

# Решение сопряженной задачи


# Теперь переходим непосредственно к проверке
# Модельная функция
raw"""
    g(n::Int, X = Xₙ, N::Int,
        m::Int, T = Tₘ, M::Int) -> Real

Модельная функция ``(1-2t)\sin(\pi x)``.
`n` — номер узла в сетке по X. `m` — номер шага в сетке по T.
`X` — сетка по X, `N` — кол-во **интервалов** сетки.
`Tₘ` — сетка по T, `M` — кол-во **интервалов** сетки.
Захватывает переменную `T` конечное время моделирования из окружения.
"""
function g(n::Int, X, N::Int,
            m::Int, Tₘ, M::Int)
    t = Tₘ[m]
    x = X[n]
    return (1 - 2t/T)*sin(π*x)
end

# Выберем априорный параметр аппроксимации дельта-функции
ω = 0.0001;
raw"""
    g_d(n::Int, Xₙ::Vector, N::Int,
        m::Int, Tₘ::Vector, M::Int,
        ε::Real, qₙ::Vector,
        u::Matrix, f1::Vector, f2::Vector) -> Real

Вычисляет невязку, т.е. рез ультат подстановки ``g`` в постановку сопряженной задачи.
`n` — номер узла в сетке по X. `m` — номер шага в сетке по T.
`X` — сетка по X, `N` — кол-во **интервалов** сетки.
`Tₘ` — сетка по T, `M` — кол-во **интервалов** сетки.
"""
function g_d(n::Int, Xₙ::Vector, N::Int,
                m::Int, Tₘ::Vector, M::Int,
                ε::Real, qₙ::Vector,
                u::Matrix, f1::Vector, f2::Vector,
               w::Real)
    x = Xₙ[n];
    t = Tₘ[m];
    out  = 2/T * sin(π * x) - ( - ε * π^2 * (1 - 2t/T) * sin(π * x)) + π * (1 - 2t/T) * cos(π * x) * u[n,m] + qₙ[n] * (1 - 2t/T) * sin(π * x) - heterogeneity(n, m, Xₙ[2:N], N, Tₘ, M, u, f1, f2, w)
    return out;
end

Uₙₘ = u;                                        # Сохраним матрицу решения прямой задачи
ψl = [0.0 for i in 1:M+1];                      # Левые  ГУ для модельной функции
ψr = [0.0 for i in 1:M+1];                      # Правые ГУ для модельной функции
ψ₀ = [g(n, Xₙ, N, M+1, Tₘ, M) for n in 1:N+1];  # Начальные условия для модельной функции
# Внимание! сопряженная задача — ретроспективная
# ``u₀ = g(x, T)``
# Модельное решение найденное с помощью известного аналитического решения
# Для его генерации используем инвертированную сетку по времени, чтобы
# массив `ψ_model` и `ψ` имели одно и тоже направление хода времени.
ψ_model = [ g(n, Xₙ, N, m, Tₘ, M) for n in 1:N+1, m in M+1:-1:1];

# Создадим функцию, которая будет вычислять вектор правой части с добавлением невязки
raw"""
    RP(y, m, Xₙ, N, Tₘ, M, ε, qₙ, u, f1, f2) -> Vector

# Return
    Вектор размера `N-1`, сеточные значения правой части уравнения,
    для которого `g(x,t)` будет являться решением.

!!! warning
    Функция не предназначена для самостоятельного использования,
    она передается в качестве аргумента в `solve`, внутри которой
    для нее сформуются аргументы нужной длины, для которых не
    потребуется смещений индексов.
"""
function RP(y, m, Xₙ, N, Tₘ, M, ε, ψl, ψr, qₙ, u, f1, f2, w)
    arp = adjointRP(y, m, Xₙ, N, Tₘ, M, ε, ψl, ψr, qₙ, u, f1, f2, w)
    d = [ g_d(n, Xₙ, N, m, Tₘ, M, ε, qₙ, u, f1, f2, w) for n in 1:N-1]
    return arp .- d
end
# Якобиан сконструируем тем же способом, что и внутри пакета в `src/adjoint.jl`.
j(y, m, Xₙ, N, Tₘ, M, ε, ψl, ψr, qₙ, U, f1, f2, w) = ForwardDiff.jacobian( z -> adjointRP(z, m, Xₙ, N, Tₘ, M, ε, ψl, ψr, qₙ, U, f1, f2, w), y)

# С использованием автоматического дифференцирования
ψ = solve_adjoint(ψ₀, Xₙ, N, Tₘ, M, ε, ψl, ψr, qₙ, Uₙₘ, f1_data, f2_data, RP, j, w = ω)
@test all(isapprox.(ψ_model, ψ, atol = 0.025))

# С использованием трехдиагонального якобиана
ψ = solve_adjoint(ψ₀, Xₙ, N, Tₘ, M, ε, ψl, ψr, qₙ, Uₙₘ, f1_data, f2_data, RP, ∂ARP_∂y, w = ω)
@test all(isapprox.(ψ_model, ψ, atol = 0.025))

(ψ, ψ_model, Xₙ, Tₘ)

Посмотрим на результат решения в сравнении с аналитическим решением

d = [missing, missing];
dd = [missing missing; missing missing];
make_gif(ψ, Xₙ, Tₘ, dd, dd, dd, d, d, ψ_model; name = "direct_check.gif")

Найдем абсолютную погрешность численного решения

err = ψ .- ψ_model
heatmap(Xₙ, Tₘ, err', xlabel=L"X_n", ylabel=L"T_m", title="Absolute Error", size=(1200, 800))

This page was generated using Literate.jl.