7 Mathematical modeling – Computational Technologies

Aleksandr G. Churbanov, Alexandr E. Kolesov, and Petr N. Vabishchevich

7 Mathematical modeling

Abstract: Examples of the numerical solution of model problems from various fields of science and engineering are presented in this chapter. All research begins with the mathematical formulation of a problem and the selection of a numerical method. The programs developed here are implemented using the GSL scientific library. The numerical results are visualized employing gnuplot.

7.1 Approximation of experimental data

The problem of constructing an approximate function for a set of one-dimensional data is formulated and solved. The polynomial approximation including the least squares fitting of polynomials is considered.

7.1.1 Problem formulation

Assume that values of a function yi = f(xi), i = 0,1, ..., n are known at points x0 < x1 < ... < xn distributed within an interval [a, b], e.g. we have the measurement data. Let us consider the problem of constructing the function f(x) which approximates this data.

Define on the interval [a, b] a system of functions , and introduce the approximating function

with real coefficients ck, k = 0, 1, ... p (p ≤ n). In the case of polynomial approximation, we have

7.1.2 Least squares method

In equation (7.1), the coefficients ck, k = 0,1,... p can be determined in various ways. Using the least squares method, these coefficients are evaluated by minimizing the sum of the squared differences between the approximating function and the experimental values:

Define a p-dimensional vector of unknown coefficients c = {ck} = {c1, c2, ..., cp} and let f = {Φ(xi)} = {Φ(x1), Φ(X2),..., Φ(xn)}. According to equation (7.1), we get

f = Rc,

where R is the rectangular matrix R = {rik} with the elements rik = Φk(xi), i = 0, 1,..., n, k = 0, 1, ..., p.

7.1.3 Program

For a polynomial approximation based on the least squares method, we apply the function gsl_multi fit_linear from GSL. The code implementing this algorithm is as follows:

Listing 7.1.

In this example, the values of the function x sin(x) at 21 points on the interval [0, 4] are treated as the observation data for approximation.

7.1.4 Computations

The results of computations (the values of the approximating function at observation points) are written in the file mnk. dat. The gnuplot script (file mnk. gnu) for visualization is also given.

This program produces the following coefficients of the approximating polynomial:

The input data and the resulting function are plotted in Figure 7.1.

Fig. 7.1. Third-degree polynomial approximation.

The second-degree polynomial approximation yields the coefficients

The corresponding plot is depicted in Figure 7.2.

Fig. 7.2. Second-degree polynomial approximation.

7.2 Predator-prey system

A numerical study of the predator-prey system based on the numerical solution of the Lotka-Volterra equations is presented below.

7.2.1 Mathematical model

We consider two species of animals that inhabit a certain area. Rabbits (prey) feed on plants, and foxes (predators) hunt for rabbits. Let v be the number of rabbits and w the number of foxes. The behavior of the system is governed by the Lotka-Volterra equations:

where t is time.

The number of prey per unit of time is increased due to the birth of new prey (the reproduction rate by the number of prey), and it is decreased due to prey being eaten. Accordingly, in equation (7.2), the coefficient α1 is the growth rate of the prey population (birthrate), and α2 stands for the coefficient of predation on a prey (the probability that a prey meets a predator and will be eaten). The growth of the predator population per unit of time is proportional to the quality of the food, and the loss is due to natural death. In equation (7.3), α3 is the coefficient of predator mortality, and α4 ensures the growth of the population by eating prey.

The dynamics of the system (7.2), (7.3) is determined by the initial conditions


7.2.2 Stationary state

We can define the stationary state of the system (7.2), (7.3). The conditions

are fulfilled if

It seems reasonable that in some range of the initial conditions the number of species will become stable. However, a significantly more complex unsteady behavior of the above-mentioned nonlinear dynamic system is also observed.

7.2.3 Numerical method

For the numerical study of the mathematical model of predator-prey, the Runge−Kutta method with an adaptive step (Runge−Kutta−Fehlberg method) from GSL is applied.

In numerically solving the Cauchy problem for the differential equations, the Runge−Kutta method with the evaluation of an integration step is used. The simplest strategy of changing an integration step is based on solving the problem with two time steps. The second one is more interesting and involves the solution of the problem with a single step, but using methods with different accuracy.

The Cauchy problem for the system of ODEs under the consideration appears to be

In vector notation, this problem may be rewritten as the Cauchy problem for the single equation

where u = {u1, u2, ..., um} is the vector of unknowns, and f = {f1, f2, ..., fm} is the right-hand side vector.

In one-step Runge−Kutta methods, the transition from the time level tn to the next level tn+1 = tn + τ may be written as


The computational implementation is based on s computations of the function f (the s-stage algorithm).

In the explicit Runge-Kutta-Fehlberg method, we have

The approximate solution is determined with fourth-order accuracy according to

To control accuracy, we employ the numerical solution , which is evaluated from

and corresponds to the method of fifth-order accuracy.

7.2.4 Program

To solve numerically the Cauchy problem for the system (7.2)-(7.4), the Runge−Kutta− Fehlberg method of fourth-fifth order accuracy mentioned above is employed. This method is implemented by means of GSL and is presented in the following code:

Listing 7.2.

We prescribe the initial time step; the program stops when the required relative or absolute accuracy is achieved.

7.2.5 Specification of the right-hand side

Here we also present the part of the program corresponding to the specification of the right-hand side for the system of equations (7.2), (7.3).

Listing 7.3.

Some predictions for the predator-prey system visualized using gnuplot are given in the following.

7.2.6 Visualization script

The numerical results are written in the file lotka. dat for further visualization via gnuplot. To use this visualization system, we construct the script written in the file lotka. gnu; it is shown as follows:

The plot obtained after data processing is saved in the file lotka. png.

7.2.7 The basic variant of input data

Assume that the coefficients of the system of equations (7.2), (7.3) are the following:

α1 = 0.05, α2 = 0.005, α3 = 0.01, α4 = 0.0002.

In this case, the stationary state is

, .

The numerical solution of the problem (7.2), (7.3) with the initial conditions v0 = 10 and w0 = 10 is shown in Figure 7.3.

Fig. 7.3. Population dynamics for the predator-prey system.

7.2.8 Variation of parameters

Using the program presented above, we can study the behavior of the system varying parameters of the problem. In particular, it is interesting to analyze the influence of the initial conditions v0, w0 on the dynamics of the population.

For example, Figure 7.4 demonstrates the time-evolution of the number of species for the initial conditions v0 = 100 and w0 = 10.

Fig. 7.4. Population dynamics for v0 = 100 and w0 = 10.

To continue this parametric study, we can also vary the coefficients α1, α2, α3, α3.

7.3 Water infiltration

In this section we consider a mathematical model describing water infiltration from the Earth’s surface into the soil. It is based on a one-dimensional quasi-linear parabolic equation. The computational algorithm for studying this problem is presented below along with the numerical results.

7.3.1 The equation of water transport

The water precipitation falls on the Earth’s surface and seeps into the soil within a certain depth. Denote by w the water content (soil moisture) and let w* be the maximum water content (0 ≤ w ≤ w*). To investigate the water infiltration in an unsaturated zone, we apply the one-dimensional equation of water transport. Assume that the axis z is directed into the ground, then the equation of water transport may be written as

The coefficient of hydraulic conductivity kw depends essentially on the water content, i.e. we have kw = kw(w). For example, the following relation (S. F. Averyanov’s formula) is typical:


where k is the flow coefficient and ϭ ranges from 2 to 4 or more. In equation (7.5), we have

where ψ is the suction height. The dependence


is widely used in predictions of problems like this. The coefficients α and β in (7.7) are determined experimentally.

Let u = w/w* be the relative water content; then, in view of (7.6), (7.7), equation (7.5) takes the form

7.3.2 Self-similar solution

For the water infiltration process, the self-similar solution is the exact solution of equation (7.8) corresponding at the asymptotic stage to the motion of the wetting boundary with some constant speed v, i.e. it is a traveling wave-like self-similar solution.

Therefore, the function ζ = zvt is selected as the self-similar variable. Some self-similar solutions of this type can be obtained in the explicit form. For instance, for ϭ = 2, 3, 5, we have

with v = k/w*. This set of solutions allows us to study the influence of the parameter ϭ on a steady water profile.

Visualization of the self-similar solutions presented above is performed for const = 0. The corresponding script (file self. gnu) is as follows:

The self-similar solutions for ϭ = 2, 3, 5 and βw*, v = 1 are depicted at time moments t = 5/3, t = 10/3 and t = 5 in Figure 7.5.

Now we shall discuss the computational algorithm for studying this problem and present some numerical results.

7.3.3 Difference scheme

To solve numerically the boundary value problem for water transport, we apply the finite-difference scheme with the coefficients taken from the previous time level. The numerical implementation is based on the solution of a system of linear equations with a tridiagonal matrix.

Fig. 7.5. Water profile dynamics for ϭ = 2 (black line), ϭ = 3 (gray) and σ = 5 (silver).

We seek the solution of the water infiltration problem where the soil is initially unsaturated. We restrict ourselves to the case with = 1, k/w* = 1. The equation (7.8) is considered for 0 < z < l with the following boundary and initial conditions:

For the numerical solution of the boundary value problem, we introduce a uniform mesh on the segment [0, l]:

where w is the set of interior nodes (i = 1, 2, ... , N − 1). Denote the approximate solution at the time moment tn = by yn, where τ is the time step. We employ the linearized difference scheme with coefficients taken from the previous time level. For the interior nodes of the spatial grid and the time moments n = 0, 1, ... , the difference scheme has the form

The boundary conditions appear as

7.3.4 Program

In the following we present the code for solving the boundary value problem for the equation of water transport. A fine enough mesh in space (N = 100) is used in calculations conducted with the coefficient ϭ = 2.

Listing 7.4.

As above, the visualization of the predictions is performed by gnuplot.

7.3.5 Visualization

The results at the specified time moments (at the interval of nWr = 50 steps) are written in the file bvpt . dat. The script for gnuplot is contained in the file bvpt . gnu and has the following form:

The plots of the solution at the individual time moments (multiples of δ t = 0.5) are written in file bvpt . png.

7.3.6 Predictions

The dynamics of the infiltration process for the basic case is shown in Figure 7.6. We can observe that during the time-evolution process the infiltration front transforms practically to the self-similar solution (see Figure 7.5 for a comparison).

Fig. 7.6. Water dynamics with the constant soaking on the Earth’s surface (u(0, t) = 1).

A more complicated situation, when soaking on the Earth’s surface is varying in time, is reflected in Figure 7.7. In this case, until the time moment 0.5T, the soaking exists and is constant, and after this point the soaking is absent. This situation is modeled by specifying the boundary condition at z = 0 in the form

Fig. 7.7. The case with soaking varying in time.

The influence of the parameter ϭ is demonstrated in Figure 7.8. As for the self-similar solutions, the increase of ϭ results in sharping the soaking front (increasing its gradient).

Fig. 7.8. Water dynamics for ϭ = 3.

7.4 Torsion of cylindrical bars

The torsion of cylindrical bars is the classical problem of the theory of elasticity. In this section, we present a computational algorithm for solving this problem and discuss the results obtained.

7.4.1 Problem formulation

The problems of the elastic torsion of cylindrical bars belong to the class of the basic and well-studied problems of solid mechanics. To study the torsion of bars with simply connected cross-sections of an arbitrary shape, it is necessary to solve the Dirichlet problem for the two-dimensional Poisson equation. These boundary value problems present the class of the basic problems of mathematical physics, which are provided with the well-known numerical methods.

We consider a cylindrical bar with a cross-section D in a plane x1, x2 (the axis x3 is directed along the axis of the bar) subjected to uniform torsion. Let α be the swirl angle per unit length of the bar. For the components of the stress tensor, we have

where u = u(x), x = (x1, x2). For elastic deformation, the stress function u(x) satisfies the equation


where G is the shear modulus of bar materiel. From (7.9), (7.10), we arrive at the Poisson equation

Assume that the section D of the bar is simply connected. Then equation (7.11) is complemented by the homogeneous Dirichlet condition on the boundary:


In processing the computational results, the focus is on the torque

Multiplying equation (7.9) by u(x) and integrating the resulting equation over the domain Ω, we get

where the function

defines local features of the loads.

7.4.2 Difference problem

To solve numerically the Dirichlet problem (7.11), (7.12), we will consider the problem in a rectangle completely containing the domain D (). Thus, we have

In , we introduce a uniform mesh with steps h1 and h2 in the corresponding direction, respectively. Let be the set of interior nodes, i.e.

and is the set of boundary nodes. The difference solution of the problem (7.11), (7.12) is denoted by y(x), x ∈ .

For the nodes outside the computational domain, in view of the boundary condition (7.12), we assume


At other nodes of the spatial grid, equation (7.11) is approximated by the following relations:


The difference scheme (7.13), (7.14) results from the piecewise linear approximation of the boundary of the computational domain D with grid nodes located on the boundary.

7.4.3 Numerical algorithm and program

To find the difference solution of the problem, the iterative successive over-relaxation (SOR) method is used.

Let yn(x) be the approximate solution at the n-th iteration and 0 < σ < 2 be a relaxation parameter; then

The program for solving the torsion problem is as follows:

Listing 7.5.

In the particular case presented below, the torsion problem is studied for the bar with the rotated L-shape cross-section, i.e. the cross-section is the rectangle (l1 = 2, l2 = 1) without the rectangular part

The main calculation parameters are the following. The uniform mesh with 201 x 101 nodes is used. The relaxation parameter σ is equal to 1.9. The iterative process stops as soon as the maximum residual becomes less then 1 x 10-7.

7.4.4 Numerical results

The numerical solution is written to the file tors . dat. In addition, the function w(x) is also calculated at the interior nodes of the mesh (the file tors1. dat). For visualization, the script (the file tors . gnu) is employed.

During the run of the program, the number of iterations and tolerance are printed:

The predicted solution is shown in Figure 7.9 as colored values of the stress function u(x) with the legend at the right. Figure 7.10 presents the calculated function w(x) in the same style.

Using this program, it is possible to investigate the influence of the mesh size N1, N2 on the accuracy of the numerical solution as well as the impact of the relaxation parameter σ on the convergence rate. We might also consider the torsion problem for bars with other cross-sections.

Fig. 7.9. Values of the stress function u(x).

Fig. 7.10. Values of the function w(x).