Arenstorf orbit =============== In this example we would like to present some adaptive time step methods in Ponio. We do that on Arenstorf orbit problem .. math:: \begin{cases} \ddot{x} = x + 2\dot{y} - \frac{1-\mu}{r_1^3}(x + \mu) - \frac{\mu}{r_2^3}(x-1+\mu) \\ \ddot{y} = y - 2\dot{x} - \frac{1-\mu}{r_1^3}y - \frac{\mu}{r_2^3}y \end{cases} with initial conditon :math:`(x,\dot{x},y,\dot{y})=(0.994, 0, 0, -2.001585106)` and :math:`r_1` and :math:`r_2` done by .. math:: r_1 = \sqrt{(x+\mu)^2 + y^2},\qquad r_2 = \sqrt{(x-1+\mu)^2 + y^2}. .. code:: ipython3 import sympy as sp The ponio library solves only first order derivative in time, to do that we need to rewrite our problem with the folowing vector .. math:: \begin{pmatrix} y_1 \\ y_2 \\ y_3 \\ y_4 \end{pmatrix} = \begin{pmatrix} x \\ y \\ \dot{x} \\ \dot{y} \end{pmatrix} the system becomes .. math:: \begin{cases} \dot{y}_1 = y_3 \\ \dot{y}_2 = y_4 \\ \dot{y}_3 = y_1 + 2y_4 - \frac{1-\mu}{r_1^3}(y_1 + \mu) - \frac{\mu}{r_2^3}(y_1-1+\mu) \\ \dot{y}_4 = y_2 - 2y_3 - \frac{1-\mu}{r_1^3}y_2 - \frac{\mu}{r_2^3}y_2 \\ \end{cases} We already present the Lawson splitting into a linear and a non-linear part: .. math:: \frac{\mathrm{d}}{\mathrm{d}t}\begin{pmatrix} y_1 \\ y_2 \\ y_3 \\ y_4 \end{pmatrix} = \underbrace{ \begin{pmatrix} 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ 1 & 0 & 0 & 2 \\ 0 & 1 &-2 & 0 \\ \end{pmatrix} }_{L} \begin{pmatrix} y_1 \\ y_2 \\ y_3 \\ y_4 \end{pmatrix} + \underbrace{ \begin{pmatrix} 0 \\ 0 \\ - \frac{1-\mu}{r_1^3}(y_1 + \mu) - \frac{\mu}{r_2^3}(y_1-1+\mu) \\ - \frac{1-\mu}{r_1^3}y_2 - \frac{\mu}{r_2^3}y_2 \end{pmatrix} }_{N(t,u)} For stability condition of Lawson methods we check the eigen values of the linear part. .. code:: ipython3 L = sp.Matrix([ [0,0,1,0], [0,0,0,1], [1,0,0,2], [0,1,-2,0], ]) t = sp.symbols(r"tau",real=True) .. code:: ipython3 for ev in L.eigenvals(): display(ev) .. math:: \displaystyle - i .. math:: \displaystyle i .. math:: \textrm{sp}(L) \subset i\mathbb{R} \implies \textrm{sp}(e^{\tau L})\subset \mathcal{C}(0,1), \forall \tau\in\mathbb{R} so the linear part of the Lawson splitting doesn’t blow up, this is a good property for the stability of the method. To be solved easily with a Lawson method we need to compute the exponential of the linear part : :math:`e^{\tau L}`, with differents values of :math:`\tau = c_i\Delta t`. .. code:: ipython3 eL = sp.exp(t*L) display(eL) .. math:: \displaystyle \left[\begin{matrix}\tau \sin{\left(\tau \right)} + \cos{\left(\tau \right)} & - \tau \cos{\left(\tau \right)} + \sin{\left(\tau \right)} & \tau \cos{\left(\tau \right)} & \tau \sin{\left(\tau \right)}\\\tau \cos{\left(\tau \right)} - \sin{\left(\tau \right)} & \tau \sin{\left(\tau \right)} + \cos{\left(\tau \right)} & - \tau \sin{\left(\tau \right)} & \tau \cos{\left(\tau \right)}\\\tau \cos{\left(\tau \right)} & \tau \sin{\left(\tau \right)} & - \tau \sin{\left(\tau \right)} + \cos{\left(\tau \right)} & \tau \cos{\left(\tau \right)} + \sin{\left(\tau \right)}\\- \tau \sin{\left(\tau \right)} & \tau \cos{\left(\tau \right)} & - \tau \cos{\left(\tau \right)} - \sin{\left(\tau \right)} & - \tau \sin{\left(\tau \right)} + \cos{\left(\tau \right)}\end{matrix}\right] We get an explicit form of :math:`e^{\tau L}`, :math:`\forall\tau\in\mathbb{R}`, so we don’t need to compute numerically exponential at each step. Now we can write a special exponential function to return this matrix at each step. In ponio, exponential function for Lawson method gets the linear part times the time step: .. code:: cppr matrix_t exp( matrix_t const& M ); with ``M = c_j*dt*L``. To get the time step of the stage, with this matrix we can get ``M(0,2)``. .. code:: cpp Eigen::Matrix exp ( Eigen::Matrix const& M ) { double tau = M(0,2); double s = std::sin(tau), c = std::cos(tau); Eigen::Matrix e; e << tau*s+c , -tau*c+s , tau*c , tau*s , tau*c-s , tau*s+c , -tau*s , tau*c , tau*c , tau*s , -tau*s+c , tau*c+s , -tau*s , tau*c , -tau*c-s , -tau*s+c ; return e; } We can also use SymPy to generate the code .. code:: ipython3 print(""" e << {}; """.format( " ,\n ".join([ " , ".join([ sp.cxxcode(eL[i,j]) for j in range(eL.cols) ]) for i in range(eL.rows) ]) )) .. parsed-literal:: e << tau*std::sin(tau) + std::cos(tau) , -tau*std::cos(tau) + std::sin(tau) , tau*std::cos(tau) , tau*std::sin(tau) , tau*std::cos(tau) - std::sin(tau) , tau*std::sin(tau) + std::cos(tau) , -tau*std::sin(tau) , tau*std::cos(tau) , tau*std::cos(tau) , tau*std::sin(tau) , -tau*std::sin(tau) + std::cos(tau) , tau*std::cos(tau) + std::sin(tau) , -tau*std::sin(tau) , tau*std::cos(tau) , -tau*std::cos(tau) - std::sin(tau) , -tau*std::sin(tau) + std::cos(tau); Or for the shake of simplicity, use Eigen functions to compute exponential of a matrix: .. code:: cpp #include // ... auto mexp = [](Eigen::Matrix const& M){ return M.exp(); }; .. code:: ipython3 %system mkdir -p arenstorf_adaptivetimestep_demo .. parsed-literal:: [] .. code:: ipython3 %%writefile arenstorf_adaptivetimestep_demo/main.cpp #include #include #include #include #include #include #include #include struct arenstorf_model { using state_t = Eigen::Matrix; using matrix_t = Eigen::Matrix; double mu; matrix_t l; // linear part matrix arenstorf_model(double m) : mu(m) { l << 0., 0., 1., 0., 0., 0., 0., 1., 1., 0., 0., 2., 0., 1.,-2., 0.; } void linear_part ( double /* t */, state_t const& y, state_t & dy ) const { dy = l*y; } void non_linear_part ( double /* t */, state_t const& y, state_t & dy ) const { double const y1 = y[0]; double const y2 = y[1]; double const y3 = y[2]; double const y4 = y[3]; double const r1 = sqrt( ( y1 + mu ) * ( y1 + mu ) + y2 * y2 ); double const r2 = sqrt( ( y1 - 1 + mu ) * ( y1 - 1 + mu ) + y2 * y2 ); dy[0] = 0.; dy[1] = 0.; dy[2] = - ( 1 - mu ) * ( y1 + mu ) / ( r1 * r1 * r1 ) - mu * ( y1 - 1 + mu ) / ( r2 * r2 * r2 ); dy[3] = - ( 1 - mu ) * y2 / ( r1 * r1 * r1 ) - mu * y2 / ( r2 * r2 * r2 ); } void operator () (double /* t */ , state_t const& y, state_t & dy ) const { double const y1 = y[0]; double const y2 = y[1]; double const y3 = y[2]; double const y4 = y[3]; double const r1 = sqrt( ( y1 + mu ) * ( y1 + mu ) + y2 * y2 ); double const r2 = sqrt( ( y1 - 1 + mu ) * ( y1 - 1 + mu ) + y2 * y2 ); dy[0] = y3; dy[1] = y4; dy[2] = y1 + 2 * y4 - ( 1 - mu ) * ( y1 + mu ) / ( r1 * r1 * r1 ) - mu * ( y1 - 1 + mu ) / ( r2 * r2 * r2 ); dy[3] = y2 - 2 * y3 - ( 1 - mu ) * y2 / ( r1 * r1 * r1 ) - mu * y2 / ( r2 * r2 * r2 ); } }; Eigen::Matrix mexp ( Eigen::Matrix const& M ) { double tau = M(0,2); double s = std::sin(tau), c = std::cos(tau); Eigen::Matrix e; e << tau*s+c , -tau*c+s , tau*c , tau*s , tau*c-s , tau*s+c , -tau*s , tau*c , tau*c , tau*s , -tau*s+c , tau*c+s , -tau*s , tau*c , -tau*c-s , -tau*s+c ; return e; } int main() { using state_t = Eigen::Matrix; using namespace ponio::observer; double mu = 0.012277471; auto model = arenstorf_model(mu); state_t yini = { 0.994, 0., 0., -2.00158510637908252240537862224 }; double dt = 2e-5; double tf = 17.0652165601579625588917206249; ponio::solve(model, ponio::runge_kutta::rk_118(), yini, {0.,tf}, dt, "arenstorf_adaptivetimestep_demo/orbit_rk_118_ref.txt"_fobs); std::cout << "explicit embedded method" << std::endl; ponio::solve(model, ponio::runge_kutta::rk54_6m().abs_tol(1e-5).rel_tol(1e-3), yini, {0.,tf}, dt, "arenstorf_adaptivetimestep_demo/orbit_rk56_6m.txt"_fobs); ponio::solve(model, ponio::runge_kutta::rk54_7m().abs_tol(1e-5).rel_tol(1e-3), yini, {0.,tf}, dt, "arenstorf_adaptivetimestep_demo/orbit_rk56_7m.txt"_fobs); ponio::solve(model, ponio::runge_kutta::rk54_7s().abs_tol(1e-5).rel_tol(1e-3), yini, {0.,tf}, dt, "arenstorf_adaptivetimestep_demo/orbit_rk56_7s.txt"_fobs); ponio::solve(model, ponio::runge_kutta::rk65_8m().abs_tol(1e-5).rel_tol(1e-3), yini, {0.,tf}, dt, "arenstorf_adaptivetimestep_demo/orbit_rk65_8m.txt"_fobs); ponio::solve(model, ponio::runge_kutta::rk87_13m().abs_tol(1e-5).rel_tol(1e-3), yini, {0.,tf}, dt, "arenstorf_adaptivetimestep_demo/orbit_rk87_13m.txt"_fobs); std::cout << "Lawson embedded method" << std::endl; auto lmodel = ponio::make_lawson_problem( model.l, [&](double t, state_t const& y, state_t& dy ) { model.non_linear_part(t, y, dy); } ); ponio::solve(lmodel, ponio::runge_kutta::lrk54_6m(mexp).abs_tol(1e-5).rel_tol(1e-3), yini, {0.,tf}, dt, "arenstorf_adaptivetimestep_demo/orbit_lrk56_6m.txt"_fobs); ponio::solve(lmodel, ponio::runge_kutta::lrk54_7m(mexp).abs_tol(1e-5).rel_tol(1e-3), yini, {0.,tf}, dt, "arenstorf_adaptivetimestep_demo/orbit_lrk56_7m.txt"_fobs); ponio::solve(lmodel, ponio::runge_kutta::lrk54_7s(mexp).abs_tol(1e-5).rel_tol(1e-3), yini, {0.,tf}, dt, "arenstorf_adaptivetimestep_demo/orbit_lrk56_7s.txt"_fobs); ponio::solve(lmodel, ponio::runge_kutta::lrk65_8m(mexp).abs_tol(1e-5).rel_tol(1e-3), yini, {0.,tf}, dt, "arenstorf_adaptivetimestep_demo/orbit_lrk65_8m.txt"_fobs); ponio::solve(lmodel, ponio::runge_kutta::lrk87_13m(mexp).abs_tol(1e-5).rel_tol(1e-3), yini, {0.,tf}, dt, "arenstorf_adaptivetimestep_demo/orbit_lrk87_13m.txt"_fobs); std::cout << "additive embedded method" << std::endl; auto amodel = ponio::make_imex_jacobian_problem( [&](double t, state_t const& y, state_t& dy ) { model.non_linear_part(t, y, dy); }, [&](double t, state_t const& y, state_t& dy ) { model.linear_part(t, y, dy); }, [&](double t, state_t const& y ) { return model.l; } ); ponio::solve(amodel, ponio::runge_kutta::ark324l2sa().abs_tol(1e-7).rel_tol(1e-5), yini, {0.,tf}, dt, "arenstorf_adaptivetimestep_demo/orbit_ark324l2sa.txt"_fobs); ponio::solve(amodel, ponio::runge_kutta::ark436l2sa().abs_tol(1e-7).rel_tol(1e-5), yini, {0.,tf}, dt, "arenstorf_adaptivetimestep_demo/orbit_ark436l2sa.txt"_fobs); ponio::solve(amodel, ponio::runge_kutta::ark548l2sa().abs_tol(1e-7).rel_tol(1e-5), yini, {0.,tf}, dt, "arenstorf_adaptivetimestep_demo/orbit_ark548l2sa.txt"_fobs); return 0; } .. parsed-literal:: Writing arenstorf_adaptivetimestep_demo/main.cpp .. code:: ipython3 %system $CXX -std=c++20 -O2 -I ../include -I ${CONDA_PREFIX}/include arenstorf_adaptivetimestep_demo/main.cpp -o arenstorf_adaptivetimestep_demo/main .. parsed-literal:: [] .. code:: ipython3 %system ./arenstorf_adaptivetimestep_demo/main .. parsed-literal:: ['explicit embedded method', 'Lawson embedded method', 'additive embedded method'] .. code:: ipython3 import numpy as np import matplotlib.pyplot as plt meths = { "rk56_6m": "RK(5,6) 6m", "rk56_7m": "RK(5,6) 7m", "rk56_7s": "RK(5,6) 7s", "rk65_8m": "RK(6,5) 8m", "rk87_13m": "RK(8,7) 13m", "lrk56_6m": "LRK(5,6) 6m", "lrk56_7m": "LRK(5,6) 7m", "lrk56_7s": "LRK(5,6) 7s", "lrk87_13m": "LRK(8,7) 13m", "ark324l2sa": "aRK(3,2,4) L2 Sa", "ark436l2sa": "aRK(4,3,6) L2 Sa", "ark548l2sa": "aRK(5,4,8) L2 Sa" } for tag,name in meths.items(): data = np.loadtxt("arenstorf_adaptivetimestep_demo/orbit_{}.txt".format(tag)) t = data[:,0] x = data[:,1] y = data[:,2] plt.plot(x,y,label=name) plt.title("Arenstorf orbit") plt.xlabel("$x$ ($y_1$)") plt.ylabel("$y$ ($y_2$)") plt.legend(loc="upper left", bbox_to_anchor=(1.01, 1.025)) plt.show() .. image:: arenstorf_adaptivetimestep_files/arenstorf_adaptivetimestep_18_0.png .. code:: ipython3 for tag,name in meths.items(): data = np.loadtxt("arenstorf_adaptivetimestep_demo/orbit_{}.txt".format(tag)) t = data[:,0] u = data[:,3] v = data[:,4] plt.plot(u,v,label=name) plt.title("Velocity") plt.xlabel("$\\dot{{x}}$ ($y_3$)") plt.ylabel("$\\dot{{y}}$ ($y_4$)") plt.legend(loc="upper left", bbox_to_anchor=(1.01, 1.025)) plt.show() .. image:: arenstorf_adaptivetimestep_files/arenstorf_adaptivetimestep_19_0.png Since :math:`y_1` and :math:`y_2` are linear in the solved problem, we observe fewer differences, contrary to the speed. The exact solution is a closed orbit (so near to the orange curve done by RK(5,6) 7m method). Ponio save also the time step of each iteration (even failed iterations). .. code:: ipython3 for tag,name in meths.items(): data = np.loadtxt("arenstorf_adaptivetimestep_demo/orbit_{}.txt".format(tag)) t = data[:,0] dt = data[:,-1] plt.plot(t,dt,"+-",label=name) plt.xlabel("time") plt.ylabel("time step: $\\Delta t$") plt.legend(loc="upper left", bbox_to_anchor=(1.01, 1.025)) plt.show() .. image:: arenstorf_adaptivetimestep_files/arenstorf_adaptivetimestep_22_0.png .. code:: ipython3 from scipy.interpolate import interp1d data_ref = np.loadtxt("arenstorf_adaptivetimestep_demo/orbit_rk_118_ref.txt") t_ref = data_ref[:,0] x_ref = data_ref[:,1] y_ref = data_ref[:,2] u_ref = data_ref[:,3] v_ref = data_ref[:,4] fx = interp1d(t_ref,x_ref,kind='cubic') fy = interp1d(t_ref,y_ref,kind='cubic') fu = interp1d(t_ref,u_ref,kind='cubic') fv = interp1d(t_ref,v_ref,kind='cubic') def error(u,v): return np.sqrt(np.square(u-v))/len(u) fig, axs = plt.subplots(2, 1) for tag,name in meths.items(): data = np.loadtxt("arenstorf_adaptivetimestep_demo/orbit_{}.txt".format(tag)) t = data[:,0] x = data[:,1] y = data[:,2] u = data[:,3] v = data[:,4] axs[0].plot(t,error(x,fx(t)) + error(y,fy(t)),"+-",label=name) axs[1].plot(t,error(u,fu(t)) + error(v,fv(t)),"+-",label=name) axs[0].set_ylabel("error on position") axs[1].set_yscale('log') axs[1].set_xlabel("time") axs[1].set_ylabel("error on velocity") axs[0].legend(loc="upper left", bbox_to_anchor=(1.01, 1.05)) plt.show() .. image:: arenstorf_adaptivetimestep_files/arenstorf_adaptivetimestep_23_0.png As we already see on the velocity curve, the method that done the lower error is the RK(5,6) 7m method (orange curve). .. code:: ipython3 n_iters = [] n_success = [] for tag,name in meths.items(): data = np.loadtxt("arenstorf_adaptivetimestep_demo/orbit_{}.txt".format(tag)) t = data[:,0] # number of iteration is done by length of of time but we get failed and successed iterations n_iters.append(len(t)) # iterations where current time are equals are failed iterations n_success.append(len(t[1:][t[:-1]-t[1:] != 0])) n_iters = np.array(n_iters) n_success = np.array(n_success) min_idx = np.argmin(n_iters) max_idx = np.argmax(n_iters) plt.bar(meths.values(),n_iters,label="Number of iterations") plt.bar(meths.values(),n_success,fill=False, hatch='//',label="Number of success iterations") plt.xticks()[1][min_idx].set_color("#2ecc71") plt.xticks()[1][max_idx].set_color("#e74c3c") plt.xticks(rotation=65) plt.legend(loc='upper left') plt.show() .. image:: arenstorf_adaptivetimestep_files/arenstorf_adaptivetimestep_25_0.png .. code:: ipython3 percent = 100*n_success/n_iters min_idx = np.argmin(percent) max_idx = np.argmax(percent) plt.bar(meths.values(),percent) plt.xticks()[1][min_idx].set_color("#e74c3c") plt.xticks()[1][max_idx].set_color("#2ecc71") plt.xticks(rotation=65) plt.ylabel("perecent of success iterations") plt.show() .. image:: arenstorf_adaptivetimestep_files/arenstorf_adaptivetimestep_26_0.png