**1. Introduction**

Mathematical models are mostly used in natural sciences (physics, chemistry, earth science) and engineering disciplines (computer science, biological science, genetic engineering), as well as in social sciences. Substantial stages of implementation (creation, design, simulation, and visualization) of a continuous mathematical model of aircraft motion in a flying simulator were described in many papers. The one-parameter control of a distributed mathematical model of motion of aircraft with all phases of implementation in a flying simulator was described in the published paper [1, 2].

#### *Aeronautics - New Advances*

Motivation of this research is using a mathematical model of an object, which in practice and the real world often can be determined analytical otherwise. For this reason, it is necessary to design mathematical models characterized for this problem statement. The more accurate the mathematical model of the object we create, the more accurate the results of its simulation.

The mathematical model of motion of aircraft can be simulated using central computer architecture, which can be based on single-processor systems. The simulation of structured mathematical model of an aircraft can be implemented using central computer architecture, a distributed memory system (DMS) using a message interface (MPI) architecture; they are identified as node computers or shared memory system (SMS) and graphics processing unit (GPU) architectures [3]. The parallel computer architecture may be based on multiprocessors; each processor has a multicore architecture. Multiprocessor systems are computationally more efficient than such systems compared with the central computer architecture, see [4].

As known from specialized literature, creation of mathematical models of motion of aircraft consists of the following phases: definition of a physical base for creation of a mathematical model, selection of notation of a mathematical model, implementation of a mathematical model and others. The Laplace transformation, computation of parameters of aircraft for a selected flight phase, determination and computation of coefficients, numerical integration, programming, and simulation on a computer, etc. are used in our problem statement and proposed solution.

According to current knowledge, the effort made to create mathematical models of motion is aimed at improvement and extension of a field of view of simulation of mathematical models on the different computer architectures. The solution of large and complex mathematical models of aircraft problems leads to considerable demand for parallel computing and strategies [5]. We can confirm the most frequent use of the following phases: creation, design, implementation, and simulation of mathematical models. The parallel computing tools guarantee an impressive reduction of computing time [6].

A block structure mathematical model of motion of aircraft can be integrated numerically using parallel computer architecture. This architecture can be motivation for using a central processing unit (CPU)-based implementation with multiprocessor systems and a graphics processing unit (GPU); see [7]. Likewise, a computer system fitted with a GPU accelerator and a software application designed and created for this architecture is also suitable for simulation of a block mathematical model of aircraft motion. From our point of view, the following three principles are important for proposed solution:


Accordingly, in the chapter, we explain how to implement a compute unified device architecture (CUDA) algorithm of computation of a block structure mathematical model of motion of aircraft so as to achieve peak speed. The efficiency is

*Simulation of a Mathematical Model of an Aircraft Using Parallel Techniques (MPI and GPU) DOI: http://dx.doi.org/10.5772/intechopen.105538*

demonstrated by discussing two-parameter control of aircraft, throttle control stick (fuel supply), and aircraft elevator control stick (angle of attack).

#### **2. Description of aircraft mathematical model**

The unstable characteristics such as aero elasticity impact, fuel density, changing geometry of aircraft, and some other parameters support complexity of our design mathematical models in simulation. The mathematical models of motion interact with intervening pilot's control of aircraft and real equipment respond to the pilot's interventions, data on equipment the flying aircraft are observed by the pilot. We can use a continuous simulation method to solve differential equations in mathematical models of motion of aircraft created this way. Aerospace engineers often use Newton's laws of motion in design and creation of a mathematical model of motion of aircraft described by differential operators.

According to Newton's third law, there is always the same opposite reaction to every action, or the interaction of two bodies on each other always leads to opposite parts. The basic system of equations has the form [8]:

$$\dot{\mathbf{x}}\_1 + f\_i(\mathbf{x}\_1, \dots, \mathbf{x}\_n; u\_1, \dots, u\_m; \xi\_1, \dots, \xi\_\gamma) = \mathbf{0}, (i = 0, 1, \dots, p) \tag{1}$$

where *x*1, … , *xn* are object coordinates, *u*1, … , *um* are elements of control, and ζ1, … , ζ*<sup>γ</sup>* are failure functions. The problems of simulation and synchronization of mathematical models on computers is know, where the accuracy of a linear motion model is used [9]. The accuracy between the transient response results of nonlinear mathematical models compared with linear mathematical models from in terms of human accuracy can be neglected [10].

#### **2.1 Equations numerical integration of mathematical model**

The proposed solution based on the presented analysis of the state model space does not contain any discontinuity in either *fi(x, u, t)* or in any of higher derivatives, *xi*(*t*) itself is a continuous function of time. Currently used methods of numerical integration are used to find numerical approximations to solutions of ordinary differential equations (ODEs). The name, also known as numerical integration, represents a broad family of algorithms for calculating the numerical value of a particular integral. The methods of numerical integration can generally be described as combining evaluations of the integrand to get an approximation to the integral. The integrand is evaluated on a finite set of points called integration points, and the weighted sum of these values is used to approximate the integral [11]. This method of calculation is often used for polynomial functions, and these are suitably used on computers. Our proposed solution deals with creation of such functions, i.e., integration of mathematical models.

Two most basic algorithms, Euler (forward and backward) and Taylor series, are motivation for this research in our solution for numerical integration. A numerical stability domain is introduced as a pillar to characterize an integration algorithm and a general procedure to find the numerical stability domain of any integration scheme, see [11]. The mathematical model in state space has the following form:

$$
\dot{\mathbf{x}}(t) = \mathbf{f}(\mathbf{x}(t), \mathbf{u}(t), t), \tag{2}
$$

where *x* is a state vector, *u* is an input vector, and *t* represents time, with a set of initial conditions:

$$\mathbf{x}\_{(t=t\_0)} = \mathbf{x}\_{\mathbf{0}}.\tag{3}$$

Let *xi(t)* represents the *i*th state trajectory as a function of simulated time *t*. As long as the state-space model does not contain any discontinuity in either *fi*(*x, u, t*) or any of higher derivatives, *xi*(*t*) is itself a continuous function of time. Such function can be approximated by any desired precision by a Taylor series expansion of any given point along its trajectory. Let *t \** denote a point in time, about which we wish to approximate the trajectory using Taylor series, and let *t \* + h* be the point in time, at which we wish to evaluate the approximation [11]. The trajectory at that point can then be given as follows:

$$\varkappa\_i(t^\* + h) = \varkappa\_i(t^\*) + \frac{d\varkappa\_i(t^\*)}{dt} \cdot h + \frac{d^2\varkappa\_i(t^\*)}{dt^2} \cdot \frac{h^2}{2!} + \cdots. \tag{4}$$

Different integration algorithms differ in how they approximate the higher-order derivatives and in the number of Taylor series members that they take into account in the approximation [10].

If the term *n +* 1 of the Taylor series is considered, the approximation accuracy of the second derivative *d2 xi*(*t \** )/*d2 t* = *dfi*(*t \** )/*dt* should be of order *n* � 2, since this factor is multiplied by *h*<sup>2</sup> . The accuracy of the third state derivatives should be of the order *n* � 3, since this factor is multiplied by *h*<sup>3</sup> , etc. In this way, the approximation is correct up to *hn* . *N* is therefore called the *approximation order* of the integration method, see [11].

Many engineering simulation applications require a global relative accuracy of approximately 0.002. If the *local integration* error is of size *el*, then the *per-unit-step integration error* assumes the value of *ep.u.s* = *el/h*. The *global integration error* is proportional to the per-unit-step integration error, as long as the integration error does not accumulate excessively across multiple steps.

In accordance with the previously made observation, this corresponds to an algorithm with an approximation order of *h*<sup>4</sup> for the local integration error. We should require, for example, a local accuracy of 0.0001 [11]. In a digital computer, a real number can only be represented with a finite precision, 32-bit numbers, and double precision for 64-bit numbers. This type of error is called *round-off error*. It occurs in one of the two main general formats that have become common and called *floating point.* The most common problems due to a rounding error occur when multiple steps are rounded in each step.

### **3. An appropriate structured mathematical model of motion in state space**

These practical requirements are determined to using linear models in a process of analyzing general processes [12]. For writing of mathematical models of motion of aircraft in a simulator, we can use state-space description. We define a linear controlled unobserved dynamic system based on the current state, see [13]:

$$\begin{aligned} \dot{\mathbf{x}}(t) &= \mathbf{A}\mathbf{x}(t) + \mathbf{B}\mathbf{u}(t), \\ \mathbf{y}(t) &= \mathbf{C}\mathbf{x}(t), \\ \mathbf{x}(t) &= \mathbf{x}\_0 \end{aligned} \tag{5}$$

*Simulation of a Mathematical Model of an Aircraft Using Parallel Techniques (MPI and GPU) DOI: http://dx.doi.org/10.5772/intechopen.105538*

where *A*, *B*, *C*, *x*, *u,* and *y* have dimensions (*nxn*), (*rxn*), (*lxn*), (*nxl*), (*mxl*), and (*rxl*) are defined matrices, see the description of items [2]. The task easier that we will focus on the control object, the first and second equation from the system (5) can have a general shape:

$$\begin{aligned} \dot{\mathbf{x}} &= \mathbf{A}\mathbf{x} + \mathbf{B}\mathbf{u}, \\ \mathbf{y} &= \mathbf{C}\mathbf{x}. \end{aligned} \tag{6}$$

Eq. (6) in matrix form represents a complex aircraft dynamic system of a mathematical model of aircraft motion comprising in a process of simulation of 11 state variable sensors of information, 18 state variables-coordinates of performing elements in the system. These are two parts and the remaining 38 state variables are divided into variables that represent unmeasured noise and sensor disturbances [9]. The differential equations of mathematical model of aircraft motion have two-parameter control the form [14]:

$$\begin{aligned} \Delta \dot{V} &+ a\_x^V \Delta V + a\_x^a \Delta a + a\_x^\theta \Delta \theta = a\_x^{\delta \mathcal{M}} \Delta \delta\_\mathcal{M} \\ \Delta \theta &+ a\_y^V \Delta V + a\_y^a \Delta a + a\_y^\theta \Delta \theta = a\_y^{\delta \mathcal{M}} \Delta \delta\_V \\ \Delta \dot{\alpha}\_z &+ a\_{mx}^V \Delta V + a\_{mx}^a \Delta a + a\_{mx}^{\dot{\alpha}} \Delta \theta + a\_{mx}^{ax} \Delta \alpha\_z = a\_{mx}^{\delta V} \Delta \delta\_V \\ \Delta \nu &= \alpha\_z, \Delta \nu = \Delta \theta + \Delta a. \end{aligned} \tag{7}$$

Four parts represent a state matrix of mathematical models of motion of aircraft expressed by the state vector *n* = 4. Coefficients *ci* and *ej* represent aerodynamic parameters, their computation, in accordance with (Eq. (7)), we get [14]:

$$
\Delta \dot{V} + a\_x^V \Delta V + a\_x^a \Delta a + a\_x^\theta \Delta \theta + a\_x^H \Delta H = a\_x^{\text{d\ $M}} \Delta \delta\_M
$ 
$$
\Delta \theta + a\_y^V \Delta V + a\_y^a \Delta a + a\_y^\theta \Delta \theta + a\_y^H \Delta H = a\_y^{\text{d\$$
M}} \Delta \delta\_V
\tag{8}
$$

$$
\Delta \dot{H} - \sin \left( \theta\_0 \right) \Delta V - \cos \left( \theta\_0 \right) \Delta V = 0; \Delta \dot{\nu} = a\_x, \Delta \nu = \Delta \theta + \Delta a.
$$

Items Δ*V*, Δ*α*, Δ*θ*, Δ*H*, Δ*δM*, and Δ*δ<sup>V</sup>* are defined in [2]. The following items: *aV x* Δ*V*, *a<sup>α</sup> <sup>x</sup>* Δ*α*, *a<sup>θ</sup> <sup>x</sup>* Δ*θ*, *aH <sup>x</sup>* Δ*H* for flight speed Δ*V*, *a<sup>V</sup> <sup>y</sup>* Δ*V*, *a<sup>α</sup> <sup>y</sup>* Δ*α*, *a<sup>θ</sup> <sup>y</sup>* Δ*θ*, *a<sup>H</sup> <sup>y</sup>* Δ*H* for thumb angle, and other parameters of aircraft motion. Transfer functions have the form [14]:

$$\begin{split} \Delta V(\mathbf{s}) &= -\mathbf{G}\_{V/\delta\mathcal{M}}(\mathbf{s})\Delta\delta\_{\mathcal{M}}(\mathbf{s}) - \mathbf{G}\_{V/\delta V}(\mathbf{s})\Delta\delta\_{V}(\mathbf{s}) \\ \Delta a(\mathbf{s}) &= -\mathbf{G}\_{a/\delta\mathcal{M}}(\mathbf{s})\Delta\delta\_{\mathcal{M}}(\mathbf{s}) - \mathbf{G}\_{a/\delta V}(\mathbf{s})\Delta\delta\_{V}(\mathbf{s}) \\ \mathbf{G}\_{V/\delta\mathcal{M}}(\mathbf{s}) &= a\_{\mathbf{x}}^{\delta\mathcal{M}}(\Delta\_{\mathbf{1}1}\Delta); \mathbf{G}\_{V/\delta V}(\mathbf{s}) = a\_{\mathbf{y}}^{\delta\mathcal{M}}(\Delta\_{\mathbf{21}}\Delta) + a\_{\mathbf{z}\mathbf{x}}^{\delta\mathcal{V}}(\Delta\_{\mathbf{31}}\Delta); \\ \mathbf{G}\_{a/\delta\mathcal{M}}(\mathbf{s}) &= a\_{\mathbf{x}}^{\delta\mathcal{M}}(\Delta\_{\mathbf{12}}\Delta); \mathbf{G}\_{a/\delta V}(\mathbf{s}) = a\_{\mathbf{y}}^{\delta\mathcal{V}}(\Delta\_{\mathbf{22}}\Delta) + a\_{\mathbf{z}\mathbf{x}}^{\delta\mathcal{V}}(\Delta\_{\mathbf{32}}\Delta). \end{split} \tag{9}$$

where Δ*ij* is the subdeterminant of the *i*th row and the *j*th column. Our problem and the proposed solution are described in the following. The first member (Eq. (9)) of the system is employed to describe the change of the aircraft speed Δ*V;* also the equation describing displacement of the speed aircraft depending on the displacement of throttle engine lever *GV/*δ*<sup>M</sup>* (*s*) and the displacement of elevator angle *GV/*δ*<sup>V</sup>* (*s*) is employed. The items of second (Eq. (9)) describe the change of aircraft attack of

angle of the system *Gα/*δ*<sup>M</sup>* (*s*) and *Gα/*δ*<sup>V</sup>* (*s*). The third row of the system (Eq. (9)) defines how to compute these changes of two parameters for a change of aircraft speed, and fourth row of the system (Eq. (9)) defines how to compute these changes of two parameters for a change of an aircraft angle of attack. For the meaning of coefficients *a<sup>δ</sup>Mx*, *aδ<sup>M</sup> <sup>y</sup>*, *aδ<sup>V</sup> <sup>y</sup>*, *aδ<sup>V</sup> mz* and *Δij*(*s*), see the paragraph after (Eq. (20)). In our problem solution, the flight of the aircraft is steady and without any random interferences (wind, storm, or other outer interferences) [12].

#### **3.1 Structure-defined mathematical model of an aircraft motion in a simulator**

At present, matrix notation of a mathematical model in the state space, which we mentioned in the previous section, is very often used. We should express the first equation from (Eq. (6)) next [15]:

$$
\begin{pmatrix}
\dot{\mathbf{x}}\_1 \\
\dot{\mathbf{x}}\_2 \\
\dot{\mathbf{x}}\_3 \\
\dot{\mathbf{x}}\_4
\end{pmatrix} = \begin{pmatrix}
a\_{11} & a\_{12} & a\_{13} & a\_{14} \\
a\_{21} & a\_{22} & a\_{23} & a\_{24} \\
a\_{31} & a\_{32} & a\_{33} & a\_{34} \\
a\_{41} & a\_{42} & a\_{43} & a\_{44}
\end{pmatrix} \begin{pmatrix}
\varkappa\_1 \\
\varkappa\_2 \\
\varkappa\_3 \\
\varkappa\_4
\end{pmatrix}. \tag{10}
$$

The described solution is relatively simple due to description of the analyzed system using equations (Eq. (10)), see [5] for more. Let us divide the mathematical model of aircraft motion into four subsystems. The first will be a subsystem of state variable sensors, the second and third will be subsystems of power elements (two systems), the fourth will measure noise and system faults. Thus, the subsystems have the order *n*<sup>1</sup> = 1, *n*<sup>2</sup> = 1, *n*<sup>3</sup> = 1, *n*<sup>4</sup> = 1. The state space is divided into four parts:

$$\mathbf{x} = (\mathbf{x}\_1, \mathbf{x}\_2, \mathbf{x}\_3, \mathbf{x}\_4)^T = (\mathbf{x}\_1 \mathbf{x}\_1, \mathbf{x}\_1 \mathbf{x}\_2, \mathbf{x}\_1 \mathbf{x}\_3, \dots, \mathbf{x}\_4 \mathbf{x}\_4),\tag{11}$$

where *xixj* represent the items of a state vector. If *i* represents order relevancy *n*, i.e., the number of the subsystem, then *j* – stands for sequential number of the item in the given subsystem. The research architecture of the discussed system matrix *A* in the state space can be formed after multiplying the following shape, see [15]:

$$\begin{aligned} \dot{\mathbf{x}}\_{11} &= a\_{11}\mathbf{x}\_{1}\mathbf{x}\_{1} + a\_{11}\mathbf{x}\_{1}\mathbf{x}\_{2} + a\_{11}\mathbf{x}\_{1}\mathbf{x}\_{3} + a\_{11}\mathbf{x}\_{1}\mathbf{x}\_{4}, \\ \dot{\mathbf{x}}\_{12} &= a\_{12}\mathbf{x}\_{1}\mathbf{x}\_{1} + a\_{12}\mathbf{x}\_{1}\mathbf{x}\_{2} + a\_{12}\mathbf{x}\_{1}\mathbf{x}\_{3} + a\_{12}\mathbf{x}\_{1}\mathbf{x}\_{4}, \\ &\vdots \\ \dot{\mathbf{x}}\_{44} &= a\_{44}\mathbf{x}\_{1}\mathbf{x}\_{1} + a\_{44}\mathbf{x}\_{1}\mathbf{x}\_{2} + a\_{44}\mathbf{x}\_{1}\mathbf{x}\_{3} + a\_{44}\mathbf{x}\_{1}\mathbf{x}\_{4}, \end{aligned} \tag{12}$$

Then we divide 16 subsystems into four parts called isolated subsystems. The above decomposition process is performed by matrix and vector operations. The sixteen blocks *A*11, *A*12, ..., *A*43, *A*<sup>44</sup> of matrix *A* are marked in the order, where:

$$\begin{aligned} \dot{\mathbf{x}}\_{11} &= A\_{11}\mathbf{x}\_{11} + A\_{11}\mathbf{x}\_{12} + A\_{11}\mathbf{x}\_{13} + A\_{11}\mathbf{x}\_{14}, \\ \dot{\mathbf{x}}\_{12} &= A\_{12}\mathbf{x}\_{11} + A\_{12}\mathbf{x}\_{12} + A\_{12}\mathbf{x}\_{13} + A\_{12}\mathbf{x}\_{14}, \\ &\vdots \\ \dot{\mathbf{x}}\_{44} &= A\_{44}\mathbf{x}\_{11} + A\_{44}\mathbf{x}\_{12} + A\_{44}\mathbf{x}\_{13} + A\_{44}\mathbf{x}\_{14}. \end{aligned} \tag{13}$$

*Simulation of a Mathematical Model of an Aircraft Using Parallel Techniques (MPI and GPU) DOI: http://dx.doi.org/10.5772/intechopen.105538*

The mathematical description of isolated subsystems has then the following form:

$$
\dot{\mathbf{x}}\_1 = A\_{11}' \mathbf{x}\_1,\\
\dot{\mathbf{x}}\_2 = A\_{22}' \mathbf{x}\_2,\\
\dot{\mathbf{x}}\_3 = A\_{33}' \mathbf{x}\_3,\\
\dot{\mathbf{x}}\_4 = A\_{44}' \mathbf{x}\_4,\tag{14}
$$

where:

$$A'\_{11} = \begin{pmatrix} A\_{11} \\ A\_{12} \\ A\_{13} \\ A\_{14} \\ A\_{14} \end{pmatrix}, \ A'\_{22} = \begin{pmatrix} A\_{21} \\ A\_{22} \\ A\_{23} \\ A\_{24} \end{pmatrix}, \ A'\_{33} = \begin{pmatrix} A\_{31} \\ A\_{32} \\ A\_{33} \\ A\_{34} \end{pmatrix}, \ A'\_{44} = \begin{pmatrix} A\_{41} \\ A\_{42} \\ A\_{43} \\ A\_{44} \end{pmatrix}. \tag{15}$$

The mutual relations between the first and second isolated subsystems are described by *l*12(*x*) meaning that the equation of the first and second isolated subsystem is:

$$\boldsymbol{I}\_{12}(\mathbf{x}) = \boldsymbol{A}\_{11}^{\prime} \begin{pmatrix} \mathbf{x}\_{1} \\ \mathbf{x}\_{2} \\ \mathbf{0} \\ \mathbf{0} \\ \mathbf{0} \end{pmatrix} = \left( \mathbf{A}\_{11}, \mathbf{A}\_{12}, \mathbf{A}\_{13}, \mathbf{A}\_{14} \right)^{T} \begin{pmatrix} \mathbf{x}\_{1} \\ \mathbf{x}\_{2} \\ \mathbf{0} \\ \mathbf{0} \end{pmatrix}. \tag{16}$$

The state of solution method is known and well described in the reference [5]. The capability of being decomposed can be calculated by means of incidental matrices that can be utilized in cases when the mathematical model of the system is defined. The units in random matrix are positioned as its elements are permuted or transformed so that they can be placed diagonally.
