Численное решение уравнения в частных производных на примере акустической волны¶
В данном примере будет рассмотрено численное решение уравнения в частных производных на примере колебания акустической волны. Будет затронута тема встраиваемых функций, чтобы показать, как можно повысить читабельность кода, не уменьшая при этом скорость его выполнения.
Рассмотрим решение уравнения колебания струны, не углубляясь в код, а лишь обратив внимание на его структуру. (Пояснения к коду будут даны в следующей главе).
@views function acoustic_wave_1D_orig()
# физический параметры
lx = 20.0
ρ, β = 1.0, 1.0
# отсчёты
nx = 200
nvis = 2
# дифференциалы
dx = lx / nx
dt = dx*√(ρ / β)
nt = nx^2 ÷ 100
xc = LinRange(dx / 2, lx - dx / 2, nx)
# инициализация массивов
Pr = @. exp(-(xc - lx / 4)^2)
Sample = copy(Pr)
Vx = zeros(Float64, nx - 1)
# цикл работы алгоритма
@gif for it = 1:nt
Vx .-= dt ./ ρ .* diff(Pr) ./ dx
Pr[2:end-1] .-= dt ./ β .* diff(Vx) ./ dx
plot(xc, [Sample, Pr]; xlims=(0, lx), ylims=(-1.1, 1.1),
xlabel="lx", ylabel="Концентрация",
title="Время = $(round(it*dt,digits=1))")
end every nvis
end
Несмотря на то, что код сопровождается комментариями, можно явно выделить этапы работы функции. Может возникнуть предположение, что внедрение функций в код под каждый этап замедлит выполнение главной функции.
Действительно, вызов функции занимает значительное время, в сравнении с арифметическими операциями. Однако компиляторы умеют оптимизировать код и делать функции встраиваемыми.
Чтобы проверить, происходит ли вызов функции - воспользуемся макросом @code_llvm, который покажет код, близкий к ассемблеру. (Для большей читабельности можно включить тёмную тему в Engee)
function fa(a,b)
return a+b
end
function fb(a,b,c)
fa(a,b)*c
end
@code_llvm fb(1,2,4)
В строчке %3 = add i64 %1, %0 происходит сложение чисел a и b.
Для того чтобы явно указать компилятору, что не нужно делать оптимизацию встраивания, воспользуемся макросом @noinline
@noinline function fa(a,b)
return a+b
end
function fb(a,b,c)
fa(a,b)*c
end
@code_llvm fb(1,2,4)
Теперь же происходит вызов функции fa, которая складывает a и b:
%3 = call i64 @j_fa_7789(i64 signext %0, i64 signext %1)
Для того чтобы убедиться в том, что встроенная функция работает быстрее, чем вызываемая - напишем следующий код:
@noinline function fa_noinline(a,b)
a + b
end
@inline function fa_inline(a,b)
a + b
end
function test_fa_noinline(N)
sum = 0
for i in 1:N
sum += fa_noinline(i,i)
end
sum
end
function test_fa_inline(N)
sum = 0
for i in 1:N
sum += fa_inline(i,i)
end
sum
end
И подключим библиотеку BenchmarkTools для оценки времени выполнения кода:
import Pkg; Pkg.add("BenchmarkTools")
using BenchmarkTools
N = 10^5
@btime test_fa_noinline($N)
@btime test_fa_inline($N)
Как видим, встроенная функция работает, быстрее.
Уравнение акустической волны¶
Уравнение акустической волны имеет вид: $\frac{\partial^2 P}{\partial t^2} = c^2 \frac{\partial^2 P}{\partial x^2}$
где $P$ - давление, $t$ - время, а $x$ - положение частиц в некоторой среде. $c$ - скорость звука в этой среде.
Введём некоторые обозначения:
$c^2 = \frac{1}{\beta \rho}$
$\frac{1}{\beta}$ - объёмный модуль сжимаемости (βulk modulus)
$\rho$ - плотность среды. Тогда исходное уравнение второго порядка: $\frac{\partial^2 P}{\partial t^2} = c^2 \frac{\partial^2 P}{\partial x^2}$ можно переписать как систему двух уравнений первого порядка:
$$\left\{ \begin{array}{cl} \frac{\partial V_x }{\partial t} = -\frac{1}{\rho} \frac{\partial P}{\partial x} \\ \frac{\partial P }{\partial t} = -\frac{1}{\beta} \frac{\partial V_x}{\partial x} \end{array} \right.$$
Численное решение уравнения акустической волны¶
Введём функцию, задающую физические параметры нашего уравнения:
@doc """
- lx - длина среды, по которой распространяется волна
- ρ - плотность
- β - модуль сжимаемости
"""
physical_parameters(;lx = 20, ρ = 1, β = 1) = (lx,ρ,β)
@doc physical_parameters
Мы её задали таким образом, чтобы продемонстрировать следующие сценарии:
physical_parameters() # хотим использовать параметры по умолчанию
physical_parameters(lx=100,β=10) # меняем какие-то параметры, остальные оставляем по умолчанию
Отметим, что кортеж - специфический тип данных, который стоит понимать, когда применять.
Однако для небольшого набора данных, которые стоит распаковать, кортеж - один из лучших вариантов.
В представленной ниже функции foo, мы можем заметить, что компилятор увидел, что "распаковка" кортежа - это просто присвоение a,b,c значений x,y и z. И не нужно делать никаких обращений к элементам кортежа.
function foo(x,y,z)
a,b,c = physical_parameters(;lx = x, ρ=y, β=z)
return a+b+c
end
@code_llvm foo(10,20,30)
Если бы physical_parameters возвращала не кортеж (lx, ρ, β), а вектор [lx, ρ, β], то ассемблерный код был бы намного больше и выполнялся дольше.
Метод конечных элементов¶
Для решения системы уравнений будем применять метод конечных разностей совместно с методом Эйлера. (Часто для решения такого рода задач применяют более точную центральную разность)
Суть метода состоит в том, чтобы представить уравнение в частных производных, через малые, но конечные разности:
$\frac{\partial V }{\partial t} = \frac{\partial P}{\partial x} \sim \frac{V_{j}^{n+1} - V_{j}^{n}}{ \Delta t } =\frac{P_{j}^{n+1} - P_{j}^{n}}{\Delta x}$
Выбор шага дискретизации¶
Имея длину "струны" lx, мы можем получить пространственную дискретизацию, разделив длину на количество разбиений, которое мы сами определяем.
nx- задаётся пользователемdx=lx / nx- пространственный шаг дискретизации- $dt = dx/\sqrt{\frac{1}{\beta \rho}}$ - временной шаг дискретизации. Его выбор влияет на то, будет ли решение устойчивым. Значения данного шага выбирается исходя из условия Куранта. В нашем случае в качестве константы выступает скорость звука в среде.
ntможет задаваться пользователем, а может зависеть от длины, чтобы всегда проходить хотя бы 1 периодxc- набор центров сетки, по которой проводится пространственная дискретизация
function differentials(;lx=20,nx=200,β=1,ρ=1)
dx=lx/nx
dt=dx*√(β*ρ)
nt = nx^2 ÷ 100
xc = LinRange(dx / 2, lx - dx / 2, nx)
return(dx=dx,dt=dt,nt=nt,xc=xc)
end
инициализация массивов¶
Зададим массивы, соответствующие давлению и скорости. Также для последующей визуализации сохраним исходный массив давлений.
function arrays_initializer(xc,lx,nx)
Vx = zeros(Float64, nx - 1)
Pr = @. exp(-(xc - lx / 4)^2)
Sample = copy(Pr)
(Vx=Vx, Pr=Pr, Sample=Sample)
end
Перед тем, как перейти к решению, подключим библиотеку для визуализации графиков решений и их анимации:
import Pkg; Pkg.add("Plots")
using Plots, Plots.Measures, Printf
# задаём параметры по умолчанию
default(size=(1200, 400), framestyle=:box, label=false, grid=false, margin=10mm, lw=6, labelfontsize=20, tickfontsize=20, titlefontsize=24)
Решение методом конечных разностей¶
a = a - 1 ~ a -= 1
Перепишем $\frac{\partial V }{\partial t} = \frac{\partial P}{\partial x} \sim \frac{V_{j}^{n+1} - V_{j}^{n}}{ \Delta t } =\frac{P_{j}^{n+1} - P_{j}^{n}}{\Delta x}$
как
$V_x^{n+1} = V_x^n - \frac{1}{\rho} \Delta t (P^{n+1}-P^n)/ \Delta x$
$V_{j}^{n}$ - скаляр, показывающий значение в момент времени $n$ в точке $j$
Подробнее - смотрите рисунок ниже.
Аналогично и для давления P.
function show_acoustic_wave((Vx,Pr,Sample),(dx,dt,nt,xc);nvis=2,ρ=1,β=1,lx=20)
@gif for it = 1:nt
Vx .-= dt ./ ρ .* diff(Pr) ./ dx
Pr[2:end-1] .-= dt ./ β .* diff(Vx) ./ dx
plot(xc, [Sample, Pr]; xlims=(0, lx), ylims=(-1.1, 1.1),
xlabel="lx", ylabel="Концентрация",
title="Время = $(round(it*dt,digits=1))")
end every nvis
end
Тогда финальный код будет иметь следующее вид:
@views function acoustic_wave_1D()
lx, ρ, b = physical_parameters() # если мы знаем порядок элементов в кортеже
nx = 200 # количество ячеек в сетке
ds = differentials() # если не знаем, то используем именованный
dx = ds.dx
dt = ds.dt
nt = ds.nt
xc = ds.xc
ars = arrays_initializer(xc,lx,nx)
show_acoustic_wave(ars,ds,nvis=2)
end
@time acoustic_wave_1D()
Сравним время выполнения с нашим исходным кодом.
@views function acoustic_wave_1D_orig()
# физический параметры
lx = 20.0
ρ, β = 1.0, 1.0
# отсчёты
nx = 200
nvis = 2
# дифференциалы
dx = lx / nx
dt = dx*√(ρ / β)
nt = nx^2 ÷ 100
xc = LinRange(dx / 2, lx - dx / 2, nx)
# инициализация массивов
Pr = @. exp(-(xc - lx / 4)^2)
Sample = copy(Pr)
Vx = zeros(Float64, nx - 1)
# цикл работы алгоритма
@gif for it = 1:nt
Vx .-= dt ./ ρ .* diff(Pr) ./ dx
Pr[2:end-1] .-= dt ./ β .* diff(Vx) ./ dx
plot(xc, [Sample, Pr]; xlims=(0, lx), ylims=(-1.1, 1.1),
xlabel="lx", ylabel="Концентрация",
title="Время = $(round(it*dt,digits=1))")
end every nvis
end
@time acoustic_wave_1D_orig()
Из чего мы можем сделать вывод, что оборачивание небольшого кода в функции не влияет существенно на скорость её выполнения. Однако можно отметить следующие преимущества разделения кода на функции:
- Документирование функций. Документирование функций - намного более эффективный и выразительный инструмент, чем комментирование кода.
- Некоторые функции могут применяться вами в разных сценариях
- Создание функций позволяет декомпозировать сложные задачи на простые, а также постепенно повышать сложность и точность вашей программы. (Сначала вставим простую функцию, но проверим работоспособность, а уже потом будет подбирать более точные и эффективные способы её реализации)