实现概述
洛伦兹吸引子是一组微分方程，我们可以用 c++ boost odeint 库来解微分方程。
2024 年 7 月 13 日
首页： 发一格 @cpp
数学上的描述
洛伦兹吸引子微分方程组
\begin{align}
&
\dot{x} = \sigma (y-x) &\\&
\dot{y} = x (\rho - z) -y &\\&
\dot{z} = xy - \beta z
&
\end{align}
它们都是时间 t 的导数
\begin{align}
&
\sigma, \rho, \beta \text{常见的取值为：} &\\&
\begin{bmatrix}
\sigma \\
\rho \\
\beta
\end{bmatrix}
=
\begin{bmatrix}
10 \\
25 \\
\frac{8}{3}
\end{bmatrix}
&
\end{align}
注：把一个点写在字母的正上方是导数在动力系统中的写法。
格式
要把表达式转化成可以对照着写 c++ 代码 (boost odeint) 的格式，减少精神内耗：
\begin{align}
&
\sigma = 10 \quad
\rho = 25 \quad
\beta = \frac{8}{3} &\\&
\begin{bmatrix}
x \\
y \\
z
\end{bmatrix}
\longrightarrow
\begin{bmatrix}
\dot{x} \\
\dot{y} \\
\dot{z}
\end{bmatrix}
=
\begin{bmatrix}
\sigma (y-x) \\
x (\rho - z) - y \\
xy - \beta z
\end{bmatrix}
&
\end{align}
c++ 代码
描述
编写 boost odeint 系统函数用于调微分方程组:
class lorenz_class
{
public:
void operator()(const state_type & in, state_type & diff, const double t) const
{
}
};
\begin{align}
&
\text{这里的 in 可以存储三个函数相关的值：} &\\&
x, y, z &\\&
\text{diff 可以存储对应的微分三个值：} &\\&
\dot{x}, \dot{y}, \dot{z} &\\&
&\\&
\dot{x}, \dot{y}, \dot{z}
\text{都是时间 t 的导数}
&
\end{align}
注意变量命名引起的混淆，如果把 in 写成 x, 那么可以理解它存储的三个数学值为 x1, x2, x3, 跟 x, y, z 对应； diff 也可以写成 dxdt, 命名不会改变计算规则。
微分方程
diff[0] = sigma * (in[1] - in[0]); // x' = sigma * (y-x) diff[1] = in[0] * (rho - in[2]) - in[1]; // y' = x * (rho - z) - y diff[2] = in[0] * in[1] - beta * in[2]; // z' = xy - beta * z
初始状态值的选择
我们可以选取空间中的某个点作为初始状态值，即可生成洛伦兹吸引子蝴蝶，但是，要避免选到吸引点，这样的点会导致无论迭代多少次，都不会逃逸吸引点，通常原点极有可能成为吸引点，所以不要选择原点作为初始状态值。完整的 c++ 代码
#include <boost/numeric/odeint.hpp>
#include <iostream>
#include <vector>
namespace odeint = boost::numeric::odeint;
using value_type = double;
using time_type = value_type;
using state_type = std::vector<value_type>;
constexpr value_type sigma = 10;
const value_type rho = 25;
const value_type beta = 8.0/3;
// 洛伦兹类
class lorenz_class
{
public:
// 微分方程
void operator()(const state_type & in, state_type & diff, const time_type t) const
{
diff[0] = sigma * (in[1] - in[0]); // x' = sigma * (y-x)
diff[1] = in[0] * (rho - in[2]) - in[1]; // y' = x * (rho - z) - y
diff[2] = in[0] * in[1] - beta * in[2]; // z' = x*y - beta * z
}
// 观察函数
void operator()(const state_type & result, const time_type t) const
{
std::cout << std::setw(10) << result[0] << std::setw(10) << result[1] << std::setw(10) << result[2]
<< std::endl;
}
};
int main()
{
auto lorenz = lorenz_class{}; // 洛伦兹对象
auto xyz = state_type{1, 1, 1}; // 初始状态
int steps = odeint::integrate_const(
odeint::runge_kutta4<state_type, value_type>{},
lorenz, // 微分方程
xyz,
time_type{0}, // 起始时间
time_type{10}, // 结束时间
time_type{0.1}, // dt 间隔
lorenz // 观察函数
);
std::cout << "steps: " << steps << std::endl;
}
编译运行输出：
$ g++ lorentz.cpp -std=c++23 -o lorentz
$ ./lorentz
1 1 1
2.06572 3.83607 1.05343
5.56918 11.1022 3.23336
13.817 23.2773 18.203
14.9366 8.52487 40.1873
5.47706 -4.22664 31.4177
-0.970417 -4.59924 23.4589
-3.2924 -5.22568 18.98
-5.28421 -7.76331 17.0037
-8.2712 -11.7431 18.9597
-11.2437 -13.0502 25.8689
-10.5808 -7.76616 30.2315
-7.04623 -3.54172 27.227
-4.58147 -3.25655 22.4171
-4.25916 -4.79911 18.6409
-5.61984 -7.76344 16.9502
-8.47063 -11.8382 19.1635
-11.3022 -12.902 26.1548
-10.4551 -7.52194 30.1788
-6.91835 -3.50379 27.0164
-4.55712 -3.32719 22.2464
-4.31587 -4.92683 18.5556
-5.74402 -7.95299 17.0076
-8.63419 -11.9909 19.4579
-11.3435 -12.7216 26.4987
-10.2971 -7.24131 30.1095
-6.76799 -3.46302 26.7702
-4.52877 -3.40915 22.0471
-4.38181 -5.07592 18.4565
-5.88888 -8.17353 17.0782
-8.82133 -12.1578 19.8055
-11.3798 -12.4973 26.8791
-10.1112 -6.92899 30.0135
-6.59986 -3.42281 26.488
-4.49941 -3.50373 21.8206
-4.45896 -5.24804 18.3461
-6.05613 -8.42591 17.1675
-9.03088 -12.3322 20.2117
-11.4049 -12.2215 27.2867
-9.89587 -6.59101 29.8833
-6.41701 -3.38733 26.1701
-4.47172 -3.61278 21.5692
-4.54949 -5.44522 18.2283
-6.24756 -8.71064 17.2823
-9.26079 -12.5056 20.682
-11.4117 -11.8869 27.7094
-9.64971 -6.2349 29.7114
-6.22313 -3.36101 25.8173
-4.44878 -3.73844 21.2956
-4.65592 -5.66984 18.1078
-6.46507 -9.02778 17.4309
-9.50793 -12.6669 21.221
-11.3923 -11.4876 28.1324
-9.37211 -5.86939 29.4907
-6.02249 -3.34834 25.4308
-4.43398 -3.88327 21.0031
-4.7812 -5.92476 17.9901
-6.71074 -9.37665 17.6231
-9.76781 -12.8022 21.8324
-11.3384 -11.0196 28.5379
-9.06323 -5.50416 29.2151
-5.81972 -3.35373 25.0121
-4.43098 -4.05029 20.6951
-4.92875 -6.21334 17.8821
-6.98668 -9.75555 17.8708
-10.0342 -12.8944 22.518
-11.2418 -10.4817 28.9061
-8.72426 -5.14935 28.8796
-5.61975 -3.38134 24.5629
-4.44365 -4.24302 20.3756
-5.10246 -6.53945 17.7919
-7.29484 -10.1611 18.1876
-10.2988 -12.9237 23.2757
-11.0942 -9.87663 29.2152
-8.35777 -4.8153 28.4808
-5.42772 -3.43499 24.085
-4.47605 -4.46543 20.0489
-5.30658 -6.90717 17.7291
-7.63651 -10.587 18.5885
-10.5503 -12.8671 24.0978
-10.8887 -9.21232 29.4424
-7.96822 -4.51223 28.0171
-5.24898 -3.51814 23.5807
-4.53237 -4.72181 19.7204
-5.54547 -7.32016 17.7052
-8.01142 -11.0226 19.0893
-10.7744 -12.701 24.9684
-10.62 -8.50304 29.5658
-7.56248 -4.25004 27.4895
-5.08924 -3.63386 23.0539
-4.61691 -5.01653 19.3966
-5.82326 -7.78079 17.7339
-8.41663 -11.4513 19.7046
-10.9535 -12.403 25.8604
-10.2857 -7.76991 29.5659
-7.15008 -4.03791 26.9024
-4.95463 -3.78505 22.5102
-4.73417 -5.35387 19.0867
-6.14344 -8.28893 17.8317
-8.8452 -11.8485 20.4454
-11.068 -11.9572 26.7351
steps: 100
输出数据
输出数据可以完整的发给 c++ sfml 图形 api 渲染成 3D 图像，可以逐步发送给 c++ sfml api，实时渲染，做成动画，绘图这部分暂不放出，下次待续。