¿Por qué mi modelo de dos cuerpos no funciona? [cerrado]

En situaciones como esta, es una buena idea ajustar el paso de tiempo en función del gradiente de la fuerza, porque todo el concepto de integración numérica es que "las cosas no cambian demasiado desde ahora hasta el siguiente paso de tiempo", y esa suposición se viola cuando te mueves rápidamente a través de una región con una fuerza que cambia rápidamente.

Esto tiene un efecto secundario arriesgado: si hace que el paso de tiempo sea "lo suficientemente pequeño" para hacer frente a la variación rápida, y la variación se vuelve "infinitamente rápida", su algoritmo puede "atascarse": puede encontrarse con la paradoja de Xeno de la Aquiles y la tortuga.

Para evitar quedarse atascado, algunas integraciones numéricas (como Runge-Kutta) son mejores para tener en cuenta la curvatura de orden superior, lo que permite un paso más grande sin perder precisión.

Sin embargo, en su caso, la aceleración para el próximo paso de tiempo completo está determinada por la posición actual, por lo que si sus partículas están cerca unas de otras al comienzo del paso, serán arrojadas lejos, y dado que su fuerza de atracción entonces se verán muy disminuidos, les será difícil volver.

Algunas críticas a tu código:

1. Tú defines particle1[4]pero solo usas dos componentes
2. Codificas el paso de tiempo en lugar de usar una variable comodt
3. Usas el método de integración más básico... por favor aprende sobre otros
4. Usted define implícitamente masa = 1 y constante de fuerza = 1; considere hacer esas variables (incluso si las establece en 1)

Con respecto al primer punto, usaría a structpara mis partículas:

typedef struct{
  double x;
  double v;
} PARTICLE;

y luego puedes acceder a sus propiedades con

PARTICLE p1, p2;
p1.v = 0.;
p1.x = 0.;
p2.v = 0.;
p2.x = 1.;

Lo cual es inmediatamente más legible... Si desea trabajar en dos dimensiones, puede escribirlas explícitamente como parte de su struct: ( double x; double y; double vx; double vy;) o convertirlas en matrices.

Pero en realidad, los pasos de tiempo variables y las interpolaciones de mayor orden ayudarían con la estabilidad de su solución.

Estás obteniendo partículas que tienen una distancia de 0 entre sí y, por lo tanto, tienen una aceleración infinita. La computadora no puede resolver esto, por lo que explota según lo cerca que esté el paso de cero.

No estoy seguro de cuál sería el comportamiento "físico", dado que las masas puntuales de distancia cero son difíciles de estudiar. Pero si agregó una coordenada y y le dio una velocidad inicial, podría obtener mejores resultados (básicamente, una elipse alargada en lugar de una línea oscilante).

Puede hacer que la velocidad y inicial sea muy pequeña, de modo que se aproxime al movimiento X. O tal vez solo agregue una pequeña compensación a su división.