8.1. Метод комплексного шага#
Метод комплексного шага (complex step method) [KW19, MSA03] относится к численным методам дифференцирования (наряду с методом конечных разностей). Однако, арифметика, стоящая за методом, роднит его с методом автоматического дифференцирования вперёд.
Заключается метод в следующем. Вычислим функцию \(f\) один раз в точке \(x + i h\) (\(i h\) и есть комплексный шаг). Разложение Тейлора для \(f(x + ih)\) имеет вид
Мнимая часть разложения
позволяет вычислить производную \(f'(x)\) с точностью \(O(h^2)\)
В свою очередь, действительная часть разложения
содержит значение \(f(x)\) с той же точностью, \(O(h^2)\)
Метод комплексного шага прост в реализации.
Метод комплексного шага
"""
diff_complexstep(f, x; h=1e-16)
Производная скалярной функции `f` в точке `x` методом комплексного шага,
размер шага `h`.
"""
diff_complexstep(f, x; h=1e-16) = imag(f(x + im * h)) / h
Как видно из реализации, метод можно применять к любой функции, вычислимой от комплекcного аргумента. При этом от программиста требуется определить только дифференцируемую функцию (скажем, уравнение состояния), а её производные можно вычислить, выполнив вычисление от комплексного числа.
В демонстрации ниже показан пример вычисления производной функции методом комплексного шага, а также сравнивается точность этого метода с парой методов конечных разностей.
Рассмотрим три численных метода дифференцирования: методы конечных разностей (прямой и центральный вариант) и метод комплексного шага.
diff_complexstep(f, x, h=1e-16) = imag(f(x + im * h)) / h
diff_forward(f, x, h=1e-8) = (f(x + h) - f(x)) / h
diff_central(f, x, h=1e-6) = (f(x + h/2) - f(x - h/2)) / h
Сравним точности методов на примере вычисления производной синуса в точке \(x = 1/3\).
f(x) = sin(x)
df(x) = cos(x)
x = 1/3
h = [2.0^i for i in -52:0]
err_complex = @. abs(1 - diff_complexstep(f, x, h)/df(x))
err_forward = @. abs(1 - diff_forward(f, x, h)/df(x))
err_central = @. abs(1 - diff_central(f, x, h)/df(x))
График зависимости невязки от размера шага \(h\) выглядит следующим образом (\(\log-\log\) шкала).
using Plots
plot(;
xaxis=:log10,
yaxis=:log10,
xlabel="Размер шага h",
ylabel="Норма относительной невязки",
leg_title="Метод", leg=:top
)
plot!(h[err_complex .> 0], err_complex[err_complex .> 0]; label="комплексный шаг")
plot!(h, err_forward; label="разность вперёд")
plot!(h, err_central; label="центральная разность")
xticks!([10.0^i for i in -14:2:0])
При этом значения невязки для метода комплексного шага при шагах менее \(1.5 \times 10^{-8}\) получились нулевыми (поэтому они не насенены на график в логарифмической шкале).
[h err_complex][err_complex .== 0, :]
27×2 Matrix{Float64}:
2.22045e-16 0.0
4.44089e-16 0.0
8.88178e-16 0.0
1.77636e-15 0.0
3.55271e-15 0.0
7.10543e-15 0.0
1.42109e-14 0.0
2.84217e-14 0.0
5.68434e-14 0.0
1.13687e-13 0.0
2.27374e-13 0.0
4.54747e-13 0.0
9.09495e-13 0.0
⋮
7.27596e-12 0.0
1.45519e-11 0.0
2.91038e-11 0.0
5.82077e-11 0.0
1.16415e-10 0.0
2.32831e-10 0.0
4.65661e-10 0.0
9.31323e-10 0.0
1.86265e-9 0.0
3.72529e-9 0.0
7.45058e-9 0.0
1.49012e-8 0.0
Так, метод комплексного шага более устойчив, чем методы конечных разностей, к ошибкам округления, возникающих в арифметике с конечной точностью. Устойчивость кроется в избегании эффекта subtractive cancellation за счёт вычисления функции один раз (вместо двух).
Метод комплескного шага гораздо лучше устойчив к численным эффектам, чем методы конечных разностей. Тем не менее, его точность аналитически ограничена (8.1). Далее мы увидим, что можно получить точное значение производной, если вычислить функцию не от комплексного, а от дуального числа.