boost 中的 Runge Kutta k 向量

Question

使用 boost::odeint 求解 ODE 系统时，我可以使用 stepper.do_step 函数及时控制每次迭代。这样，我可以存储 u(t)。

我正在研究一种算法，该算法需要在每次迭代时存储 s 阶段 Runge-Kutta 方法的 k 向量。例如，对于 RK4，我需要存储 k1、k2、k3 和 k4 的值。

我是否可以从 boost::odeint 库中获取这些值？我认为没有内置的方法可以做到这一点。

我尝试在 generic_rk_algorithm.hpp 中搜索以某种方式更改库的方法。我发现在 do_step 函数中，我想要 calculate_stage 结构的 F 属性的取消引用值，但我不确定如何更改代码来获取它。

Answer 1

刚落地。 这里是 Chua 电路吸引子实现的代码。 如果我正确理解了这个问题，那么这个例子应该对你有帮助。 在此示例中，存在三个变量，其中每个变量的 k 值都存储在每个时间间隔。希望，这有帮助。

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# parameters
alpha = 9.0
beta = 14.286
a = 2./7.
b = -1./7.
c = 1

def chua(u):
    x, y, z = u  
    f_x = b*x + 0.5*(a-b)*(abs(x+c)-abs(x-c))
    dudt = [alpha*(y - f_x), x - y + z, -beta*y]
    return dudt

# time discretization
t_0 = 0
tau = 0.01
t_final = 600.0
t = np.linspace(t_0, t_final, int(t_final/tau))
k1 = np.zeros((int(t_final/tau),3), dtype='float')
k2 = np.zeros((int(t_final/tau),3), dtype='float')
k3 = np.zeros((int(t_final/tau),3), dtype='float')
k4 = np.zeros((int(t_final/tau),3), dtype='float')

# initial conditions
u0 = [0.1,0.0,0.0]
# integrate ode system
sol = np.zeros((int(t_final/tau),3))
sol[0] = u0
for i in range(1, int(t_final/tau)):
    k1[i-1] = tau * np.asarray(chua(sol[i-1]))
    k2[i-1] = tau * np.asarray(chua(sol[i-1] + 0.5 * k1[i-1]))
    k3[i-1] = tau * np.asarray(chua(sol[i-1] + 0.5 * k2[i-1]))
    k4[i-1] = tau * np.asarray(chua(sol[i-1] + k3[i-1]))
    sol[i] = sol[i-1] + (1.0 / 6.0)*(k1[i-1] + 2*(k2[i-1] + k3[i-1]) + k4[i-1])

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
ax.plot(sol[:,0],sol[:,1],sol[:,2])
plt.show()