This is a classic example of a system for which an adaptive step size method is

Question

This is a classic example of a system for which an adaptive step size method is useful because for the large periods of time when the comet is moving slowly we can use long time steps, so that the program runs quickly, but short time-steps are crucial in the brief but fast-moving period close to the sun. The differential equation obeyed by a comet is straightforward to derive. The form between the sun, with mass M at the origin, and a cornet of mass m with positions vector r is GMm/r^2 in direction -r/r (i.e., the direction towards the Sun), and here Newton's second law tells us that m d^2r/dt^2 = -(GMm/r^2) r/r, Canceling the m and taking the x component we have d^2x/dt^2 = -GM x/r^3, and similarly for the other two coordinates. We can, however, throw out one of the coordinates because the comet stays in a single plane as it orbits. If we orient our axes so that this plane is perpendicular to the z-axis, we can forget about the z coordinate and we are left with just two second-order equations to solve: d^2x/dt^2 = -GM x/r^3, d^2y/dt^2 = -GM y/r^3 where r = squareroot x^2 + y^2 a) Turn these two second-order equations into four first-order equations, using the methods you have learned b) Write a program to solve your equations using the fourth-order Runge-Kutta method with a fixed step size. You will need to look up the mass of the Sun and Newton's gravitational constant G. As an initial condition, take a comet at coordinates x = 4 billion kilometers and y = 0 (which is somewhere out around the orbit of Neptune) with initial velocity v_x = 0 and v_y = 500ms^-1. Make a graph showing the trajectory of the comet (i.e., a plot of y against x) Choose a fixed step size h that allows you to accurately calculate at least two full orbits of the comet. Since orbits are periodic, a good indicator of an accurate calculation is that successive orbits of the comet lie on top of one another on your plot. If they do not then you need a smaller value of h. Give a short description d)

Explanation / Answer

taking a single 4th order Runge-Kutta step in y with step size

dy= - y(tn + dt)

#include <algorithm>
#include <cmath>
#include <cstdlib>
#include <fstream>
#include <iostream>
#include <string>
using namespace std;

#include "diffeq.hpp"
using namespace cpt;

const double pi = 4 * atan(1.0);
const double G_m1_plus_m2 = 4 * pi * pi;

bool step_using_y = false; // to interpolate to y = 0

// Derivative vector for Newton's law of gravitation
Matrix<double,1> equations(Matrix<double,1>& trv) {
double t = trv[0], x = trv[1], y = trv[2], vx = trv[3], vy = trv[4];
double r = sqrt(x*x + y*y);
double ax = - G_m1_plus_m2 * x / (r*r*r);
double ay = - G_m1_plus_m2 * y / (r*r*r);
Matrix<double,1> flow(5);
flow[0] = 1;
flow[1] = vx;
flow[2] = vy;
flow[3] = ax;
flow[4] = ay;
if (step_using_y) // change independent variable from t to y
for (int i = 0; i < 5; i++)
flow[i] /= vy;
return flow;
}

void print_data(Matrix<double,1> values)
{
char names[][6] = { "t", "x", "y", "v_x", "v_y" },
units[][6] = { "yr", "AU", "AU", "AU/yr", "AU/yr" };
for (int i = 0; i < 5; i++)
cout << " " << names[i] << " " << values[i]
<< " " << units[i] << " ";
}

void write_data(ofstream& file, Matrix<double,1> values)
{
for (int i = 0; i < values.size(); i++)
file << values[i] << " ";
file << " ";
}

void integrate(Matrix<double,1>& trv, double dt, double t_max,
double accuracy=1e-6, bool adaptive=false)
{
string file_name("comet_fixed.data"), f_or_a("fixed");
if (adaptive) {
file_name = "comet_adapt.data";
f_or_a = "adaptive";
}
ofstream file(file_name.c_str());

cout << " Initial conditions: ";
print_data(trv);
cout << " Integrating with " << f_or_a << " step size ...";
int step = 0;
double dt_min = dt, dt_max = dt;

while (true) {
write_data(file, trv);
double y_save = trv[2];
step_using_y = false;
if (adaptive) {
dt = RK4_adaptive_step(trv, dt, equations, accuracy);
dt_min = min(dt, dt_min);
dt_max = max(dt, dt_max);
} else
RK4_step(trv, dt, equations);
double t = trv[0], x = trv[1], y = trv[2], vx = trv[3], vy = trv[4];
if (x > 0 && y * y_save < 0) {
step_using_y = true;
RK4_step(trv, -y, equations);
write_data(file, trv);
break;
}
if (t > 10 * t_max) {
cerr << " t too big, quitting ..." << endl;
break;
}
step += 1;
}

file.close();
cout << " Number of " << f_or_a << " steps = " << step << endl;
if (adaptive)
cout << " Minimum dt = " << dt_min << " "
<< " Maximum dt = " << dt_max << endl;
print_data(trv);
cout << " Trajectory data in " << file_name << endl;
}

int main()
{
cout << " comet orbit using fixed and then adaptive Runge-Kutta "
<< " Enter aphelion distance in AU: ";
double r_aphelion, eccentricity, dt, accuracy;
cin >> r_aphelion;
cout << " Enter eccentricity: ";
cin >> eccentricity;
double a = r_aphelion / (1 + eccentricity);
double T = pow(a, 1.5);
double vy0 = sqrt(G_m1_plus_m2 * (2 / r_aphelion - 1 / a));
cout << " Semimajor axis a = " << a << " AU "
<< " Period T = " << T << " yr "
<< " v_y(0) = " << vy0 << " AU/yr "
<< " Enter step dt: ";
cin >> dt;
cout << " Enter desired accuracy for adaptive integration: ";
cin >> accuracy;

Matrix<double,1> trv(5);
trv[0] = 0; trv[1] = r_aphelion; trv[2] = 0; trv[3] = 0; trv[4] = vy0;
integrate(trv, dt, T, accuracy, false);
trv[0] = 0; trv[1] = r_aphelion; trv[2] = 0; trv[3] = 0; trv[4] = vy0;
integrate(trv, dt, T, accuracy, true);

}

Related Questions

Navigate

• Browse (All)
• Browse T
• Subjects

---

---

Integrity-first tutoring: explanations and feedback only — we do not complete graded work. Learn more.