Movimiento de proyectil de arrastre cuadrático

$\newcommand{\dd}{\mathrm{d}}$Básicamente tienes dos ODE para resolver:

$$\begin{align} \frac{\dd v^\mu}{\dd t}&=\frac{1}{m}F(x^\mu,v^\mu) \tag{1} \\ \frac{\dd x^\mu}{\dd t}&=v^\mu\tag{2} \end{align}$$ que es más o menos el caso de la mayoría de las fuerzas en la mecánica newtoniana. Para resolver esto numéricamente, desea discretizar el espacio y el tiempo . Con un sistema como (1) y (2), realmente solo debemos preocuparnos por dividir el tiempo.

Una de las rutinas más estables no es en realidad RK4, sino una variación de la integración de salto llamada velocidad verlet . Esto convierte (1) y (2) en un proceso de varios pasos:

$$\begin{align} a_1^\mu &= F\left(x^\mu_i\right)/m \\ x_{i+1}^\mu &= x^\mu_i + \left(v_i + \frac{1}{2}\cdot a_1^\mu\cdot\Delta t\right)\cdot\Delta t \\ a_2^\mu & =F\left(x^\mu_{i+1}\right)/m \\ v_{i+1}^\mu &= v^\mu_i + \frac{1}{2}\left(a_1^\mu+a_2^\mu\right)\cdot\Delta t \end{align}$$ que en realidad es un poco fácil de implementar numéricamente, literalmente solo llama a la función para la fuerza y ​​luego actualiza un par de matrices ( x,y,vx,vy).

Donde tu problema difiere es que$a^\mu=a^\mu\left(x^\mu,\,v^\mu\right)$, lo que hace que calcular la segunda aceleración sea un poco complicado ya que$a_2$depende de$v_{i+1}^\mu$y viceversa. Esta respuesta en GameDev (definitivamente vale la pena leer algunos aspectos numéricos del problema) sugiere que puede usar el siguiente algoritmo

$$\begin{align} a_1^\mu &= F\left(x^\mu_i,\,v_i^\mu\right)/m \\ x_{i+1}^\mu &= x^\mu_i + \left(v_i + \frac{1}{2}\cdot a_1^\mu\cdot\Delta t\right)\cdot\Delta t \\ v_{i+1}^\mu &= v_i^\mu + \frac{1}{2}\cdot a_1^\mu\cdot\Delta t \\ a_2^\mu & =F\left(x^\mu_{i+1},\,v_{i+1}^\mu\right)/m \\ v_{i+1}^\mu &= v^\mu_i + \frac{1}{2}\cdot\left(a_2^\mu-a_1^\mu\right)\cdot\Delta t \end{align}$$ aunque el autor de esa publicación afirma,

No es tan preciso como Runge-Kutta de cuarto orden (como cabría esperar de un método de segundo orden), pero es mucho mejor que Euler o Verlet de velocidad ingenua sin la estimación de velocidad intermedia, y aún conserva la propiedad simpléctica de normal velocidad Verlet para fuerzas conservadoras no dependientes de la velocidad.

Dado que este es un movimiento de proyectil,$x=y=0$es probablemente una elección natural para las condiciones iniciales, con$v_y=v_0\sin(\theta)$y$v_x=v_0\cos(\theta)$como es normal