Примеры
На этой странице мы делаем include примеров, которые лежат в корневом каталоге пакета в папке "examples/".
Literate.jl включает текст этих примеров и в md файл.
Пример №1, подробный, статическая сетка
Сопряженная задача определяется для сеточного решения некоторой прямой задачи, и некоторых априорных данных. Поэтому, сначала необходимо сгенерировать априорные данные и некоторую матрицу u. Сделаем это на увеличенном числе интервалов, т.к. решение сопряженной задачи менее гладкое.
using NonLinearReactionAdvectionDiffusionWithFrontData
using NonLinearReactionAdvectionDiffusionWithFrontData: phidetermination, Φ;
using NonLinearReactionAdvectionDiffusionWithFrontData: f1, f2;
u_l(t) = -8
u_r(t) = 4
q(x) = 4*sin(3 * π * x); # Коэффициент линейного усиления, который в обратнойзадаче необходимо определить, но при генерации априорной информации мы задаем некоторый коэффициент, который, собственно, после имея априорную информацию и будем определять.
ε = 0.2; # Малый параметр при старшей производной
a, b = 0, 1; # Область по X
t₀, T = 0, 0.28; # Область по T
N, M = 50, 80; # Увеличенное Кол-во разбиений по 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ₙ = q.(Xₙ); # Сеточные значения коэффициента линейного усиления
ulₘ= u_l.(Tₘ); # Сеточные значения левого ГУ
urₘ= u_r.(Tₘ); # Сеточные значения правого ГУ
y₀ = u_init.(Xₙ); # Начальные условияu, XX, TP = solve(y₀, 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); # Нужно подать инвертированную сетку
ϕr = reverse(ϕr, dims=1); # А после — инвертировать решение по X
ϕ = Φ(ϕl, ϕr, N, M); # Серединный корень
f1_data = f1(ϕ, u, Xₙ, N, M); # Положение переходного слоя
f2_data = f2(f1_data, u, Xₙ, N, M); # Значение функции на переходном слое
# Решение сопряженной задачи
Uₙₘ = u; # Сохраним старую матрицу
y₀ = [0.0 for i in 1:N+1]; # Нулевые начальные условия
ψl = [0.0 for i in 1:M+1]; # Нулевые ГУ
ψr = [0.0 for i in 1:M+1]; # Нулевые ГУ
ψ = solve_adjoint(y₀, Xₙ, N, Tₘ, M, ε, ψl, ψr, qₙ, Uₙₘ, f1_data, f2_data, w = 0.0005)Визуализация
Мы не можем строить большие анимации на стороне Travis-a. Если мы на CI, то будем рисовать только 10 кадров. Если мы генерируем документацию локально, то рисуем 80 кадров.
isTravis = in("Travis", keys(ENV))
ftw = isTravis ? range(1, stop = M+1, length=9) : [1; 2:div(M+1, 80):M+1];На отрисовку, решение сопряженной задачи передадим в инвертированном времени. Передадим свою подпись к графиками с помощью keyword label="\\psi", не забыл про экранировку спец символа
make_gif(reverse(ψ, dims=2), XX, Tₘ;
label="\\psi",
name="solution_adjoint_ex1.gif",
frames_to_write = ftw)Результат должен быть около нулевой, ведь в качестве текущего приближения q мы взяли искомое, а при нем — градиент должен обнуляться.
Пример №3, подробный, динамическая сетка
# Сопряженная задача определяется для сеточного решения некоторой прямой задачи,
# и некоторых априорных данных. Поэтому, сначала необходимо сгенерировать
# априорные данные и некоторую матрицу `u`.
# Сделаем это на увеличенном числе интервалов, т.к. решение сопряженной задачи
# имеет более резкие фронты.
using NonLinearReactionAdvectionDiffusionWithFrontData
using NonLinearReactionAdvectionDiffusionWithFrontData: f1, f2;
using NonLinearReactionAdvectionDiffusionWithFrontData: shishkin_mesh;
using NonLinearReactionAdvectionDiffusionWithFrontData: phidetermination, Φ;
using NonLinearReactionAdvectionDiffusionWithFrontData: apply_on_dynamic_mesh;
u_l(t) = -8; # ГУ
u_r(t) = 4; #
qf(x) = 4*sin(3 * π * x); # Коэффициент линейного усиления
ε = 0.03; # Малый параметр при старшей производной
a, b = 0, 1; # Область по X
t₀, T = 0, 0.20; # Область по T
x_tp = 0.22; # Положение переходного слоя
M = 500; # Кол-во разбиений по X, T
τ = (T-t₀)/M; # шаг по T
Tₘ = [t₀ + m*τ for m in 0:M]; # Сетка по Т
ulₘ= u_l.(Tₘ); # Сеточные значения левого ГУ
urₘ= u_r.(Tₘ); # Сеточные значения правого ГУ
# Новое поведение будет реализоваться с помощью передачи функции создания сетки внутрь `solve` в качестве аргумента.
# Эта функция должна принимать только один аргумент — положение переходного слоя и формировать соответствующую сетку.
# В пакете есть формирование кусочно-равномерной сетки со сгущениями на границе и на переходном слое.
# Создадим замыкание этой функции, которое будет иметь нужную сигнатуру.
mshfrm(x_tp) = shishkin_mesh(a, b, x_tp, ε, 40, 1.0, 1.0, 0.2, 1.0);
Xₙ = mshfrm(x_tp); ; # Сетка по Х
N = length(Xₙ) - 1 # Примем за N длину сетки, что получилась в итоге.
qₙ = 1.0 * qf.(Xₙ); # Сеточные значения коэффициента линейного усиления
# чтобы посмотреть на поведение сопряженной задачи
# когда прямая задача не соответствует наблюдаемым данным,
# измените этот массив. Например qₙ = 0.2 * qf.(Xₙ);
u₀ = u_init.(Xₙ, ε=ε, x_tp = x_tp); # Начальные условия
#
u, XX, TP = solve(u₀, Xₙ, N, Tₘ, M, ε, ulₘ, urₘ, qₙ, create_mesh = mshfrm);
# Генерация априорной информации
ϕl = phidetermination(qₙ, ulₘ, Xₙ, N, Tₘ, M); # Левый вырожденный корень
ϕr = phidetermination(qₙ, urₘ,Xₙ, N, Tₘ, M, reverseX = true); # Нужно подать инвертированную сетку и q
ϕr = reverse(ϕr, dims=1); # А после — инвертировать решение по X
ϕ = Φ(ϕl, ϕr, N, M); # Полуразность вырожденных корней
ϕ = apply_on_dynamic_mesh(ϕ, XX, N, M); # Аппроксимация на переменную сетку
ϕl = apply_on_dynamic_mesh(ϕl, XX, N, M); # Аппроксимация на переменную сетку
ϕr = apply_on_dynamic_mesh(ϕr, XX, N, M); # Аппроксимация на переменную сетку
f1_data = f1(ϕ, u, XX, N, M); # Положение переходного слоя
f2_data = f2(f1_data, u, XX, N, M); # Значение функции на переходном слоеПро эмпирический параметр и его выбор смотри следующий подраздел.
w = 0.00015; # Эмпирический параметр апроксимации дельта-функции
# Решение сопряженной задачи
Uₙₘ = u; # Сохраним старую матрицу
ψ₀ = [0.0 for i in 1:N+1]; # Нулевые начальные условия
ψl = [0.0 for i in 1:M+1]; # Нулевые ГУ
ψr = [0.0 for i in 1:M+1]; # Нулевые ГУ
ψ = solve_adjoint(ψ₀, XX, N, Tₘ, M, ε, ψl, ψr, qₙ, Uₙₘ, f1_data, f2_data, w = w)Результат должен быть около нулевой, ведь в качестве текущего приближения q мы взяли искомое, а при нем — градиент должен обнуляться.
Контроль параметра аппроксимации дельта-функции
Проведем контроль выбора априорного параметра в аппроксимации дельта-функции. На точных данных, эта неоднородность должна быть нулевой почти везде исходя из следующих соображений:
- сопряженная задача — ретроспективная,
- её начальные условие — нулевые,
- если из уравнения вычеркнуть эту неоднородность, то решением с такими ГУ и НУ будет тривиальное нулевое решение.
- исходя из формулировки градиента,
- нулевое решение сопряженной задачи даст нам нулевой градиент, а значит мы нашли решение обратной задачи.
Выбор этого параметра значительно влияет на вычисление градиента, а значит и на решение, к которому сойдется итерационный процесс.
using NonLinearReactionAdvectionDiffusionWithFrontData: heterogeneity;
using LaTeXStrings, Plots;
Uₙₘ = u[2:N, :];
X = XX[2:N, :];
H = [ - heterogeneity(n, m, X[:, m], N, Tₘ, M, Uₙₘ, f1_data, f2_data, w) for n in 1:N-1, m in 1:M+1]
heatmap(H', title=L"-2 δ( x - f_1(t))(u(f1(t), t) - f2(t) ")Визуализация
ftw = isTravis ? range(1, stop = M+1, length=9) : [1; 2:div(M+1, 80):M+1];
make_gif(reverse(ψ, dims=2), XX, Tₘ;
label="\\psi",
name="solution_adjoint_ex3.gif",
frames_to_write = ftw)This page was generated using Literate.jl.