A Convenient Approximation of the Projectile Motion with Quadratic Air Resistance

A convenient approximated analytic solution is proposed for the problem of the motion of a body under a resistive force, acting in the magnitude of the squared velocity of the body. This solution is an explicit function of time, that keeps a good behavior both near the initial state and far from the initial state. To obtain a general analytic solution, we firstly used a reduction principle to be able to manipulate scalar objects, and we analyzed limit behaviors, both near the initial state and far from the initial state. Secondly, we proposed an approximated analytic solution with heuristics based on the built knowledge. Finally, a robust and stable integration scheme is proposed, based on the obtained analytic solution. We compared the scheme with other standard integration schemes.


Introduction
One of the classic problems of physics concerns an object moving through a resistive medium. We focus on two particular laws of the resistive force, either linear or quadratic in the relative velocity in respect with the medium. Both laws are approximation in specific regimes of the motion, as detailed in any textbook on hydrodynamics. Solutions with a linear and a quadratic resistive force have already been obtained. The case of a motion with a linear resistive force is explicitly solvable, by meaning that the exact solution can be expressed by elementary functions of time. The case with a quadratic resistive force is not explicitly solvable. We call 'elementary functions' all functions that can be locally expressed with convergent power series.
Our aim is to provide a stable integration scheme of the projectile motion, where 'stable' means that the scheme is unconditionally stable in regards with the time increment (also called delta-time or time-step). This scheme could be used for real-time applications for which the time-step is unbounded. Also, it could be used to handle pre-warmed simulations, where simulations start at a non-zero time. (For example, loading a game level with smoke stacks already fully grown).

Previous work
Implicit analytic solutions have already been proposed. The solution described in [1] parameterizes the velocity magnutide with direction angle of the velocity. Even when the proposed solutions can be expressed by elementary functions, they are expressed with parametric variables which are built from trajectory parameters. Therefore, such solutions cannot be used to implement time-based integration schemes.
Explicit solutions have already been proposed using various approximation methods. The solutions found by [2] are valid for a nearly flat trajectory, and only near the initial state or far from the initial state. More general solutions have also been found in the form of Taylor series near the initial state, by [3]. Lastly, [4] provides solutions as a ratio of two series expansions near the initial state. However, time-based integration schemes implemented built from series expansions won't be unconditionally stable.

Equation of motion
From the first Newton law, the time derivative of the motion quantum is equal to the resultant of the forces applied on the body. We consider a resultant of the forces that is composed of a constant and uniform force, and a resistive force (also called drag force) that depends on the velocity. We also consider that the mass of the body is constant, and we reduce the body to its center of inertia and its mass. We obtain the following well-known motion equation: where m is the non-null mass of the body, v its velocity, Fc a constant and uniform force and F drag the drag force. Fc often includes the gravity action, and it is replaced by Q such as Fc = m Q. We obtain the general form of the equation, where the varying terms are placed on the left-hand side and the constant terms on the right-hand side:

Analytic solution with a linear drag
In this section, we recall common results when the drag force follows a linear law. The insertion of linear drag force (2) into (1) gives: As the wind is constant and uniform, this equation can be rewritten with the relative velocity u = v− v wind : given the initial condition u(t 0 ) = u 0 .
Solution of the degenerated case With α = 0, the equation admits a solution u(t) = u 0 + (t − t 0 ) Q. By removing the drag force, the motion is the well-known motion of an uniformly accelerated body in a Newtonian frame.
As the velocity is a function that can be analytically integrated, the position can be expressed as a function of time. The exact analytic solution can be used in the integration scheme. Consequently, the scheme is exact, and it obviously does not diverge when the time increment becomes large. When α tends toward zero, the position can be expressed as a Taylor expansion in order to prevent numerical instabilities.
Scheme (integration scheme with a linear drag) Given the position xn and the velocity vn at the time tn, the new position and velocity at time t n+1 = tn + δt are: Unlike with the linear drag, an exact analytic solution that can be expressed with elementary functions of time does not exist, to our knowledge so far. We will first deduce some properties of the solution. Secondly, we will express the solution in an implicit form, by meaning that the solution is expressed with the solution itself. The equation of motion (1), with the quadratic drag force (3), becomes: As the wind is constant and uniform, this equation can be rewritten with the relative velocity u = v − v wind : given the initial condition u(t 0 ) = u 0 .

Properties
In this section, the properties of the solution of the differential equation (5) will be detailed. This will help to restrict approximations, and to guarantee their behavior. The motion without drag (α = 0) admits a solution, already described in section 1.4. So we enforce α > 0. Firstly, we can restrict the relative velocity u to be in a 2D space. Indeed, the variation of the relative velocity is a weighted sum of the vectors Q and itself. Consequently, the velocity stays in the plane formed by Q and u 0 . Moreover, if Q is null, then the velocity will stay in a 1D space. This case is detailed in the section 3. Then, we consider here that Q is not null. Secondly, we assume that the solution exists, is unique, and admits a finite limit on infinite time.
Stationary solution : The equation admits a stationary solution: We define γ = α Q .
Asymptotic properties of the solution around the initial time : Around the initial state, we can trivially write: 1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65 Asymptotic properties of the solution on infinite time : The solution converges toward u∞ on infinite time. Let's introduce a decomposition of u such as: where f is a non-null function of time that tends toward zero on infinite time, and W is an unit vector ( W (t) = 1). As W is an unit vector, its derivative is perpendicular to W . Because the velocity belongs to a 2D-space, we decompose W such as: W (t) = cos(θ(t)) n∞ + sin(θ(t)) t∞ (8) where θ is a function of time, n∞ the unit vector with the same direction as u∞, and t∞ an unit vector, perpendicular to n∞. We inject (7) into (5): We apply a first order Taylor expansion at +∞, and the constant terms vanish by using (6): We project on W ′ ( W terms vanish), and we factorize by f : We are using (8), and by excluding the equation θ ′ (t) = 0, we obtain: This ordinary differential equation admits a limit on infinite time, and 2 stationary points (within [0, π[). However, only the stationary point θ = 0 is stable. Finally, we can conclude that, either it stays on the unstable stationary point: either it converges toward the stable stationary point:

Implicit form of the solution
An implicit solution can be found by using the reduction principle for certain coupled systems of ordinary differential equations introduced by [4]. Let P denote the scalar function P : t → α u(t) . The equation (5) can be written in the canonized form: u ′ + P u = Q (10) The solution of (10) can be written as: Proof The first step is to find the integral solution of the homogeneous equation u ′ + P u = 0. The second step is to write the solution of the origin equation as the product of the homogeneous solution and an unknown function. The last step is to remind the solution exists and is unique.  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65 Let ϕ : t → e t t 0 P (r)dr denote a scalar function of time. The dependency on the variable t 0 is implicit. Let φ denote the integral of ϕ: φ : t → t t0 ϕ(r)dr. Then the solution (11) can be expressed as: Differential equation on φ : The scalar function φ is a solution of the second-order differential equation: Proof By injecting (12) into (5), or by following the deduction detailed in [4]. ⊓ ⊔ Asymptotic properties of φ at the initial time : From the asymptotic properties of u at the initial time and the equation (13), we can trivially determine the derivatives at t 0 .
Asymptotic properties of φ on infinite time : Given the asymptotic properties of u on infinite time and the expression of u in (12), we can deduce that with γ = α Q and a a constant.
First-all, the equation (19) admits a stationary solution: The analysis of the equation (19) is not detailed here, but we remind important results: the equation admits a unique solution (given an initial condition), either a solution is the stationary solution, either the solution never reaches or crosses u∞. Solutions are monotone and bounded, and they always converge to u∞. For Q = 0, we have the homogeneous equation; the related motion's regime is called the homogeneous regime. For Q > 0, as u tends towards u∞ which is greater than zero, and by the limit definition, we obtain the following statement: From this, we deduct that u will become and stay positive after a certain time. If u 0 is negative, knowing that the solution is monotone, then u will stay negative within a limited time, and will become positive. We call down-to-up regime the time interval in which u is negative, and terminal regime the time interval in which u is positive.

Homogeneous regime
When Q = 0, the solution is: From this solution, we obtain an analytic expression of φ: This solution fits with the limit of the terminal solution when Q → 0. The proof can be done with a Taylor expansion from the expressions of u and φ in the terminal regime.

Degenerated regime
When α = 0, the solution is u(t) = Q(t−t 0 )+u 0 and φ(t) = t. Like the homogeneous regime, the degenerated regime is a limit of the terminal regime when α → 0. We can obtain a truncated Taylor expansion of φ from its expression in the terminal regime:

Approximation of the motion of a point under a quadratic drag
Given that the velocity can be expressed with the function φ thanks to the equation (12), we will propose an approximation of φ. One of our main outcome is to write φ as a sum of exponential terms of time. Moreover, the approximation can match with the exact solution when the motion is in the one-dimensional terminal-regime. We can have up to 7 parameters to build the expression of φ since we gathered 7 relations: 5 relations near the initial state with the equations (14), (15), (16) and 1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65 (17), plus 2 relations far from the initial state with the equation (18). The general form is proposed as: With this form, the previous properties of φ can be more detailed: -Equations (14) and (15) are essential as they are deducted from the definition of φ. It leads to the following system: -Equation (16) allows that the velocity obeys to the equation of motion at the initial state. It leads to: -Equation (17) allows even more accurate motion near the initial state, but the coefficients will be more sensitive to the initial state. It leads to: (18) is in fact composed of 2 relations: the exponential coefficient value allows that the velocity converges to the terminal velocity far from the initial state, while the constant coefficient value impacts the way the velocity converges to the terminal velocity. It leads to:

Identification of the parameters
Without loss of generality, we fix Q = −q ey and u 0 = u 0 (sin(θ 0 ) ex − cos(θ 0 ) ey) where ex and ey are 2 perpendicular unit vectors and q is a positive scalar. There is multiple ways to identity the parameters of φ G . We have enough expressions to directly solve the non-linear system formed by equations from (21) to (25), but this strategy leads to spurious results as it tries to equalize an approximated form with asymptotic constrains. Therefore, we used error-minimization techniques.
To do so, we used a curve-fitting method to fit φ G with reference data obtained from a numeric integration (that we call φ ref ). Equality-constrains are implemented using Lagrange-multipliers technique, while inequality-constrains are implemented using interior penalization technique. We also choose γ 3 < γ 2 to distinguish both parameters .   1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63   From the obtained parameters values, as shown in figures 1 and 2, we can either create a 3-dimensional tabulation, or fit the data with functions. In this paper, we decided to find functions that can approximate the values of γ 2 and γ 3 , and to deduct other parameters from the equations. To build the functions, we did some observations: firstly, the exponential coefficients γ 2 and γ 3 are close to be proportional to γ 1 . Secondly, with help of the knowledge of the one-dimensional terminal-regime (in section 3.2), we can identify one of the two exponential coefficients to −γ 1 . Thirdly, the exponential coefficients must converge to a finite value when u 0 tends toward infinity. We finally end with the following empirical   1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65 expressions of the exponential coefficients: Then, other parameters (a 1 , a 2 , a 3 ) can be deducted and computed from the linear-system:   1 1 1 We can recognize here the Vandermonde matrix, which is inversible when γ 1 , γ 2 and γ 3 have distinct values. The function φ G , approximation of φ, is now fully determined. From the implicit solution (12), the approximated velocity is: The figure 3 shows the velocity path during the motion, from different initial velocity. Hopefully, the velocity converges towards the same value, and we remark also that the velocity converges with a certain tangent. We proved above that this tangent is perpendicular with the constant right-hand side vector Q (so ey here). While we obtained good approximation at both asymptotic limits, the in-between approximation is less accurate. The inaccuracy is getting stronger when the motion spends more time in the down-to-up regime (so when cos(θ 0 ) is closer to −1).

Degenerated motions
The present method will fail when λ → 0. This scenario can happen when either α → 0, either Q → 0. In the first case, the terminal velocity diverges. In the second case, the motion tends to the one-dimensional motion, where the velocity path does no longer tends toward the terminal velocity with the same tangent. In both limit cases, the function φ folds into a polynomial in t. We can find a truncated Taylor expansion on φ, with help of the first derivatives of φ from equations (14) to (17):

Approximation of the position
Whereas the velocity has an analytic expression, its integral (the position) cannot be expressed by elementary functions of time. A first approach is to use numeric integral schemes, that will still benefit from our stable approximation of the velocity. Fortunately, as shown in the figure 4, our approximation is still relevant in terms of accuracy, combined with multi-step integral schemes. We can illustrate such schemes with the following mixed scheme: Scheme (mixed integration scheme with a quadratic drag (Ap-E2)) Given the position xn and the velocity vn at the time tn, the new position and velocity at time t n+1 = tn + δt are: .. (as defined by (28)) Our reference integration scheme is using a 4th order Runge-Kutta integration for the velocity, and a 2nd order Euler integration for the position: Scheme (reference integration scheme with a quadratic drag (RK4-E2)) Given the position xn and the velocity vn at the time tn, the new position and velocity at time t n+1 = tn + δt are: 2   1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65 One of our goals is to provide a stable integration, and to minimize the drift of the position when the time increment becomes large. To do so, we will approximate again the velocity to be able to have an analytical expression of the position. The expression (28) is not analytically integrable because of the denominator form. The denominator is a sum of three exponential terms, in which one becomes dominant compared to the two others. Then, it is approximated by: With this approximation, the integral of u G can be analytically expressed. Finally, the proposed integration scheme is the following: Scheme (integration scheme with a quadratic drag (Ap-VP)) Given the position xn and the velocity vn at the time tn, the new position and velocity at time t n+1 = tn + δt are: v n+1 = v wind + u + Q(a 0 + a 1 e γ1t + a 2 e γ2t + a 3 e γ3t ) (a 1 γ 1 e γ1t + a 2 γ 2 e γ2t + a 3 γ 3 e γ3t ) We compare the different schemes by comparing the given position of the body after 10 seconds, after proceeding by successive time-steps with a fixed delta-time .   1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65 As shown in the figure 4, our scheme (Ap-VP) is able to give the final position in one step (delta-time is 10s), with roughly the same precision than the precision obtained by proceeding 15 time-steps with the reference scheme (RK4-E2). Even though the proposed scheme (Ap-VP) gives a good accuracy regardless the number of integration steps, it has a heavier computational complexity. Usage of tabulations may be a way to lower the computational complexity.

Conclusions
From a detailed analysis of the equation limits, both near the initial sate and far from the initial state, we were able to build a good knowledge of the behavior of the motion. We preferred to keep a good balance between the analysis depth of the initial state and terminal state, rather than going for more precision near the initial state. Furthermore, it was demonstrated that the velocity converges toward the terminal velocity with a particular tangent.
Using the reduction principle was a key element towards the ability to work with a single auxiliary variable, so we successfully built an approximation of this auxiliary variable, that leads to deduct an approximated velocity. Several explorations of the present work are now possible, with different forms of the auxiliary variable and different ways to compute its coefficients. Also, as the approximation form is similar to the analytic form in the one-dimensional terminal regime (when the velocity goes in the same direction as the terminal velocity), the obtained solution is less accurate with starting regimes.