以声波为例的偏微分方程数值解

在这个例子中，我们将使用声波振荡的例子来考虑偏微分方程的数值解。
嵌入式函数的主题将被触及，以展示如何在不降低代码执行速度的情况下提高代码的可读性。

让我们考虑字符串振荡方程的解，而不深入代码，而只关注其结构。
（代码将在下一章中解释。)
``'茱莉亚
@views函数acoustic_wave_1D_orig()
#物理参数
lx=
20.0r，β=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)
样本=复制(Pr)
Vx=零(Float64,nx-1)
#算法周期
@gif为它=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="浓度",
标题="时间= $(round(itdt,digits=1))")
end every nvis
end


Несмотря на то, что код сопровождается комментариями, можно явно выделить этапы работы функции.
Может возникнуть предположение, что внедрение функций в код под каждый этап замедлит выполнение главной функции.

Действительно, вызов функции [занимает значительное время](https://habr.com/ru/companies/otus/articles/343566/), в сравнении с арифметическими операциями. Однако компиляторы умеют оптимизировать код и делать функции встраиваемыми. 

Чтобы проверить, происходит ли вызов функции - воспользуемся макросом `@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)

;  @ In[2]:5 within `fb`
define i64 @julia_fb_11486(i64 signext %0, i64 signext %1, i64 signext %2) #0 {
top:
;  @ In[2]:6 within `fb`
; ┌ @ In[2]:2 within `fa`
; │┌ @ int.jl:87 within `+`
    %3 = add i64 %1, %0
; └└
; ┌ @ int.jl:88 within `*`
   %4 = mul i64 %3, %2
   ret i64 %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)

;  @ In[3]:5 within `fb`
define i64 @julia_fb_11528(i64 signext %0, i64 signext %1, i64 signext %2) #0 {
top:
;  @ In[3]:6 within `fb`
  %3 = call i64 @j_fa_11530(i64 signext %0, i64 signext %1)
; ┌ @ int.jl:88 within `*`
   %4 = mul i64 %3, %2
   ret i64 %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

test_fa_inline (generic function with 1 method)

И подключим библиотеку BenchmarkTools для оценки времени выполнения кода:

Pkg.add("BenchmarkTools")

using BenchmarkTools
N = 10^5
@btime test_fa_noinline($N)

  157.812 μs (0 allocations: 0 bytes)

10000100000

@btime test_fa_inline($N)

  2.875 ns (0 allocations: 0 bytes)

10000100000

Как видим, встроенная функция работает, быстрее.

Уравнение акустической волны

Уравнение акустической волны имеет вид:

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

Введём некоторые обозначения:

- объёмный модуль сжимаемости (βulk modulus)

- плотность среды.
Тогда исходное уравнение второго порядка:

можно переписать как систему двух уравнений первого порядка:

声波方程的数值解

让我们介绍一个定义方程物理参数的函数:

@doc """
- lx - длина среды, по которой распространяется волна
- ρ  - плотность
- β  - модуль сжимаемости
"""
physical_parameters(;lx = 20, ρ = 1, β = 1) = (lx,ρ,β)
@doc physical_parameters

我们以这样的方式设置它，以演示以下场景:

physical_parameters() # хотим использовать параметры по умолчанию

(20, 1, 1)

physical_parameters(lx=100,β=10) # меняем какие-то параметры, остальные оставляем по умолчанию

(100, 1, 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)

;  @ In[11]:1 within `foo`
define i64 @julia_foo_12020(i64 signext %0, i64 signext %1, i64 signext %2) #0 {
top:
;  @ In[11]:4 within `foo`
; ┌ @ operators.jl:587 within `+` @ int.jl:87
   %3 = add i64 %1, %0
   %4 = add i64 %3, %2
; └
  ret i64 %4
}

如果只 physical_parameters 我没有归还车队 (lx, ρ, β)，和向量 [lx, ρ, β] 然后汇编程序代码会更大，需要更长的时间来执行。

有限元法

要解决方程组，我们将使用[有限差分]方法（https://en.wikipedia.org/wiki/Finite_difference_method ）与欧拉方法相结合。
（通常，要解决此类问题，[应用]（https://en.wikipedia.org/wiki/FTCS_scheme ）更准确[中心差]（https://ru.wikipedia.org/wiki/Finite 差异））

该方法的本质是用小但有限的差异来表示偏微分方程。:

选择采样步骤

具有"字符串"的长度 lx 我们可以通过将长度除以我们自己确定的分区数量来获得空间离散化。

nx -由用户设置
dx = lx / nx -空间采样步骤
-采样的时间步长。他的选择会影响解决方案是否可持续。此步骤的值为[selected](https://d-arora.github.io/Doing-Physics-With-Matlab/mpDocs/em_swe.pdf#page=5 ）基于[Courant条件]（https://en.wikipedia.org/wiki/Courant–Friedrichs–Lewy_condition )。在我们的例子中，介质中的声速充当常数。
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

differentials (generic function with 1 method)

初始化数组

让我们设置对应于压力和速度的数组。

我们还将保存初始压力数组以供后续可视化。

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

arrays_initializer (generic function with 1 method)

在继续讨论解决方案之前，我们将连接一个库，用于可视化解决方案的图形并对它们进行动画处理。:

# Pkg.add("Plots")
using Plots, Plots.Measures, Printf
# задаём параметры по умолчанию
default(size=(1500, 600), framestyle=:box, label=false, grid=false, margin=10mm, lw=6, labelfontsize=20, tickfontsize=20, titlefontsize=24)

有限差分解

a = a - 1 ~ a -= 1

让我们重写它

如何

-标量，显示给定时间的值在点

有关更多信息，请参阅下图。

同样，对于压力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

show_acoustic_wave (generic function with 1 method)

然后最终的代码将如下所示:

@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()

 28.463060 seconds (18.00 M allocations: 1.377 GiB, 1.71% gc time, 47.50% compilation time: 4% of which was recompilation)

让我们将执行时间与我们的源代码进行比较。

@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()

 14.281645 seconds (3.29 M allocations: 400.516 MiB, 0.73% gc time, 0.65% compilation time)

由此我们可以得出结论，在函数中包装小代码不会显着影响其执行速度。但是，可以注意将代码划分为函数的以下优点:

功能文件。记录函数是比注释代码更有效和更具表现力的工具。
您可能会在不同的场景中使用某些功能
创建函数允许您将复杂的任务分解为简单的任务，以及逐渐增加程序的复杂性和准确性。（首先，我们将插入一个简单的函数，但我们将检查其可操作性，只有这样我们才会选择更准确和有效的方法来实现它）