以声波为例的偏导数方程数值解法

在本例中，我们将以声波振荡为例，研究偏导数方程的数值解法。
其中将涉及内联函数的主题，以说明如何在不降低执行速度的情况下提高代码的可读性。

我们将在不深入研究代码的情况下考虑弦振荡方程的解法，只关注其结构。
(代码将在下一章解释）。
``朱莉娅
@views 函数 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="Concentration"、
title="时间 $(round(it*dt,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[1]:5 within `fb`
define i64 @julia_fb_7009(i64 signext %0, i64 signext %1, i64 signext %2) #0 {
top:
;  @ In[1]:6 within `fb`
; ┌ @ In[1]: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[2]:5 within `fb`
define i64 @julia_fb_7058(i64 signext %0, i64 signext %1, i64 signext %2) #0 {
top:
;  @ In[2]:6 within `fb`
  %3 = call i64 @j_fa_7060(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 для оценки времени выполнения кода:

# import Pkg; Pkg.add("BenchmarkTools")
using BenchmarkTools
N = 10^5
@btime test_fa_noinline($N)

  156.660 μs (0 allocations: 0 bytes)

10000100000

@btime test_fa_inline($N)

  3.161 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[9]:1 within `foo`
define i64 @julia_foo_7547(i64 signext %0, i64 signext %1, i64 signext %2) #0 {
top:
;  @ In[9]: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, ρ, β) ，汇编代码就会变得更大，执行时间也会更长。

有限元方法

为了求解方程组，我们将结合欧拉方法使用有限差分方法。
(通常使用更精确的中心差分来解决这类问题）。

该方法的本质是通过微小但有限的差分，用偏导数表示方程：

选择采样步骤

有了 "字符串 "的长度lx ，我们就可以通过用长度除以我们自己定义的分区数来获得空间离散度。

-nx - 用户定义
-dx =lx / nx - 空间离散化步骤
- - 时间离散化步骤。它的选择会影响解是否稳定。这一步的值是根据库朗条件选择的。在本例中，常数为介质中的声速。
-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()

 12.091866 seconds (2.74 M allocations: 334.784 MiB, 0.59% gc time)

[ Info: Saved animation to /user/start/examples/math_and_optimization/pde_acoustic_wave/tmp.gif

让我们比较一下源代码的执行时间。

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

 13.480746 seconds (2.74 M allocations: 334.462 MiB, 0.89% gc time)

[ Info: Saved animation to /user/start/examples/math_and_optimization/pde_acoustic_wave/tmp.gif

由此我们可以得出结论，将小代码封装在函数中并不会明显影响其执行速度。不过，我们可以注意到将代码拆分为函数的以下优点：
1.1. 函数的文档化。与注释代码相比，记录函数是一种更有效、更有表现力的工具。
2.某些函数可以在不同的场景中使用
3.创建函数可以将复杂任务分解为简单任务*，逐步提高程序的复杂性和准确性。(让我们先插入一个简单的函数，但要检查它是否有效，然后我们再寻找更精确、更有效的方法来实现它）