ODE求解器C程序系统的初值问题 [英] Initial value problem for a system of ODEs solver C program
问题描述
所以我想用C程序实现月球绕地球的路径。
我的问题是您知道月球在Apogee和Perigee的速度和位置。
因此我开始从Apogee解决它，但是我不知道如何将第二个速度和位置添加为初始值。我使用 if 进行了尝试，但结果之间没有任何区别。谢谢您的帮助！
这是我的代码：
＃包括< 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在近地点
的位置for（int i = 0; i< n; i ++）
{
f（p，0，k，dk）;
for（int j = 0; j< N; j ++）
{
if（k [0] == Rp）//这是我在评论中提到的 if
// x近地点
{
k [1] = 0; // y坐标在近地点
k [2] = 0; // x近地点
k [3]的速度分量= 1076; //近地点的y速度分量
}
kn [j] = k [j] + h * dk [j];
printf（％f，kn [j]）;
k [j] = kn [j];
}
printf（ \n）;
}
}
voidgravity_equation（double * p，double t，double * k，double * dk）
{
//地球在（0，0）
double G = p [0]； //引力常数
double m = p [1]; //地球质量
double x = k [0]; // x在Apogee
double y = k [1]处的坐标； //在Apogee
处的y坐标double Vx = k [2]; // x在Apogee
处的速度分量double Vy = k [3]; //在Apogee
dk [0] = Vx处的y速度分量；
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; //有多少个等式
double initial_values [N];
initial_values [0] = 4.055 * pow（10,8）; // x在Apogee
initial_values [1] = 0处的位置； // y在Apogee处的位置
initial_values [2] = 0; // x在Apogee
处的速度分量initial_values [3] =（-1）* 964; // y近地点的速度分量
int p = 2; //有多少个参数
double参数[p];
参数[0] = 6.67384 * pow（10，-11）; //引力常数
参数[1] = 5.9736 * pow（10，24）; //地球质量
double h = 3600; //步长
int n = 3000; //步数
euler（& gravity_equation，parameters，0，initial_values，h，n，N）;
}
int main（）
{
run_gravity_equation（）;
返回0；
}
您的界面是
euler（odefun，params，t0，y0，h，n，N）
其中
N =状态空间的大小
n =执行步骤的数量
h =步骤大小
t0，y0 =初始时间和值
此过程的预期功能似乎是在数组 y0 内部返回更新后的值。没有理由插入一些hack来迫使状态具有一些初始条件。初始条件作为参数传递。就像在 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 。这具有解决单次或多次射击的边值问题的高级任务，此外，其上限是可变的。
您可能想使用 Kepler定律来进一步了解这个问题，这样您就可以通过正向集成解决它。 第一个开普勒定律的开普勒椭圆具有公式（按 phi = 0 的Apogee缩放， phi =的Perigee缩放pi ）
r = R /（1-E * cos（phi））
因此
R /（1-E）= 4.055e8和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 =常量= sqrt（R * G * m）
给出Apogee的切向速度（ 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 与Euler步骤的集成给出了 h = 3600 的下表：
n [x [n]，y [n]] for欧德特/ lsode for欧拉
328 [-4.05469444 e + 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欧德特/ lsode for欧拉
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]
对于Euler方法而言，它稍微接近一点，但效果却更好。 p>
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("\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;
}
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.
这篇关于ODE求解器C程序系统的初值问题的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
