MATHEMATICAL MODELING AND COMPUTER SIMULATION OF A BASIC PROBLEM OF TUBE ARTILLERY EXTERNAL BALLISTICS BY MEANS OF THE MATHCAD SOFTWARE

The paper presents mathematical modeling and computer simulation (MM&CS) in the area of numerical solving of the basic problem of external ballistics for tubed artillery. Five-stage MM&CS scheme for conducting a ballistic simulation is developed. It is shown that a formal mathematical procedure allowing to solve the basic problem of external ballistics is a numerical solution of the Cauchy problem for a system of ballistic differential equations. The trajectory of a projectile flight for the 57-mm ZIS-2 anti-tank cannon is estimated. A solving algorithm and the Mathcad program code are given.The numerical solution for a system of four first order ballistic differential equations is a five-dimensional space.The possibility of visual presentation for a numerical solution was proposed in the form of a square matrix. The boundaries of each subspace are determined. A procedure based on spline functions is developed for checking the correctness of the numerical solution. As a result of such verification, the effects of a light increase in the error at the edges of the integration interval are observed. A comparison of the numerical solution of the basic ballistics problem is conducted by means of “soft“ and “stiff“ solver-functions.The trajectory parameters estimated by “soft“ and “stiff“ methods are the same up to the fifth decimal place.


Introduction
As a physical phenomenon, a tubed artillery shot belongs to the class of fast processes.One of the ways of studying them is based on the method of mathematical modeling and computer simulation which replaces the real physical process by its formalized mathematical description.A problem is thus dealt with a significant amount of calculation, but it can be partially solved by using computer technology and automating calculations due to a program environment.As such, the Mathcad software was selected.
In view of the insignificant number of printed publications related to the technology of conducting external ballistic calculations in Mathcad1 , as well as their non-systemic nature, the author wanted to contribute to filling this gap.To this goal, basic program modules have been developed for making it possible to carry out the MM&CS process for an external-ballistic trajectories assessment and to realize the calculations of flight path parameters.
In this part, the paper is a continuation of (Khaikov, 2018), where a method for estimating a projectile initial velocity as a two-point boundary value problem was proposed.
The aim of this paper is to solve the following subtasks: -the development of a program module for the Cauchy problem solution due to various numerical methods for "soft" systems of ordinary differential equations (ODEs); -the proposal of a method for the correctness verification of the obtained numerical solution; -the application of numerical methods for the solution of the classical basic problem of external ballistics as probably "stiff" ODEs.
All the subtasks must be determined in the Mathcad software.
MM&CS in the area of external ballistics MM&CS (mathematical modeling and computational experiment2 ) in the ballistics of tubed weapons is regarded here as a research methodology and technology, based on applied mathematics, mechanics, programming and computer technologies for: -a deeper understanding of a shot as a physical phenomenon; -developing new and improving existing designs of small firearms and artillery tubed weapons; -developing new and improving existing ammunition; -developing measuring instruments for experimental ballistics goals.
With regard to ballistics problems, the process of MM&CS can be conditionally divided into five stages (Figure 1).Тheir names and essence are listed in Table 1.In the fifth stage, researchers can go back to the previous stages and make the necessary changes (for example, they can change the program-code or refine the initial data).The ways to perform feedback are shown by a dashed line in Figure 1.
Further on, the paper will describe a sequence of stages for assessing unguided projectile external-ballistic trajectories.
Ballistic trajectories and their elements as a result of numerical integrating of a system of ODEs throughthe Mathcad solver-functions The following will be demonstrated: 1) how to calculate a trajectory and evaluate its parameters on the basis of a numerical solution, 2) how to verify the correctness of the obtained numerical solution, and 3) how to solve the basic ballistics (sometimes the word "main" is used) problem of external ballistics by various numerical methods.
The realization of all five stages of the MM&CS process will be shown below by solving problem (1).

Numerical calculation of an external-ballistic flight path and its parameters MM&CS Stage 1. Mathematical model of a projectile external-ballistic flight path.
The longitudinal motion of an artillery projectile in the Earth's atmosphere can be described by the system of ODEs with an independent argument time (t) (Burlov et al, 2006, p.249).This type of mathematical expressions belongs to the type of the Point-Mass Trajectory Model: here 3 x -the abscissa of the trajectory (m); y -the ordinate of the trajectory (m); θ -the angle of a velocity vector relative to the base of the trajectory (degrees or radians); v -the instantaneous projectile velocity (m/s); t -the time of flight (s); g -the acceleration of gravity (m/s 2 ); J -the acceleration of the air resistance force (m/s 2 ).
In order to assess an external-ballistic trajectory, we will perform an example for the flight of the 57-mm Armor Piercing Tracer Solid Projectile 53-BR-271SP (weight 3.14 kg) of the anti-tank gun, ZIS-2 1943 model4 .The ZIS-2 gun is a so-called direct fire weapon and provides shooting with a flat trajectory.

MM&CS Stage 2. Projectile flight path calculation algorithm.
A calculation algorithm is a simple chain of performing computational operations.It is presented in Figure 2. The chain consists of eight elements.In boxes 1, 2, and 6, the values of the constants are determined.Boxes 3 and 4 define individual functions, and box 7 sets a functional matrix-column.The value of a ballistic coefficient is calculated in box 5.
The numerical solution of the Cauchy problem is done in box 8.The difference between box 8 and boxes 1-7 is that box 8 is called the built-in Mathcad solver-function.Data from boxes 1-8 is necessary for the solverfunction for its correct functioning.In the block diagram, Figure 2, it is implied that after box 8 the results (errors) are output in a tabular or graphical form.In accordance with Figure 1, a process for program-code preparing is divided into two parts: the preparation of a database and the development of computational Mathcad procedures.The structure of the initial DB is given in Table 2.When using system (1), 12 constants (variables) and 2 functions are necessary for a numerical solution of the Cauchy problem.
The Mathcad-code program is commented below.Using formula 1 and Table 1, one starts programming the elements of the Cauchy problem.
Let us determine the caliber for an artillery projectile (57 mm = 0.057 m), the weight (3.14 kg) and the value of the i coefficient (according to the quadratic drag law): d:=0.057q:=3.14 i:=4.83 (2) An angle of departure (in radians) is calculated as a set of angular degrees, minutes and seconds according to the formula:

Gradus:=0
Min:=10 Sec:=0 (3) At a point of departure, the value of the acceleration coefficient of gravity is determined as 9.18 m/s 2 .Further on, it is necessary to determine the time interval of integration, its boundaries and the total number of points: As a law of air resistance, we will use the quadratic dependence of the type: where B si = 0.00002694.This simplification is made in order to avoid the complicated mathematical formulas that describe the law of air resistance.
The value of the ballistic coefficient is estimated from formula (3) (Germershausen, 1982, p.159).It is equal to 0.5 m 2 /kg.The acceleration of the air resistance force J will be characterized by form6 The initial conditions (for (1)) are determined as a matrix-column y, which will contain their known numerical values: In view of the fact that the initial velocity of the 57 mm armor-piercing solid tracer projectile 53-BR-271SP is 990 m/s (Burlov et al, 2006, p.425), the matrix-column y will look like: . 990 10 909 . 2 The matrix-column D is the right-hand part of the system of ODEs (1).It includes variables with the following notation: the abscissa of the trajectory x-y 0 ; the ordinate of the trajectory y-y 1 ; the angle of the inclination of the tangent θ -y 2 , and the instantaneous projectile velocity v -y 3 :  Result1:= rkfixed(y, t beg , t end , n point , D).
For solving differential equations, the solver-function rkfixed implement the fourth order Runge-Kutta methods with a fixed step.
The variable Result1 is the matrix containing the results of the numerical solution of the system of ODEs (1).In this case, the matrix has a size of 5 by 1001 elements and contains 5005 numbers.Five rows of the Result1 matrix are the independent argument time (t) and the elements of the matrix y (or D).The first column (argument time) counts 1001 numbers.This number includes the initial conditions plus the numbers of the points in the integration interval.
A special feature of the basic problem of external ballistics is that the system of differential equations (1) contains the quantity J, which does not have an exact analytic expression.The fact is that J is a empirical function which can be described as where C -a ballistic coefficient estimated by the Siacci mathematical expression (Germershausen, 1982, p.159), m 2 /kg; H(y) -a function characterizing the dependence of air density vs altitude; F(v) -a function characterizing the dependence of the air resistance vs the instantaneous projectile velocity.
The functions H(y) and F(v) can be described by means of: -an approximate analytical function; -a piecewise function; -a spline function.
For example, the 1943 air drag law7 can be described by a piecewise function (Konovalov, 1979, p.84 By using the drag functions C x (M) the quantity J will be calculated as where 0 a -the speed of sound constant (about 340 m/s).

MM&CS Stage 4. Calculation.
An example of the numerical solution of ODEs (1) as a Result1 matrix is shown in Figure 4.Note that the Mathcad software independently enumerates the columns and rows.The vertical column and the horizontal row with numbers and gray background should not be confused with numerical data.
The first row contains the initial conditions for (1).The leftmost column is the independent argument (in this case, the time of flight).The values of the first column are limited by the boundaries of the integration interval.
In order to use the found values independently from the matrix Result1 will create 6 variables; their purpose are explained in Table 3.If necessary, the variables can be combined into a new matrix.
Usually, the output of tabular data is accompanied by graphs.Their coordinates type and scale must be convenient for further analyses.
An angle of the velocity vector relative to the base of a trajectory in radians 5 θ Gr An angle of the velocity vector relative to the base of a trajectory in degrees 6

V
An instantaneous projectile velocity (m/s)

MM&CS Stage 5. Analysis of results.
It is known that the analysis of any two quantities can be expressed analytically, in a tabular form and graphically.The numerical solution of the system of ODEs (1) enables ananalysis of tabular and graphic data.The graphs of functions performed in the Cartesian coordinates for the system of ODEs (1) are shown in Figure 4.
Figure 4 shows the complete set of dependencies for the five quantities which characterize the system of ODEs (1).This image is constructed in the form of a square matrix.The quantities indicated in the upper horizontal line (X, Y,θ Gr , V, T) of Figure 4 are plotted respectively on the ordinate axis of each graph.The quantities showed in the left vertical column of this figure are plotted respectively on the abscissa axis.The total number of dependencies is 25.Five graphs along the matrix diagonal are functions of the magnitude itself.Thus, there are 20 dependencies.From this number, 10 are direct dependencies and 10 are reverse ones.The functions Y-X, V-X, and V-T are the most widespread ones in practice.
Sometimes for the value of Y as a height of the trajectory and as an argument of function, there are 2 points (ordinates) on the graphs (for one Y-value).Such ambiguity characterizes the ascending and descending parts of the trajectory.
A starting point in the form of a square is added to each graph.This is done to explain the dynamics of the projectile flight and it displays the initial state from the point of view of the physical process.
The domain of the numerical solution for the system of ODEs ( 1) is the five-dimensional space (5D) of X, Y, θ Gr , V and T quantities.The boundaries of each of the quantities are shown in Table 4.
If the quantities X, Y, θ Gr and V are the parametric equations from the argument T, then, in the three-dimensional space, they can be expressed by two types of 3D-graphs: X, Y, θ Gr X, Y and V.A parametric equation and the three-dimensional coordinate space give an opportunity to create a better graphic image.
The horizontal range (X, m) 0 1000 2 The height of the trajectory (Y, m) 0 1.1 3 The angle of the velocity vector relative to the horizontal plane, in degrees -0.5 2.909•10 -3 4 The instantaneous projectile velocity (m/s) 820 990 5 The time of flight (t, s) 0 1.1

Verification of the correctness for external-ballistics trajectory parameters obtained as result of a numerical solution
In view of the fact that the existing solver-functions are nothing more than "black boxes", and the information about their implementation algorithm is characterized by incompleteness, it may be necessary to verify the correctness of the obtained solution by the substitution method.
At the first level of this process, we create a matrixline of discrete time values: At the second level, we perform the spline interpolation of the four obtained quantities (X, Y, θ Gr , V) that are the result of the numerical solution of the system of ODEs, and express them as a spline function with the argument t in .Thus, we define the four new functions: X in (t in ), Y in (t in ), θ in (t in ) and V in (t in ).
In this case, a linear spline function is used, but M allows the use of quadratic and cubic spline functions.
For the first differential equation ( 1), namely express the left part and right one in the form of the spline functions.For the left-hand side, we perform numerical differentiation and obtain the function L X (t in )


The right-hand side of R X (t in ) is obtained by multiplying One visualizes the results (both the right and the left part) in the form of a graph in the Cartesian coordinates The graph of the right and left part of the first differential equation of system ( 1) is presented as a decreasing function of time.Figure 5 shows its form for the integration interval 0-1.1 s.In order to find a magnitude of the relative error of their difference, one will look for (in percent): To improve the image, the value of Delta X (t in ) is summed with the constant a = 1•10 -10 .The graph in Figure 6 shows that at the ends of the integration interval there are "edge" effects and the error of numerical solving increases.In a similar way, the left and right parts for the second, third and fourth differential equations of system (1) are compared (Figure 7).
There is a coincidence of the function graphs.The right-hand sides of the second and third equations of the system of ODEs (1) are monotonically decreasing time functions.Their shape is close to linear dependence.In contrast, the right-hand side of the fourth equation is an increasing function.By the adopted initial conditions, their domain is negative.
A visual presentation of the errors of the right-hand and left-hand sides of equations 2-4 of the system of ODEs (1) is shown in Figure 7.
A visual presentation of the relative error of the right and left parts for the second, third and fourth differential equation expressed in percent (integration interval 0-1.1 s.) is shown in Figure 8.The conclusion for paragraph 2.2 is that the verification of the numerical solution can be carried out using spline functions with the help of which the approximation of the numerical solution is performed and the right parts of the systems of differential equations are calculated.The lefthand parts of the system of ODEs are calculated by numerical differentiation.At the final stage, we compared the right and left parts of the chosen interval of integration.For system (1), the error in the numerical solution does not exceed 1•10 -7 %.
Numerical solving of the basic external-ballistics problem by numerical methods for "soft" and "stiff" systems of ODEs From a practical point of view, it is interesting to consider the following problems: I. Comparison of the calculation results obtained through various numerical methods.
II. Determination of the "stiff" properties of the system of ODEs characterizing the external ballistician problems.
Regarding the first problem: doing research for comparing the usefulness of various numerical methods by solving a system of ballistic differential equations will possibly allow an increase in a number of methods that can be conducted for numerical solving procedures.During the research, the level of errors, computational costs and the average solution time will be established.
The Mathcad solver-functions for solving the ODEs for "soft" and "stiff" systems are collected in Table 5.There are 7 solver-functions.The computational experiments conducted by the author showed that the calculations of system (1) performed by "soft" ODEs such as rkfixed, Rkadapt and Bulstoer give identical results up to the 5th digit after the decimal point (there were no experiments with the solverfunction Adams).
Regarding the second problem: a parallel solution of the same ballistic problem by several numerical methods will make it possible to detect the appearance of probable "stiff" properties of the system of ODEs.
System (1) was solved with numerical methods for solving "stiff" ODEs: for instance, the Stiffb solver-function which implements the Bulirsch-Stoer method and the Stiffr solver-function -the Rosenbrock method.
To evaluate a trajectory using "stiff" methods, it is necessary to know the matrix J k that has a size n*(n+1): n rows and n+1 columns.The first column J k contains the partial derivatives dF/dt, while the remaining columns and rows represent the Jacobi matrix dF/dy.The structure of the matrix J k applied to (1) has the following form For D(t,y) (8), the matrix J k (t,y) of size 4*5 (4 rows and 5 columns) The result of testing solutions (1) by the methods of solving a "stiff" system of ODEs Stiffb, Stiffr gives a similar result as "soft" methods.For such solutions, there are 5 digits after the decimal point.

Conclusion
In this article, it is proposed to consider that a numerical solution of the basic problem of external ballistics can be considered as a five-stage process of MM&CS.
It is shown that a formal mathematical procedure allowing to solve the basic problem of external ballistics is a numerical solution of the Cauchy problem for a system of ballistic differential equations.
The flight trajectory of a projectile for the 57-mm ZIS-2 anti-tank cannon is estimated and visualized.A solving algorithm and the Matcad program-code are given.
The numerical solution for a system of four first order ballistic differential equations is a five-dimensional space.The possibility of a 2D visual representation for a numerical solution was proposed in a form of a square matrix.The boundaries of each subspace are determined.
A procedure based on spline functions is developed for checking the correctness of the numerical solution.For verification purposes, spline functions of three kinds can be used: a linear spline, a quadratic spline and a cubic spline.As a result of such verification, the effects of a light increase in the error at the edges of the integration interval are observed.
A comparison of the numerical solution of the basic ballistics problem is conducted by means of "soft" and "stiff" solver-functions.The article considers a possibility of solving a ballistic problem using five types of the Mathcad solver-functions.The trajectory parameters estimated by "soft" and "stiff" methods are the same up to the fifth decimal place.Konovalov, A.A., Nikolayev, Yu.V. 1979

Figure 1 -
Figure 1 -A generalized presentation of the MM&CS process Рис. 1 -Обобщенное представление процесса матмоделирования Слика 1 -Уопштена представа процеса MM&KSTable 1 -Name of the stage and the essence of the performed actions Таблица 1 -Характеристика этапов математического моделирования Табела 1 -Називи фаза и кратки описи предузетих активности № Name of the stage The essence of the performed actions 1 2 3 1 Development of a mathematical ballistic model The construction of a mathematical model as a formal mathematical description of a physical phenomenon or object. 2 Development of a calculation algorithm Selection of numerical methods for the calculation and development of a computational algorithm.3 Choice of the programming environment and the development of a program code Implementation of the developed algorithm in the selected software for engineering calculation.An important element in the development of program code is a definition of the existence area for the initial database (DB) and the development of management tools.
the solver-function rkfixed for numerical solution, (1) is performed by the command:

Table 2 -
Characteristics of the initial database ): Table is constant and it is determined as (in milliseconds)

Table 4 -
Characterization of the boundaries of the quantities (X, Y, θ Gr