Задача начального значения для системы программы ODE solver C

Question

Итак, я хотел реализовать путь Луны вокруг Земли с помощью программы на Си.Моя проблема в том, что вы знаете скорость и положение Луны в Апогее и Перигее.Так что я начал решать это из Апогея, но я не могу понять, как я мог бы добавить вторую скорость и позицию в качестве «начального значения» для него.Я попробовал это с if, но я не вижу никакой разницы между результатами.Любая помощь приветствуется!

Вот мой код:

#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("\n");
    }
}

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

Ваш интерфейс

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

, где

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

Предполагается, что предполагаемая функция этой процедуры заключается в том, что обновленные значения возвращаются в массиве y0.Нет причин вставлять некоторые хаки, чтобы заставить состояние иметь некоторые начальные условия.Начальное условие передается в качестве аргумента.Как вы делаете в void run_gravity_equation().Процедура интеграции должна оставаться независимой от деталей физической модели.

Крайне маловероятно, что вы достигнете того же значения во k[0] == Rp во второй раз.Что вы можете сделать, это проверить изменения знака в Vx, то есть k[1], чтобы найти точки или отрезки экстремальной x координаты.

---

Попытка интерпретировать ваше описание ближе, чтоВы хотите сделать, это решить краевую задачу, где x(0)=4.055e8, x'(0)=0, y'(0)=-964 и x(T)=-3.633e8, x'(T)=0.Это имеет расширенные задачи для решения краевой задачи с одиночной или многократной съемкой и, кроме того, что верхняя граница является переменной.

---

Возможно, вы захотите использовать законы Кеплера , чтобы получить более полное представление о параметрах этой проблемы, чтобы вы могли решить ее просто с помощью прямой интеграции.Эллипс Кеплера первого закона Кеплера имеет формулу (масштабируется для Апогея в phi=0, Перигея в phi=pi)

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

, так что

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

, что дает

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 

Кроме того, угловая скорость задается вторым законом Кеплера

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

, который дает тангенциальные скорости в Апогее (r=R/(1-E))

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

и Перигея (r=R/(1+E))

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

, который действительно воспроизводит константы, которые вы использовали в своем коде.

Площадь эллипса Кеплера равнаpi/4 раз произведение наименьшего и наибольшего диаметра.Наименьший диаметр можно найти на cos(phi)=E, наибольший - это сумма радиуса апогея и перигея, так что площадь равна

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

. В то же время это интеграл по 0.5*phi*r^2 пополный период 2*T, таким образом, равный

sqrt(R*G*m)*T

, что является третьим законом Кеплера .Это позволяет вычислять полупериод как

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

. При h = 3600 половина должна быть достигнута между n=329 и n=330 (n=329.395).Интеграция с scipy.integrate.odeint против шагов Эйлера дает следующую таблицу для 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]

То же самое для 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]

, что немногоближе к методу Эйлера, но не намного лучше.