Coding Runge-Kutta 4 in C++ for a Force proportional to 1/r^2 outputs a trajectory different from scipy.optimise.solve_ivp

Question

I am trying to code an RK4 solver for forces Fx = -x/r^3 and Fy = -y/r^3 where r=sqrt(x^2+y^2) is the distance between the origin and the position. I am comparing my results to that of scipy.optimise.solve_ivp.

For initial values (x, y, vx, vy) = (1, 0, 0, 1), scipy outputs a trajectory of a satellite slowly falling towards the source of the force (origin). However, my code outputs a stable trajectory. Any solutions? I have attached a graph showing the two outputted trajectories, highlighting how they differ.

#include <iostream>
#include <vector>
#include <tuple>
#include <cmath>
#include <fstream>

//double const constantForce_x = 1e-3, constantForce_y = 1e-3;
double const constantForce_x = 0., constantForce_y = 0.;

// Function to compute the central force components
std::pair<double, double> central_force(double x, double y) {
    double r = sqrt(x*x + y*y);
    double Fx = -x / (r*r*r) + constantForce_x;
    double Fy = -y / (r*r*r) + constantForce_y;
    return std::make_pair(Fx, Fy);
}

// Function to perform one step of the RK4 method
std::tuple<double, double, double, double> rk4_step(double x, double y, double vx, double vy,     double dt) {
    auto k1 = central_force(x, y);
    double k1x = k1.first, k1y = k1.second;
    double k1vx = vx, k1vy = vy;

    auto k2 = central_force(x + k1vx * dt / 2, y + k1vy * dt / 2);
    double k2x = k2.first, k2y = k2.second;
    double k2vx = vx + k1x * dt / 2, k2vy = vy + k1y * dt / 2;

    auto k3 = central_force(x + k2vx * dt / 2, y + k2vy * dt / 2);
    double k3x = k3.first, k3y = k3.second;
    double k3vx = vx + k2x * dt / 2, k3vy = vy + k2y * dt / 2;

    auto k4 = central_force(x + k3vx * dt, y + k3vy * dt);
    double k4x = k4.first, k4y = k4.second;
    double k4vx = vx + k3x * dt, k4vy = vy + k3y * dt;

    x += (k1vx + 2*k2vx + 2*k3vx + k4vx) * dt / 6;
    y += (k1vy + 2*k2vy + 2*k3vy + k4vy) * dt / 6;
    vx += (k1x + 2*k2x + 2*k3x + k4x) * dt / 6;
    vy += (k1y + 2*k2y + 2*k3y + k4y) * dt / 6;

    return std::make_tuple(x, y, vx, vy);
}

// Function to simulate trajectory
std::vector<std::pair<double, double>> simulate_trajectory(double x0, double y0, double vx0,     double vy0, double dt, int steps) {
    double x = x0, y = y0, vx = vx0, vy = vy0;
    std::vector<std::pair<double, double>> trajectory;
    trajectory.push_back(std::make_pair(x, y));

    for (int i = 0; i < steps; ++i) {
        std::tie(x, y, vx, vy) = rk4_step(x, y, vx, vy, dt);
        trajectory.push_back(std::make_pair(x, y));
    }

    return trajectory;
}

int main() {
    // Parameters
    double x0 = 1.0;    // initial x position
    double y0 = 0.0;    // initial y position
    double vx0 = 0.0;   // initial x velocity
    double vy0 = 1.0;   // initial y velocity
    double dt = 0.001;   // time step
    int steps = 50000;   // number of steps

    // Simulate trajectory
    auto trajectory = simulate_trajectory(x0, y0, vx0, vy0, dt, steps);

    // Write trajectory to file
    std::ofstream outfile("trajectory.dat");
    for (auto point : trajectory) {
        outfile << point.first << " " << point.second << "\n";
    }
    outfile.close();

    return 0;
}

I have tried to code the RK45 method used my scipy to test if that is the difference, but there was no change.

Answer 1

solve_ivp will solve this problem perfectly well, but you will have to include the optional argument rtol (or atol) to set the intended relative (or absolute) tolerances expected so that the automatic adjustment of timestep that it uses will be adequate.

In the code below I have set rtol=1.0e-6. However, the default value (if you don't include this argument) is 1.0e-3 and, if you try it, you will find that the solution with that is rubbish (and spirals in slowly).

import math
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp


def deriv( t, Y ):
    x, y, vx, vy = Y
    rcubed = ( x ** 2 + y ** 2 ) ** 1.5      # note: r/r^3 --> 1/r^2, so inverse-square law
    return np.array( [ vx, vy, -x / rcubed, -y / rcubed ] )


r = 1.0
v = 1.0 / math.sqrt( r )               # for a circular orbit
period = 2 * math.pi / ( v / r )       # ditto
tmax = 10 * period
Y0 = np.array( ( r, 0, 0, v ) )        # initial x, y, vx, vy
sol = solve_ivp( deriv, [ 0.0, tmax ], Y0, rtol=1.0e-6 )

x = np.array( sol.y[0] )
y = np.array( sol.y[1] )
plt.plot( x, y )
plt.axis( 'scaled' )
plt.xlabel( 'x' );   plt.ylabel( 'y' )
plt.show()