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

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

我们将在不深入研究代码的情况下考虑弦振荡方程的解法，只关注其结构。 (代码将在下一章解释）。 ``朱莉娅 @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)

В строчке %3 = add i64 %1, %0 происходит сложение чисел a и b. Для того чтобы явно указать компилятору, что не нужно делать оптимизацию встраивания, воспользуемся макросом @noinline

Теперь же происходит вызов функции fa, которая складывает a и b:

%3 = call i64 @j_fa_7789(i64 signext %0, i64 signext %1)

Для того чтобы убедиться в том, что встроенная функция работает быстрее, чем вызываемая - напишем следующий код:

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

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

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

Уравнение акустической волны имеет вид: $\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.$$

声波方程的数值解法¶

让我们引入一个指定方程物理参数的函数：

我们将其设置为演示以下情况：

请注意，元组是一种特殊的数据类型，值得了解何时应用。

不过，对于值得解包的小数据集来说，元组是最佳选择之一。 在下面的函数foo 中，我们可以看到编译器认为 "解包 "一个元组就是简单地将a,b,c 赋值给x 、y 和z 。而且不需要对元组中的元素进行任何调用。

如果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 是进行空间离散化的网格中心集

数组初始化¶

让我们设置与压力和速度相对应的数组。 我们还将保存压力的初始数组，以便稍后进行可视化。

在继续求解之前，让我们连接求解图形可视化和动画库：

有限差分解法¶

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

这样，最终代码的形式就如下所示：

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

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