differential equations - Initial value problem for a system of ODEs solver C program

Question

So I wanted to implement the path of the Moon around the Earth with a C program. My problem is that you know the Moon's velocity and position at Apogee and Perigee. So I started to solve it from Apogee, but I cannot figure out how I could add the second velocity and position as "initial value" for it. I tried it with an if but I don't see any difference between the results. Any help is appreciated!

Here is my code:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

typedef void (*ode)(double* p, double t, double* k, double* dk);

void euler(ode f, double *p, double t, double* k, double h, int n, int N)
{
    double kn[N];
    double dk[N];
    double Rp = - 3.633 * pow(10,8); // x position at Perigee

    for(int i = 0; i < n; i++)
    {
        f(p, 0, k, dk);
        for (int j = 0; j < N; j++)
        {
            if (k[0] == Rp)     // this is the "if" I mentioned in my comment
                                // x coordinate at Perigee
            {                    
                k[1] = 0;   // y coordinate at Perigee
                k[2] = 0;   // x velocity component at Perigee
                k[3] = 1076; // y velocity component at Perigee
            }
            kn[j] = k[j] + h * dk[j];
            printf("%f ", kn[j]);
            k[j] = kn[j];
        }
        printf("
");
    }
}

void gravity_equation(double* p, double t, double* k, double* dk)
{
    // Earth is at the (0, 0)

    double G = p[0]; // Gravitational constant
    double m = p[1]; // Earth mass
    double x = k[0]; // x coordinate at Apogee
    double y = k[1]; // y coordinate at Apogee
    double Vx = k[2]; // x velocity component at Apogee
    double Vy = k[3]; // y velocity component at Apogee
    dk[0] = Vx;
    dk[1] = Vy;
    dk[2] = (- G * m * x) / pow(sqrt((x * x)+(y * y)),3);
    dk[3] = (- G * m * y) / pow(sqrt((x * x)+(y * y)),3);

}

void run_gravity_equation()
{
    int N = 4;  // how many equations there are

    double initial_values[N];
    initial_values[0] = 4.055*pow(10,8); // x position at Apogee
    initial_values[1] = 0; // y position at Apogee
    initial_values[2] = 0; // x velocity component at Apogee
    initial_values[3] = (-1) * 964; //y velocity component at Perigee

    int p = 2; // how many parameters there are

    double parameters[p];
    parameters[0] = 6.67384 * pow(10, -11); // Gravitational constant
    parameters[1] = 5.9736 * pow(10, 24); // Earth mass


    double h = 3600; // step size
    int n = 3000; // the number of steps

    euler(&gravity_equation, parameters, 0, initial_values, h, n, N);
}

int main()
{
    run_gravity_equation();
    return 0;
}

Answer 1

Your interface is

euler(odefun, params, t0, y0, h, n, N)

where

N = dimension of state space
n = number of steps to perform
h = step size
t0, y0 = initial time and value

The intended function of this procedure seems to be that the updated values are returned inside the array y0. There is no reason to insert some hack to force the state to have some initial conditions. The initial condition is passed as argument. As you are doing in void run_gravity_equation(). The integration routine should remain agnostic of the details of the physical model.

It is extremely improbable that you will hit the same value in k[0] == Rp a second time. What you can do is to check for sign changes in Vx, that is, k[1] to find points or segments of extremal x coordinate.

---

Trying to interpret your description closer, what you want to do is to solve a boundary value problem where x(0)=4.055e8, x'(0)=0, y'(0)=-964 and x(T)=-3.633e8, x'(T)=0. This has the advanced tasks to solve a boundary value problem with single or multiple shooting and additionally, that the upper boundary is variable.

---

You might want to to use the Kepler laws to get further insights into the parameters of this problem so that you can solve it just with a forward integration. The Kepler ellipse of the first Kepler law has the formula (scaled for Apogee at phi=0, Perigee at phi=pi)

 r = R/(1-E*cos(phi))

so that

R/(1-E)=4.055e8  and  R/(1+E)=3.633e8, 

which gives

R=3.633*(1+E)=4.055*(1-E)  
==>  E = (4.055-3.633)/(4.055+3.633) = 0.054891,
     R = 3.633e8*(1+0.05489)           = 3.8324e8 

Further, the angular velocity is given by the second Kepler law

phi'*r^2 = const. = sqrt(R*G*m)

which gives tangential velocities at Apogee (r=R/(1-E))

y'(0)=phi'*r = sqrt(R*G*m)*(1-E)/R =  963.9438 

and Perigee (r=R/(1+E))

-y'(T)=phi'*r = sqrt(R*G*m)*(1+E)/R = 1075.9130

which indeed reproduces the constants you used in your code.

The area of the Kepler ellipse is pi/4 times the product of smallest and largest diameter. The smallest diameter can be found at cos(phi)=E, the largest is the sum of apogee and perigee radius, so that the area is

pi*R/sqrt(1-E^2)*(R/(1+E)+R/(1-E))/2= pi*R^2/(1-E^2)^1.5

At the same time it is the integral over 0.5*phi*r^2 over the full period 2*T, thus equal to

sqrt(R*G*m)*T

which is the third Kepler law. This allows to compute the half-period as

T = pi/sqrt(G*m)*(R/(1-E^2))^1.5 = 1185821

With h = 3600 the half point should be reached between n=329 and n=330 (n=329.395). Integration with scipy.integrate.odeint vs. Euler steps gives the following table for h=3600:

 n      [ x[n], y[n] ] for odeint/lsode              for Euler

 328  [ -4.05469444e+08,   4.83941626e+06]    [ -4.28090166e+08,   3.81898023e+07]
 329  [ -4.05497554e+08,   1.36933874e+06]    [ -4.28507841e+08,   3.48454695e+07]
 330  [ -4.05494242e+08,  -2.10084488e+06]    [ -4.28897657e+08,   3.14986514e+07]

The same for h=36, n=32939..32940

 n      [ x[n], y[n] ] for odeint/lsode              for Euler

32938 [ -4.05499997e+08   5.06668940e+04]    [ -4.05754415e+08   3.93845978e+05]
32939 [ -4.05500000e+08   1.59649309e+04]    [ -4.05754462e+08   3.59155385e+05]
32940 [ -4.05500000e+08  -1.87370323e+04]    [ -4.05754505e+08   3.24464789e+05]
32941 [ -4.05499996e+08  -5.34389954e+04]    [ -4.05754545e+08   2.89774191e+05]

which is a little closer for the Euler method, but not much better.