Alexei Yurievich Vinogradov
Numerical methods of solving stiff and non-stiff
boundary value problems
Propositions: Improvement of S.K.Godunov’s method of orthogonal sweep, 3 methods for nonstiff cases of boundary value problems, 2 methods for stiff cases of boundary value problems, 1
method for calculating composite shells and with frames, a C++ program for the best method
proposed.
Monograph
2019 Moscow, Russia
AlexeiVinogradov@yandex.ru,
+7(963)991-05-10, +7(977)810-55-23 (WhatsApp, Viber)
2
Table of contents
Table of contents.
2
Introduction.
4
Chapter 1. Known formulas of the theory of matrices for ordinary
differential equations.
10
Chapter 2. Improvement of S.K.Godunov’s method of orthogonal sweep
for solving boundary value problems with stiff ordinary differential
equations.
12
2.1. The formula for the beginning of the calculation by
S.K.Godunov’s sweep method.
12
2.2. The second algorithm for the beginning of the calculation by
S.K.Godunov’s sweep method.
16
2.3. The replacement of the Runge-Kutta’s numerical integration
method in S.K.Godunov’s sweep method.
17
2.4 Matrix-block realizations of algorithms for starting calculation
by S.K.Godunov’s sweep method.
2.5.
Conjugation of parts of the
17
integration interval
for
S.K.Godunov’s sweep method.
20
2.6. Properties of the transfer of boundary value conditions in
S.K.Godunov’s sweep method.
2.7. Modification of S.K.Godunov’s sweep method.
22
23
Chapter 3. The method of "transferring of boundary value conditions"
(the direct version of the method) for solving boundary value problems
with non-stiff ordinary differential equations.
25
Chapter 4. The method of "additional boundary value conditions" for
solving boundary value problems with non-stiff ordinary differential
equations.
26
Chapter 5. The method of "half of the constants" for solving boundary
value problems with non-stiff ordinary differential equations.
29
Chapter 6. The method of "transferring of boundary conditions" (stepby-step version of the method) for solving boundary value problems
with stiff ordinary differential equations.
6.1. The method of "transfer of boundary value conditions" to any
31
31
3
point of the interval of integration.
6.2. The case of "stiff" differential equations.
33
6.3. Formulas for computing the vector of a particular solution of
inhomogeneous system of differential equations.
6.4. Applicable formulas for orthonormalization.
35
39
Chapter 7. The simplest method for solving boundary value problems
with stiff ordinary differential equations without orthonormalization the method of "conjugation of sections of the integration interval",
which are expressed by matrix exponents.
41
Chapter 8. Calculation of shells of composite and with frames by the
simplest method of "conjugation of sections of the integration interval".
8.1. The variant of recording of the method for solving stiff
boundary value problems without orthonormalization - the method of
"conjugation of sections, expressed by matrix exponents "- with positive
directions of matrix formulas of integration of differential equations.
43
8.2. Composite shells of rotation.
44
8.3. Frame, expressed not by differential, but algebraic equations.
47
8.4. The case where the equations (of shells and frames) are
expressed not with abstract vectors, but with vectors, consisting of
specific physical parameters.
51
Appendix. Computational experiments (a C++ program).
55
List of published works.
64
4
Introduction.
Relevance of the problem:
The solution to the problem of weight reduction of structures is associated with their
complication and the use of thin-walled elements. Even the simplest variant method of
constructive optimization requires parametric studies on a computer using numerical methods for
solving boundary value problems. The most famous among them are:
- finite-difference methods for constructing approximate solutions of differential equations
using finite-difference approximations of derivatives;
- various modifications of the finite element method, the Bubnov-Galerkin method, the
Rayleigh-Ritz method, which are based on approximations of solutions of differential equations
by finite linear combinations of given functions: polynomials, trigonometric functions, etc .;
- methods for the numerical determination of the integrals of the ordinary differential
equations Runge-Kutta, Volterra, Picard, etc.
The main success of the methods of finite differences and finite elements is that on their
basis universal algorithms are constructed and packages of applied programs for calculating
complex structures are created. The constructed computing facilities are able to detect the flow
of forces in the structure and, therefore, the most strained elements of it. Nevertheless, they
require unjustifiably high costs of the programmer's efforts and powerful computing tools when
the task is to determine the stresses in the places of their concentration.
The most obvious effectiveness of the methods of numerical integration of ordinary
differential equations consists in calculating individual parts of complex spatial structures and
their individual thin-walled elements with the refinement of the stress-strain state in the places of
its rapid change. Efficiency is determined by the small expenses of the programmer's work, by
the small expenditure of computer time and the computer's operational memory.
5
Thus, increasing the effectiveness of known numerical methods, constructing their
modifications and constructing new methods, is an urgent research task.
The proposed scientific novelty consists in the following:
1.
S.K.Godunov’s method of orthogonal sweep has been improved,
2.
A method of "transferring of boundary value conditions" (a direct version of the
method) is proposed for solving boundary value problems with non-stiff ordinary
differential equations,
3.
A method of "additional boundary value conditions" is proposed for solving boundary
value problems with non-stiff ordinary differential equations,
4.
A "half of constants" method is proposed for solving boundary value problems with
non-stiff ordinary differential equations,
5.
A method of "transferring of boundary value conditions" (step-by-step version of the
method) is proposed for solving boundary value problems with stiff ordinary
differential equations,
6.
The simplest method for solving boundary value problems with stiff ordinary
differential equations without orthonormalization is proposed - the method of
"conjugation of sections of the integration interval", which are expressed by matrix
exponentials,
7.
The simplest method for calculating the shells of composite and with frames is
proposed.
6
Some of the works on which the methods are based are published jointly with Dr.Sc.
Professor Yu.I.Vinogradov.
Contribution of Dr.Sc. Professor Yu.I. Vinogradov in these joint publications consisted
either of 1) in the discussion of the results of verification calculations of those formulas and
methods proposed by A.Yu. Vinogradov, or that 2) in addition to A.Yu.Vinogradov’s methods
Yu.I.Vinogradov proposed statement that the Cauchy matrix can be computed not only in the
form of matrix exponentials, but in addition there is the possibility of calculating them in the
sense of Cauchy-Krylov functions, using for this purpose the analytical solutions of the systems
of differential equations of the structural mechanics of plates and shells obtained by someone,
that 3) Yu.I. Vinogradov proposed his own, different from A.Yu.Vinogradov’s formula, the
formula for computing the vector of the particular solution of an inhomogeneous system of
ordinary differential equations, which, however, looks much more complicated than the simple
A.Yu.Vinogradov’s formula.
Also, in the co-authors of some articles, Yu.A.Gusev and Yu.I.Klyuyev. Their
contribution to the publication material consisted in performing multivariate verification
calculations in accordance with the formulas, algorithms and methods proposed by A.Yu.
Vinogradov in his Ph.D. thesis. The Ph.D. thesis was defended in 1996.
In addition, we can say that on the basis of the material from A.Yu.Vinogradov’s Ph.D.
thesis completed two more candidate's physical and mathematical dissertations under the
direction of Yu.I.Vinogradov, whose material consists mainly of a multivariate application and
verification by calculations of what was proposed by A.Yu.Vinogradov in his Ph.D. thesis, in
application to various concrete problems of construction mechanics of thin-walled shells with the
identification and analysis of properties of formulas, algorithms and methods from
A.Yu.Vinogradov’s Ph.D. thesis.
Here are the data of these two dissertations:
Year: 2008 Petrov, Vitaliy Igorevich "Reduction of boundary value problems to initial
7
problems and investigation of stress concentration in thin-walled constructions by the
multiplicative method"
Scientific degree: Candidate of Physical and Mathematical Sciences
Specialty code of VAC: 05.13.18 Specialty: Mathematical modeling, numerical methods
and program complexes.
Year: 2003 Gusev, Yuri Alekseevich "Multiplicative algorithms for the transfer of
boundary conditions in problems of the mechanics of shell deformation"
Scientific degree: Candidate of Physical and Mathematical Sciences
Specialty code of VAC: 01.02.04 Specialty: Mechanics of a deformable solid.
In addition, in accordance with the modern capabilities of the Internet and in the presence
of new theses in the open access it was revealed the use of materials from A.Yu.Vinogradov
Ph.D. thesis with relevant references to A.Yu.Vinogradov’s relevant articles in the following
candidate and doctoral dissertations of technical and physical and mathematical sciences:
Year: 2005 Russian academy of sciences, Siberian Branch, Institute of Computational
Technologies, Novosibirsk, Yurchenko Andrey Vasilevich “Numerical solution of
boundary problems of elastic deformation of composite shells of rotation”
Scientific degree: Candidate of Physical and Mathematical Sciences
Specialty code of VAC: 05.13.18 - mathematical modeling, numerical methods and
program complexes.
Year: 2003 Kochetov, Sergey Nikolaevich, Moscow "Methods and algorithms for
determining the stress-strain state of thin-walled reinforced rotation structures from
nonlinear elastic material"
Scientific degree: Candidate of Physical and Mathematical Sciences
Specialty code of VAC: 05.23.17 - construction mechanics.
Year: 1998 Chekanin, Alexander Vasilievich, Moscow "Development of the method of
superelements applied to the problems of statics and dynamics of thin-walled spatial
systems"
Scientific degree: Doctor of Technical Sciences
Specialty code of VAC: 05.23.17 - construction mechanics.
8
Year: 2005 Golushko, Sergey Kuzmich, Novosibirsk "Direct and inverse problems of the
mechanics of elastic composite plates and shells of revolution"
Scientific degree: Doctor of Physics and Mathematics
Specialty code of VAC: 01.02.04 - mechanics of deformable solids.
Year: 2003 Gazizov, Khatib Sharifzyanovich, Ufa "Development of theory and methods
for calculating the dynamics, stiffness and stability of composite shells of revolution"
Scientific degree: Doctor of Technical Sciences
Specialty code of VAC: 01.02.04 - mechanics of deformable solids.
Year:
2001 Shlenov,
Alexey Yurievich, Moscow "Dynamics of structurally
inhomogeneous shell structures with allowance for the elastic-plastic properties of the
material"
Scientific degree: Candidate of physical and mathematical sciences
Specialty code of VAC: 01.02.04 - mechanics of deformable of a solid.
Year: 1996 Rogov, Anatoly Alekseevich, Moscow "Dynamics of the pipeline after the
break"
Scientific degree: Candidate of Physical and Mathematical Sciences
Specialty code of VAC: 01.02.04 - mechanics of deformable solids.
A.Yu.Vinogradov’s articles were published in VAC magazines:
-
Reports of the Academy of Sciences of the Russian Federation - 2 articles
-
Mechanics of the Solid Body of the Russian Academy of Sciences - 2 articles
-
Journal of Computational Mathematics and Mathematical Physics, Russian Academy of
Sciences - 1 article
-
Mathematical Modeling of the Russian Academy of Sciences - 2 articles
-
Fundamental research - 1 article
-
Modern problems of science and education - 1 article
-
Modern high technology - 1 article.
9
In general, the formulas for solving boundary value problems for linear ordinary
differential equations borrowed in this paper were taken from only 4 sources:
1. Abramov A.A. On the transfer of boundary conditions for systems of linear differential
equations (variant of the sweep method) // Journal of Computational Mathematics and
Mathematical Physics, 1961. - T.I. -N3. -C.542-545.
2. Berezin I.S., Zhidkov N.P. Methods of calculation. M.: Fizmatgiz, 1962.-T.2. - 635 p.
3. Gantmakher F.R. Theory of matrices. – М.: Science, 1988. - 548 p,
4. Godunov S.K. On the numerical solution of boundary value problems for systems of
linear ordinary differential equations. Uspekhi Matematicheskikh Nauk, 1961. -T. 16, no.
3, (99). - p.171-174.
10
Chapter 1. Known formulas of the theory of matrices for ordinary differential equations.
All methods are given by the example of a system of differential equations of the
cylindrical shell of a rocket - a system of ordinary differential equations of the 8th order (after
the separation of partial derivatives by Fourier’s method).
The system of linear ordinary differential equations has the form:
Y ( x) AY ( x) F ( x) ,
where Y (x) – the required vector-valued function of the dimension 8х1, Y (x) – the derivative
of the required vector-valued function of the dimension 8х1, A – square matrix of coefficients
of a differential equation of dimension 8х8, F (x) – vector-function of external action on the
system of dimension 8х1.
Hereinafter, vectors are denoted by boldface instead of dashes over letters.
The boundary value conditions have the form:
UY (0) u,
VY (1) v ,
where Y (0) – the value of the required vector-function on the left-hand side x = 0 of dimension
8x1, U – rectangular horizontal matrix of coefficients of boundary conditions of the left edge of
dimension 4х8, u – vector of external influences on the left edge of dimension 4х1,
Y (1) – the value of the required vector-function on the right-hand side x = 1 of dimension
8x1, V – rectangular horizontal matrix of coefficients of the boundary conditions of the right
edge of dimension 4х8, v – vector of external actions on the right edge of dimension 4х1.
In the case when the system of differential equations has a matrix with constant
coefficients A = const, the solution of Cauchy’s problem has the form [Gantmakher]:
x
Y ( x) e A( x x 0)Y ( x0 ) e Ax e At F (t )dt ,
x0
11
where e A( x x 0) E A( x x0 ) A 2 ( x x0 ) 2 / 2! A3 ( x x0 ) 3 / 3!... , where E - is the unit
matrix.
The matrix exponent can still be called Cauchy’s matrix and can be written as:
K ( x x0 ) K ( x x0 ) e A( x x 0 ) .
Then the solution of Cauchy’s problem can be written in the form:
Y ( x) K ( x x0 )Y ( x0 ) Y ( x x0 ) ,
x
where Y ( x x0 ) e Ax e At F (t )dt is the vector of a particular solution of an inhomogeneous
x0
system of differential equations.
From the theory of matrices [Gantmakher], the property of multiplication of matrix
exponentials (Cauchy’s matrices) is known:
K ( xi x0 ) K ( xi xi 1 ) K ( xi 1 xi 2 ) ... K ( x2 x1 ) K ( x1 x0 )
In the case when the system of differential equations has a matrix with variable
coefficients A A(x ) , the solution of Cauchy’s problem is proposed, as is known, to be sought
using the property of multiplication of Cauchy’s matrices. That is, the interval of integration is
divided into small sections and in small parts Cauchy’s matrix approximately calculated by the
formula for the constant matrix in the exponent. And then Cauchy’s matrices calculated on small
sections are multiplied:
K ( xi x0 ) K ( xi xi 1 ) K ( xi 1 xi 2 ) ... K ( x2 x1 ) K ( x1 x0 ) ,
where Cauchy’s matrices are approximately computed by the formula:
K ( xi 1 xi ) e A( xi ) xi exp( A( xi ) xi ) , где xi xi 1 xi .
12
Chapter 2. Improvement of S.K.Godunov’s method of orthogonal sweep for solving
boundary value problems with stiff ordinary differential equations.
2.1. The formula for the beginning of the calculation by S.K.Godunov’s sweep method.
Let us consider S.K.Godunov’s sweep method problem.
Suppose that we consider the shell of the rocket. This is a thin-walled tube. Then the
system of linear ordinary differential equations will be of the 8th order, the matrix A of
coefficients will have the dimension 8x8, the required vector-function Y (x ) will have the
dimension 8x1, and the matrices of the boundary conditions will be rectangular horizontal
dimensions 4x8.
Then in S.K.Godunov’s method for such a problem the solution is sought in the following
form [Godunov]:
Y ( x) Y1 ( x)c1 Y2 ( x)c 2 Y3 ( x)c3 Y4 ( x)c 4 Y ( x) ,
or it can be written in the matrix form:
Y ( x) Yматрица ( x)c Y ( x) ,
where vectors Y1 ( x), Y2 ( x), Y3 ( x), Y4 ( x) are linearly independent vector-solutions of the
homogeneous system of differential equations, and the vector Y (x) is a vector of a particular
solution of the inhomogeneous system of differential equations.
Here Yматрица ( x) Y1 ( x), Y2 ( x), Y3 ( x), Y4 ( x) is the matrix of dimension 8x4, and c is
the corresponding vector of dimension 4x1 with the required constants c1 , c2 , c3 , c4 .
But in general, the solution for such a boundary-value problem with dimension 8 (outside
the framework of S.K.Godunov's method) can consist not of 4 linearly independent vectors
Y1 ( x), Y2 ( x), Y3 ( x), Y4 ( x) , but entirely of all 8 linearly independent solution vectors of the
homogeneous system of differential equations:
13
Y ( x) Y1 ( x)c1 Y2 ( x)c2 Y3 ( x)c3 Y4 ( x)c4
Y5 ( x)c5 Y6 ( x)c6 Y7 ( x)c7 Y8 ( x)c8 Y ( x).
And just the difficulty and problem of S.K.Godunov’s method is that the solution is
sought with only half the possible vectors and constants, and the problem is that such a solution
with half the constants must satisfy the conditions on the left edge (the starting edge for the
sweep) for all possible values of the constants, in order to find these constants from the
conditions on the right edge.
That is, in S.K.Godunov’s method, there is a problem of finding such initial values
Y1 (0), Y2 (0), Y3 (0), Y4 (0), Y (0) of the vectors Y1 ( x), Y2 ( x), Y3 ( x), Y4 ( x), Y ( x) , so that you can
start the run from the left edge x = 0, that is, that the conditions UY (0) u on the left edge are
satisfied for any values of the constants c1 , c2 , c3 , c4 .
Usually this difficulty is "overcome" by the fact that differential equations are written not
through functionals, but through physical parameters and consider the simplest conditions on the
simplest physical parameters so that the initial values Y1 (0), Y2 (0), Y3 (0), Y4 (0), Y (0) can be
guessed. That is, problems with complex boundary conditions can not be solved in this way: for
example, problems with elastic conditions at the edges.
Below we propose a formula for the initiation of computations by S.K.Godunov’s
method.
We perform the line orthonormalization of the matrix equation of the boundary
conditions on the left edge:
UY (0) u ,
where the matrix U is rectangular and horizontal dimension 4x8.
As a result, we obtain an equivalent equation of boundary conditions on the left edge, but
already with a rectangular horizontal matrix U орто of dimension 4x8, which will have 4
orthonormal rows:
14
U ортоY (0) uорто ,
where, as a result of orthonormalization of the matrix U , the vector u is transformed into the
vector uорто .
How to perform line orthonormation of systems of linear algebraic equations can be
found in [Berezin, Zhidkov].
We complete the rectangular horizontal matrix U орто to a square non-degenerate matrix
W:
W
U орто
M
,
where a matrix M of dimension 4х8 must complete the matrix U орто to a non-degenerate square
matrix W of dimension 8х8.
As matrix M rows, we can take those boundary conditions, that is, expressions of those
physical parameters that do not enter the parameters of the left edge or are linearly independent
with them. This is quite possible, since for boundary value problems there are as many linearly
independent physical parameters as the dimensionality of the problem, that is, in this case there
are 8 of them, and if 4 are given on the left edge, then 4 can be taken from the right edge.
We complete the orthonormalization of the constructed matrix W , that is, we perform the
line orthonormalization and obtain a matrix Wорто of dimension 8x8 with orthonormal rows:
Wорто
U орто
M орто
.
We can write down that
T
Yматрица (0) (M орто ) транспонированная M орто
.
Then, substituting in the formula of S.K. Godunov’s method, we get:
Y (0) Yматрица (0)c Y (0)
15
or
T
Y (0) M орто
c Y (0) .
We substitute this last formula into the boundary conditions U ортоY (0) uорто and
obtain:
T
U орто[M орто
c Y (0)] uорто .
From this, we obtain that on the left-hand side the constants c no longer influence
T
anything, since U ортоM орто
0 and it remains only to find the vector Y (0) from the
expression:
U ортоY (0) uорто .
But the matrix U орто has a dimension of 4x8 and it must be supplemented to a square
non-degenerate one in order to find the vector Y (0) from the solution of the corresponding
system of linear algebraic equations:
U орто
M орто
Y (0)
uорто
,
0
where 0 is any vector, including a vector of zeros.
Hence we obtain by means of the inverse matrix:
Y (0)
U орто
1
uорто
M орто
0
.
Then the formula for starting the computation by S.K. Godunov's method is as follows:
Y (0) M
T
орто
c
U орто
M орто
1
uорто
0
.
From the theory of matrices [Gantmakher] it is known that if the matrix is orthonormal,
then its inverse matrix is its transposed matrix. Then the last formula takes the form:
16
Y ( 0) M
T
орто
c
T
T
Y (0) M орто
c U орто
U орто
T
uорто
M орто
T
M орто
,
0
uорто
0
,
T
T
T
Y (0) M орто
c U орто
uорто M орто
0,
T
T
Y (0) M орто
c U орто
uорто .
T
T
The column vectors of the matrix M орто
and the vertical convolution vector U орто
uорто
are linearly independent and satisfy the boundary condition UY (0) u .
2.2. The second algorithm for the beginning of the calculation by S.K.Godunov’s sweep
method.
This algorithm requires the addition of a matrix of boundary conditions U to a square
non-degenerate one:
U
.
M
The initial values Y1 (0), Y2 (0), Y3 (0), Y4 (0), Y (0) are found from the solution of the
following systems of linear algebraic equations:
U
u
Y (0)
,
M
0
1
0
U
0
Yi (0) , где i ,
0
M
i
0
where 0 is a vector of zeros of dimension 4х1.
0
1
,
0
0
0
0
,
1
0
0
0
,
0
1
17
The column vectors Yi (0) and the column vector Y (0) are linearly independent and,
taking part in the formation of the vector Y (0) , satisfy the boundary condition UY (0) u .
2.3.
The replacement
of the Runge-Kutta’s numerical integration
method
in
S.K.Godunov’s sweep method.
In S.K.Godunov's method, as shown above, the solution is sought in the form:
Y ( x) Yматрица ( x)c Y ( x) .
At each specific section of S.K.Godunov's method of sweeping between points of
orthogonalization one can use the theory of matrices instead of Runge-Kutta’s method and
perform the calculation through Cauchy’s matrix:
Y матрица ( x j ) K ( x j xi )Y матрица ( xi ) .
So perform calculations faster, especially for differential equations with constant
coefficients, since in the case of constant coefficients it is sufficient to calculate once Cauchy’s
matrix in a small section and then only multiply by this once computed Cauchy’s matrix.
Similarly, through the theory of matrices, we can also calculate the vector Y (x) of a
particular solution of an inhomogeneous system of differential equations. Or, for this vector,
Runge-Kutta’s method can be used separately, that is, one can combine the theory of matrices
and Runge-Kutta’s method.
2.4 Matrix-block realizations of algorithms for starting calculation by S.K.Godunov’s
sweep method.
We consider a system of linear algebraic equations expressing the boundary conditions
for x=0
UY (0) u
18
Suppose that there is a constructed quadratic nondegenerate matrix G
Similarly, we write a vector g
U
.
U*
u
where the introduced vector u * is unknown.
*
u
We write the system of linear algebraic equations
GY (0) g
or in block form
U
u
Y (0) * .
*
U
u
It follows that
U
Y (0) *
U
1
u
G 1 g Ng .
*
u
Imagine N T T * .
Then
Y (0)
U
U*
1
u
G 1 g Ng T
*
u
T*
u
Tu T * u * .
*
u
At the same time, we remember that the solution of the boundary value problem is sought
in the form
Y ( x) Yматрица ( x)c Y ( x)
.
Comparing
Y (0) T * u* Tu and
Y (0) Yматрица (0)c Y (0)
given that here the vector of unknown constants is u* c , we obtain the initial values of the
vectors for the beginning of integration in S.K.Godunov’s method.:
Yматрица (0) T и
Y (0) Tu .
Another matrix derivation can be stated in the following form.
19
We transform the system
U
u
Y (0) *
*
U
u
by line orthonormalization to an equivalent system with orthonormal rows
W
w
Y (0) * .
*
W
w
Then we can write
Y (0)
W
W*
1
w
W
*
w
W*
T
w
WT
*
w
W *T
w
W T w W *T w * .
*
w
Making a comparison of two expressions:
Y (0) Yматрица (0)c Y (0)
Y (0) W *T w * W T w
and from what c w * is a vector of unknown constants, we get:
Yматрица (0) W T и
Y (0) W T w .
Note that another matrix-block derivation of the formulas is possible.
Transition from the system
U
u
Y (0) *
*
U
u
to the system
W
w
Y (0) *
*
W
w
can be realized by another method, replacing the line orthonormation of GY (0) g by the
following orthonormal decomposition of the matrix G
G T JL
where the matrix J has orthonormal columns, and the matrix L is upper triangular.
20
Then, taking into account the rule of transposition of matrices, we can write
G LT J T .
As a result, we get
GY (0) g , LT J T Y (0) g , J T Y (0) ( LT ) 1 g .
Here the rows of the matrix are orthonormal.
Comparing
W
w
Y (0) *
*
W
w
and
J T Y (0) ( LT ) 1 g
we get
w
W
u
J T , * ( LT ) 1 g ( LT ) 1 * .
*
W
w
u
Thus, we again obtain the orthonormal initial values of the unknown vector-valued
functions of the solution.
2.5. Conjugation of parts of the integration interval for S.K.Godunov’s sweep method.
To automate the computational process on the entire integration interval, which is
composed for conjugate shells with different physical and geometric parameters, the deformation
of which is described by different functions, it is necessary to have the procedures of conjugation
of the corresponding functions.
In the general case, the resolving functions of various parts of the integration interval of
the problem have no physical meaning, and the physical parameters of the problem are expressed
in various ways through these functions and their derivatives. At the same time, conjugation of
adjacent sections must satisfy kinematic and force conditions at the point of conjugation.
21
Solve the problem of conjugating parts of the interval of integration in the following way.
The vector P containing the physical parameters of the problem is formed using the
matrix M of coefficients and the required vector-function Y(x):
P MY (x)
where M is a square nondegenerate matrix.
Then at the conjugation point x = x* we can write expression
M Y- ( x * ) P M Y ( x * ) ,
where P is the vector corresponding to a discrete change in the physical parameters
when passing through the conjugation point from the left to the right; index "-" means "to the left
of the conjugation point", and the index "+" means "to the right of the conjugation point".
In S.K.Godunov’s method, the vector-function Y (x ) of the problem on each section is
sought in the form
Y ( x) Yматрица ( x)c Y ( x)
Suppose that the conjugation point does not coincide with the point of the orthogonal
transformation. Then the expression for conjugation conditions of adjacent sections
M Y- ( x * ) P M Y ( x * )
will take the form
M ( Y- матрица ( x)c - Y ( x)) P M ( Y матрица ( x)c Y ( x)) .
If now demand
c- c
then in the direct course of the sweep method, the integration can be continued from the left to
the right in the following expressions:
1
Y матрица ( x) M M Y- матрица ( x) ,
22
1
Y* ( x) M ( M Y-* ( x) P ) .
2.6. Properties of the transfer of boundary value conditions in S.K.Godunov’s sweep
method.
When solving a boundary-value problem for a system of "stiff" linear ordinary
differential equations by S.K.Godunov’s method says that discrete orthogonalization is carried
out by Gramm-Schmidt’s method with respect to vector-functions forming the variety of
solutions of the given problem in order to overcome the tendency of degeneration of these
vector-functions into linearly dependent ones.
At the same time, in the implementation of S.K.Godunov’s method, the boundary
conditions from the initial boundary are also transferred to the other edge. Let us show the
properties of this transfer.
Previously recorded
Y (0) Yматрица (0)c Y (0)
Y (0) W *T w * W T w
Yматрица (0) W T и
Y (0) W T w .
Then we can say that:
-
the vector w*, which is unknown, is a vector of constants c,
-
at the same time, the vector w* has the physical meaning of an external influence on the
deformed system unknown at the edge x=0,
-
the matrix W* is the matrix of boundary conditions unknown on the boundary x=0.
It follows from the formulated propositions that the transfer of boundary conditions in
S.K.Godunov’s method has the following meaning.
23
The continuation of integration, beginning with the vector Y (0) W T w , means the
transfer of the "convolution" W T w of the matrix equation of the boundary conditions at x=0 to
the right edge x=1.
The continuation of the integration, beginning with the vectors in the matrix Y матрица (0) ,
means that the matrix of the boundary conditions W*, which are unknown at the edge x=0, is
carried to the edge x=1.
Integration of differential equations is carried out with the goal of transferring the vector
c to the edge x=1, and hence the vector w*, which expresses the conditions unknown at the edge
x=0.
The transfer of the matrix W* and the vector w* means that the matrix equation
W *Y (0) w * of the boundary conditions, which are unknown at the edge x=0, is carried to the
edge x=1.
2.7. Modification of S.K.Godunov’s sweep method.
The solution in S.K.Godunov’s method is sought, as written above, in the form of the
formula
Y ( x) Yматрица ( x)c Y ( x) .
We can write this formula in two versions - in one case the formula satisfies the boundary
conditions of the left edge (index L), and in the other - the conditions on the right edge (index
R):
YL ( x) Y матрица L ( x)c L Y L ( x)
,
YR ( x) Y матрица R ( x)c R Y R ( x)
.
At an arbitrary point we have
24
YL ( x) YR ( x)
.
Then we obtain
Y матрица L ( x)c L Y L ( x) Y матрица R ( x)c R Y R ( x)
,
Y матрица L ( x)c L Y матрица R ( x)c R Y R ( x) Y L ( x)
,
Y матрица L ( x)
Y матрица R ( x)
cL
Y R ( x) Y L ( x)
.
cR
That is, a system of linear algebraic equations of the traditional kind with a square matrix
of coefficients Y матрица L ( x)
cL
is obtained.
cR
Y матрица R ( x) for the computation of the vectors of constants
25
Chapter 3. The method of "transferring of boundary value conditions" (the direct version
of the method) for solving boundary value problems with non-stiff ordinary differential
equations.
It is proposed to integrate by the formulas of the theory of matrices [Gantmakher]
immediately from some inner point of the interval of integration to the edges:
Y (0) K (0 x)Y ( x) Y (0 x) ,
Y (1) K (1 x)Y ( x) Y (1 x) .
We substitute the formula for Y (0) in the boundary conditions of the left edge and
obtain:
UY (0) u ,
U [ K (0 x)Y ( x) Y (0 x)] u ,
UK (0 x)Y ( x) u - UY (0 x) .
Similarly, for the right boundary conditions, we obtain:
VY (1) v ,
V [ K (1 x)Y ( x) Y (1 x)] v ,
VK (1 x)Y ( x) v VY (1 x) .
That is, we obtain two matrix equations of boundary conditions transferred to the point x
under consideration:
[UK (0 x)] Y ( x) u - UY (0 x) ,
[VK (1 x)] Y ( x) v VY (1 x) .
These equations are similarly combined into one system of linear algebraic equations
with a square matrix of coefficients to find the solution Y (x ) at any point x under
consideration:
26
UK (0 x)
u UY (0 x)
.
Y ( x)
VK (1 x)
v VY (1 x)
Chapter 4. The method of "additional boundary value conditions" for solving boundary
value problems with non-stiff ordinary differential equations.
Let us write on the left edge one more equation of the boundary conditions:
MY (0) m .
As matrix M rows, we can take those boundary conditions, that is, expressions of those
physical parameters that do not enter into the parameters of the boundary conditions of the left
edge U or are linearly independent with them. This is entirely possible, since for boundary value
problems there are as many independent physical parameters as the dimensionality of the
problem, and only half of the physical parameters of the problem enter into the parameters of the
boundary conditions.
That is, for example, if the problem of the shell of a rocket is considered, then on the left
edge 4 movements can be specified. Then for the matrix M we can take the parameters of forces
and moments, which are also 4, since the total dimension of such a problem is 8.
The vector m of the right side is unknown and it must be found, and then we can assume
that the boundary value problem is solved, that is, reduced to Cauchy’s problem, that is, the
vector Y (0) is found from the expression:
U
u
Y (0)
,
M
m
that is, the vector Y (0) is found from the solution of a system of linear algebraic equations with
a square non-degenerate coefficient matrix consisting of blocks U and M .
Similarly, we write on the right edge one more equation of the boundary conditions:
27
NY (1) n ,
where the matrix N is written from the same considerations for additional linearly independent
parameters on the right edge, and the vector n is unknown.
For the right edge, too, the corresponding system of equations is valid:
V
v
Y (1)
.
N
n
We write Y (1) K (1 0)Y (0) Y (1 0) and substitute it into the last system of linear
algebraic equations:
V
v
[K (1 0)Y (0) Y (1 0)]
,
N
n
V
v
V
K (1 0)Y (0)
Y (1 0) ,
N
n
N
V
v V Y (1 0)
,
K (1 0)Y (0)
N
n N Y (1 0)
V
s
K (1 0)Y (0)
.
N
t
We write the vector Y (0) through the inverse matrix:
U
Y (0)
M
1
u
m
and substitute it in the previous formula:
V
U
K (1 0)
N
M
1
u
s
m
t
Thus, we have obtained a system of equations of the form:
B
u
s
,
m
t
where the matrix B is known, the vectors u and s are known, and the vectors m and t are
unknown.
28
We divide the matrix B into 4 natural blocks for our case and obtain:
B11
B21
B12
u
s
,
B22 m
t
from which we can write that
B11u B12 m s,
B21u B22 m t.
Consequently, the required vector m is calculated by the formula:
m B121 ( s B11u)
And the required vector n is calculated through the vector t :
t B21u B22 m ,
n t N Y (1 0) .
29
Chapter 5. The method of "half of the constants" for solving boundary value problems
with non-stiff ordinary differential equations.
In this method we use the idea proposed by S.K.Godunov to seek a solution in the form
of only one-half of the possible unknown constants, but a formula for the possibility of starting
such a calculation and further application of matrix exponents (Cauchy’s matrices) are proposed
by A.Yu.Vinogradov.
The formula for starting calculations from the left edge with only one half of the possible
constants:
T
T
Y (0) U орто
uорто M орто
c,
T
Y (0) U орто
T
M орто
uорто
c
.
Thus, a formula is written in the matrix form for the beginning of the calculation from the
left edge, when the boundary conditions are satisfied on the left edge.
Then write VY (1) v and Y (1) K (1 0)Y (0) Y (1 0) collectively:
V [K (1 0)Y (0) Y (1 0)] v ,
VK (1 0)Y (0) v VY (1 0)
and substitute in this formula the expression for Y(0):
T
VK (1 0) U орто
T
M орто
uорто
c
v VY (1 0)
or
T
VK (1 0) U орто
T
M орто
uорто
Thus, we have obtained an expression of the form:
c
p.
30
D
uорто
c
p,
where the matrix D has a dimension of 4x8 and can be naturally represented in the form of two
square blocks of 4x4 dimension:
D1
D2
uорто
c
p.
Then we can write:
D1 uорто D2 c p .
Hence we obtain that:
c D21 ( p D1uорто ) .
Thus, the required constants are found.
31
Chapter 6. The method of "transferring of boundary conditions" (step-by-step version of
the method) for solving boundary value problems with stiff ordinary differential equations.
6.1. The method of "transfer of boundary value conditions" to any point of the interval of
integration.
The complete solution of the system of differential equations has the form
Y ( x) K ( x x0 )Y ( x0 ) Y ( x x0 ) .
Or you can write:
Y (0) K (0 x1 )Y ( x1 ) Y (0 x1 ) .
We substitute this expression for Y (0) into the boundary conditions of the left edge and
obtain:
UY (0) u ,
U [ K (0 x1 )Y ( x1 ) Y (0 x1 )] u ,
UK (0 x1 )Y ( x1 ) u UY (0 x1 ) .
Or we get the boundary conditions transferred to the point x1 :
U1Y ( x1 ) u1 ,
U1 UK (0 x1 )
where
and
u1 u UY (0 x1 ) .
Further, we write similarly
Y ( x1 ) K ( x1 x2 )Y ( x2 ) Y ( x1 x2 )
And substitute this expression for Y ( x1 ) into the transferred boundary conditions of the
point x1 :
U1Y ( x1 ) u1 ,
32
U1[ K ( x1 x2 )Y ( x2 ) Y ( x1 x2 )] u1 ,
U1 K ( x1 x2 )Y ( x2 ) u1 U1Y ( x1 x2 ) .
Or we get the boundary conditions transferred to the point x 2 :
U 2Y ( x2 ) u2 ,
where
U 2 U1 K ( x1 x2 )
and
u2 u1 U 1Y ( x1 x2 ) .
And so we transfer the matrix boundary condition from the left edge to the point x and
in the same way transfer the matrix boundary condition from the right edge.
Let us show the steps of transferring the boundary conditions of the right edge.
We can write:
Y (1) K (1 x n 1 )Y ( x n 1 ) Y (1 x n 1 )
We substitute this expression for Y (1) in the boundary conditions of the right edge and
obtain:
VY (1) v ,
V [ K (1 x n 1 )Y ( x n 1 ) Y (1 x n 1 )] v ,
VK (1 x n 1 )Y ( x n 1 ) v VY (1 x n 1 )
Or we get the boundary conditions of the right edge, transferred to the point x n 1 :
Vn1Y ( xn1 ) v n1 ,
where
Vn1 VK (1 xn1 )
and
v n 1 v VY (1 x n 1 ) .
Further, we write similarly
Y ( x n 1 ) K ( x n 1 x n 2 )Y ( x n 2 ) Y ( x n 1 x n 2 )
And substitute this expression for Y ( xn1 ) in the transferred boundary conditions of the
point x n 1 :
Vn1Y ( xn1 ) v n1 ,
33
Vn 1 [ K ( x n 1 x n 2 )Y ( x n 2 ) Y ( x n 1 x n 2 )] v n 1 ,
Vn 1 K ( x n 1 x n 2 )Y ( x n 2 ) v n 1 Vn 1Y ( x n 1 x n 2 ) .
Or we get the boundary conditions transferred to the point xn2 :
Vn2Y ( xn2 ) v n2 ,
where
Vn2 Vn1 K ( xn1 xn2 )
and
v n 2 v n 1 Vn 1Y ( x n 1 x n 2 ) .
And so at the inner point x of the integration interval we transfer the matrix boundary
condition, as shown, and from the left edge and in the same way transfer the matrix boundary
condition from the right edge and obtain:
U Y (x ) u ,
V Y (x ) v .
From these two matrix equations with rectangular horizontal coefficient matrices, we
obviously obtain one system of linear algebraic equations with a square matrix of coefficients:
U
u
.
Y
(x
)
V
v
6.2. The case of "stiff" differential equations.
In the case of "stiff" differential equations, it is proposed to apply a line
orthonormalization of the matrix boundary conditions in the process of their transfer to the point
under consideration. For this, the orthonormalization formulas for systems of linear algebraic
equations can be taken in [Berezin, Zhidkov].
That is, having received
U1Y ( x1 ) u1
we apply a line orthonormation to this group of linear algebraic equations and obtain an
equivalent matrix boundary condition:
34
U 1ортоY ( x1 ) u1орто .
And in this line orthonormal equation is substituted
Y ( x1 ) K ( x1 x2 )Y ( x2 ) Y ( x1 x2 ) .
And we get
U1орто[ K ( x1 x2 )Y ( x2 ) Y ( x1 x2 )] u1орто ,
U1орто K ( x1 x2 )Y ( x2 ) u1орто U1ортоY ( x1 x2 ) .
Or we get the boundary conditions transferred to the point x 2 :
U 2Y ( x2 ) u2 ,
where
U 2 U 1орто K ( x1 x 2 )
and
u2 u1орто U1ортоY ( x1 x2 ) .
Now we apply linear orthonorming to this group of linear algebraic equations and obtain
an equivalent matrix boundary condition:
U 2ортоY ( x 2 ) u2орто
And so on.
And similarly we do with intermediate matrix boundary conditions carried from the right
edge to the point under consideration.
As a result, we obtain a system of linear algebraic equations with a square matrix of
coefficients, consisting of two independently stepwise orthonormal matrix boundary conditions,
which is solved by Gauss’ method with the separation of the main element for obtaining the
solution Y ( x ) at the point x under consideration:
U орто
Vорто
Y (x )
uорто
v орто
.
35
6.3. Formulas for computing the vector of a particular solution of inhomogeneous system of
differential equations.
Instead of the formula for computing the vector of a particular solution of an
inhomogeneous system of differential equations in the form [Gantmaher]:
x
Y ( x x0 ) e Ax e At F (t )dt
x0
it is proposed to use the following formula for each individual section of the integration interval:
xj
Y ( x j xi ) Y ( x j xi ) K ( x j xi ) K ( xi t ) F (t )dt .
xi
The correctness of the above formula is confirmed by the following:
xj
Y ( x j xi ) exp( A( x j xi )) exp( A( xi t )) F (t )dt ,
xi
xj
Y ( x j xi )
exp( A( x
j
xi )) exp( A( xi t )) F (t )dt ,
xi
xj
Y ( x j xi )
exp( A( x
j
xi xi t )) F (t )dt ,
xi
xj
Y ( x j xi )
exp( A( x
j
t )) F (t )dt ,
xi
xj
Y ( x j xi ) exp( Ax j ) exp( At ) F (t )dt ,
xi
x
Y ( x xi ) exp( Ax) exp( At ) F (t )dt ,
xi
36
which was to be confirmed.
The calculation of the vector of a particular solution of a system of differential equations
is performed using the representation of Cauchy’s matrix under the integral sign in the form of a
series and integrating this series elementwise:
xj
Y ( x j xi ) Y ( x j xi ) K ( x j xi ) K ( xi t ) F (t )dt
xi
xj
K ( x j xi ) ( E A( xi t ) A 2 ( xi t ) 2 / 2!...) F (t )dt
xi
xj
xj
xj
K ( x j xi )( E F (t )dt A ( xi t ) F (t )dt A / 2! ( xi t ) 2 F (t )dt ...).
xi
xi
xi
2
This formula is valid for the case of a system of differential equations with constant
coefficient matrix A =const.
Let us consider the variant, when the steps of the integration interval are chosen
sufficiently small, which allows us to consider the vector F (t ) in the region ( x j xi )
approximately as a constant F ( xi ) constant , which allows us to remove this vector from the
signs of the integrals:
xj
xj
xj
Y ( x j xi ) K ( x j xi )( E dt A ( xi t )dt A 2 / 2! ( xi t ) 2 dt ...) F (t ).
xi
xi
xi
It is known that when T=(at+b) we have T n dt
In our case, we have
(b - t)
n
dt
xj
Then we obtain
(x
xi
i
t ) n dt
1
T n 1 const (при n -1).
a(n 1)
1
(b - t) n 1 const (при n -1).
(-1)(n 1)
1
( xi x j ) n 1 .
n 1
37
Then we obtain a series for computing the vector of a particular solution of an
inhomogeneous system of differential equations on a small section ( x j xi ) :
Y ( x j xi ) K ( x j xi ) ( E A( xi x j ) / 2! A2 ( xi x j ) 2 / 3!...) ( x j xi ) F ( xi ).
For the case of differential equations with variable coefficients, an averaged matrix
Ai A( xi ) of the coefficients of the system of differential equations can be used for each
section.
If the considered section of the integration interval is not small, then the following
iterative (recurrent) formulas are proposed.
We give the formulas for computing the vector of a particular solution, for example,
Y ( x3 x0 ) on the considered section ( x3 x0 ) through the vectors of the particular solution
Y ( x1 x 0 ) , Y ( x 2 x1 ) , Y ( x3 x 2 ) , corresponding to the subsections ( x1 x0 ) ,
( x2 x1 ) , ( x3 x2 ) .
We have Y ( x) K ( x x0 )Y ( x0 ) Y ( x x0 ) .
We also have a formula for a separate subsection:
xj
Y ( x j xi ) Y ( x j xi ) K ( x j xi ) K ( xi t ) F (t )dt .
xi
We can write:
Y ( x1 ) K ( x1 x0 )Y ( x0 ) Y ( x1 x0 ) ,
Y ( x2 ) K ( x2 x1 )Y ( x1 ) Y ( x2 x1 ) .
We substitute Y ( x1 ) in Y ( x2 ) and get:
Y ( x 2 ) K ( x 2 x1 )[ K ( x1 x0 )Y ( x0 ) Y ( x1 x0 )] Y ( x 2 x1 )
K ( x 2 x1 ) K ( x1 x0 )Y ( x0 ) K ( x 2 x1 )Y ( x1 x0 ) Y ( x 2 x1 ) .
Let us compare the expression obtained with the formula:
38
Y ( x 2 ) K ( x 2 x0 )Y ( x0 ) Y ( x 2 x0 )
and we get, obviously, that:
K ( x2 x0 ) K ( x2 x1 ) K ( x1 x0 )
and for the particular vector we obtain the formula:
Y ( x 2 x0 ) K ( x 2 x1 )Y ( x1 x0 ) Y ( x 2 x1 ) .
That is, the subsector vectors Y ( x1 x0 ), Y ( x 2 x1 ) are not simply add with each
other, but with the participation of Cauchy’s matrix of the subsection.
Similarly, we write down Y ( x3 ) K ( x3 x 2 )Y ( x 2 ) Y ( x3 x 2 ) and substitute the
formula for Y ( x2 ) and get:
Y ( x3 ) K ( x3 x 2 )[ K ( x 2 x1 ) K ( x1 x0 )Y ( x0 ) K ( x 2 x1 )Y ( x1 x0 ) Y ( x 2 x1 )]
Y ( x3 x 2 ) K ( x3 x 2 ) K ( x 2 x1 ) K ( x1 x0 )Y ( x0 )
K ( x3 x 2 ) K ( x 2 x1 )Y ( x1 x0 ) K ( x3 x 2 )Y ( x 2 x1 ) Y ( x3 x 2 ).
Comparing the expression obtained with the formula:
Y ( x3 ) K ( x3 x0 )Y ( x0 ) Y ( x3 x0 )
obviously, we get that:
K ( x3 x0 ) K ( x3 x2 ) K ( x2 x1 ) K ( x1 x0 )
and together with this we get the formula for a particular vector:
That is, in this way a particular vector is calculated - the vector of the particular solution
of the inhomogeneous system of differential equations, that is, for example, a particular vector
Y ( x3 x 0 )
vectors
on the considered section
( x3 x 0 )
Y ( x1 x 0 ) , Y ( x 2 x1 ) , Y ( x3 x 2 )
, ( x2 x1 ) , ( x3 x2 ) .
is calculated through the computed partial
corresponding to the sub-sections
( x1 x0 )
39
6.4. Applicable formulas for orthonormalization.
Taken from [Berezin, Zhidkov]. Let there be given a system of linear algebraic equations
of order n:
A x =b .
Here, above the vectors, we draw dashes instead of their designation in boldface.
We will consider the rows of the matrix A of the system as vectors:
a i =( a i1 , a i2 ,…, a in ).
We orthonormalize this system of vectors.
n
a12k .
The first equation of the system A x = b we divide by
k 1
In doing so, we get:
с11 x1 + с12 x 2 +…+ с1n x n = d 1 , c1 =( c11 , c12 ,…, c1n ),
where
c1k =
a1k
,
n
n
b1
d1=
n
k 1
a12k
a12k
k 1
c12k =1.
,
k 1
The second equation of the system is replaced by:
с 21 x1 + с 22 x 2 +…+ с 2n x n = d 2 , c 2 =( c 21 , c 22 ,…, c 2n ),
/
where
c 2k =
c 2k
n
/
,
d 2=
n
,
c 2/ 2k
c 2/ 2k
k 1
/
d2
k 1
/
c 2k = a 2k -( a 2 , c1 ) c1k , d 2 = b 2 -( a 2 , c1 ) d 1 .
Similarly we proceed further. The equation with the number i takes the form:
40
с i1 x1 + с i2 x 2 +…+ с in x n = d i , ci =( c i1 , c i2 ,…, c in ),
/
where
c ik =
cik
n
cik/ 2
k 1
/
,
di
di=
n
,
cik/ 2
k 1
/
cik = a ik -( a i , c1 ) c1k -( a i , c 2 ) c 2k -…-( a i , ci 1 ) ci 1, k ,
/
d i = bi -( a i , c1 ) d 1 -( a i , c 2 ) d 2 -…-( a i , ci 1 ) d i 1 .
Here ( a i , c j ) is the scalar multiplication of vectors.
The process will be realized if the system of linear algebraic equations is linearly
independent.
As a result, we come to a new system C x d , where the matrix C will be with
orthonormal rows, that is, it has the property C C T E , where E is the identity matrix.
41
Chapter 7. The simplest method for solving boundary value problems with stiff ordinary
differential equations without orthonormalization - the method of "conjugation of sections
of the integration interval", which are expressed by matrix exponents.
The idea of overcoming the difficulties of computation by dividing the interval of integration
into conjugate areas belongs to Dr.Sc. Professor Yu.I.Vinogradov (his doctoral thesis was
defended including on this idea), and the simplest realization of this idea through the formulas of
the theory of matrices belongs to the Ph.D. A.Yu.Vinogradov.
We divide the interval of integration of the boundary value problem, for example, into 3
sections. We will have points (nodes), including edges:
x0 , x1 , x2 , x3 .
We have boundary conditions in the form:
UY ( x0 ) u,
VY ( x3 ) v.
We can write the matrix equations of conjugation of sections:
Y ( x0 ) K ( x0 x1 )Y ( x1 ) Y ( x0 x1 ) ,
Y ( x1 ) K ( x1 x2 )Y ( x2 ) Y ( x1 x2 ) ,
Y ( x 2 ) K ( x 2 x3 )Y ( x3 ) Y ( x 2 x3 ) .
We can rewrite it in a form more convenient for us further:
EY ( x0 ) K ( x0 x1 )Y ( x1 ) Y ( x0 x1 ) ,
EY ( x1 ) K ( x1 x2 )Y ( x2 ) Y ( x1 x2 ) ,
EY ( x 2 ) K ( x 2 x3 )Y ( x3 ) Y ( x 2 x3 ) .
where E is the identity matrix.
42
Then in a combined matrix form we obtain a system of linear algebraic equations in the
following form:
U
E
0
0
0
0
0
0
u
Y ( x0 )
K ( x0 x1 )
0
0
Y ( x0 x1 )
Y ( x1 )
Y ( x1 x 2 ) .
E
K ( x1 x 2 )
0
Y ( x2 )
0
E
K ( x 2 x3 )
Y ( x 2 x3 )
Y ( x3 )
0
0
V
v
This system is solved by Gauss’ method with the separation of the main element.
At points located between nodes, the solution is to be solved by solving Cauchy’s
problems with the initial conditions in the i-th node:
Y ( x) K ( x xi )Y ( xi ) Y ( x xi ) .
It is not necessary to apply orthonormalization for boundary value problems for stiff
ordinary differential equations, since on each section of the integration interval the calculation of
each matrix exponent is fulfilled independently and from the initial orthonormal identity matrix,
which makes it unnecessary to use orthonormalization, unlike S.K.Godunov’s method, which
greatly simplifies programming in comparison with S.K.Godunov's method.
It is possible to calculate Cauchy’s matrices not in the form of matrix exponents, but with
Runge-Kutta’s methods from the starting identity matrix, and the vector of the particular solution
of the inhomogeneous system of differential equations can be calculated on each site by RungeKutta’s methods from the starting zero vector. In the case of Runge-Kutta’s methods, error
estimates are well known, which means that calculations can be performed with a known
accuracy.
43
Chapter 8. Calculation of shells of composite and with frames by the simplest method of
"conjugation of sections of the integration interval".
8.1. The variant of recording of the method for solving stiff boundary value problems
without orthonormalization - the method of "conjugation of sections, expressed by matrix
exponents "- with positive directions of matrix formulas of integration of differential
equations.
We divide the interval of integration of the boundary value problem, for example, into 3
sections. We will have points (nodes), including edges:
x0 , x1 , x2 , x3 .
We have boundary conditions in the form:
UY ( x0 ) u,
VY ( x3 ) v.
We can write the matrix equations of conjugation of sections:
Y ( x1 ) K ( x1 x0 )Y ( x0 ) Y ( x1 x0 ) ,
Y ( x2 ) K ( x2 x1 )Y ( x1 ) Y ( x2 x1 ) ,
Y ( x3 ) K ( x3 x 2 )Y ( x 2 ) Y ( x3 x 2 ) .
We can rewrite it in a form more convenient for us further:
EY ( x1 ) K ( x1 x0 )Y ( x0 ) Y ( x1 x0 ) ,
EY ( x2 ) K ( x2 x1 )Y ( x1 ) Y ( x2 x1 ) ,
EY ( x3 ) K ( x3 x 2 )Y ( x 2 ) Y ( x3 x 2 ) .
where E is the identity matrix.
As a result, we obtain a system of linear algebraic equations:
44
U
0
0
0
u
Y ( x0 )
K ( x1 x0 )
E
0
0
Y ( x1 x0 )
Y ( x1 )
Y ( x 2 x1 ) .
0
K ( x 2 x1 )
E
0
Y ( x2 )
0
0
K ( x3 x 2 ) E
Y ( x3 x 2 )
Y ( x3 )
0
0
0
V
v
This system is solved by Gauss’ method with the separation of the main element.
It turns out that it is not necessary to apply orthonormalization, since sections of the
integration interval are chosen so long that the computation on them is stable.
At points near the nodes, the solution is found by solving the corresponding Cauchy’s
problems with the origin at the i-th node:
Y ( x) K ( x xi )Y ( xi ) Y ( x xi ) .
8.2. Composite shells of rotation.
Let us consider the conjugation of segments of the composite shell of rotation.
Suppose we have 3 sections, where each section can be expressed by its differential
equations and the physical parameters can be expressed differently - different formulas on
different sections:
In the general case (for the example of section 12), the physical parameters of the section
(vector P12 ( x) ) are expressed in terms of the required parameters of the system of ordinary
differential equations of this section (through the vector Y12 ( x) ) as follows:
P12 ( x) M12Y12 ( x) ,
where the matrix M12 is a square non-degenerate matrix.
45
With the transition of the conjugation point, we can write in a general form (but using the
conjugation point x1 ):
P01 ( x1 ) P0112 L0112 P12 ( x1 ) ,
where P0112 is the discrete increment of physical parameters (forces, moments) during the
transition from the "01" section to the "12" section, and the square nondegenerate matrix L0112 is
diagonal and consists of units and minus ones on the main diagonal to establish the correct
correspondence between the positive directions of forces, angular momenta, displacements and
angles when going from "01" to "12", which may be different (in different differential equations
of different conjugate regions) in the equations to the left of the conjugation point and in the
equations to the right of the conjugation point.
The last two equations combine to form the equation:
M 01Y01 ( x1 ) P0112 L0112 M 12Y12 ( x1 ) .
At the conjugation point x 2 , we similarly obtain the equation:
M 12Y12 ( x2 ) P1223 L1223 M 23Y23 ( x2 ) .
If the shell consisted of identical parts, then we could write in a combined matrix form a
system of linear algebraic equations in the following form:
U
0
0
0
u
Y ( x0 )
K ( x1 x0 )
E
0
0
Y ( x1 x0 )
Y ( x1 )
Y ( x 2 x1 ) .
0
K ( x 2 x1 )
E
0
Y ( x2 )
0
0
K ( x3 x 2 ) E
Y ( x3 x 2 )
Y ( x3 )
0
0
0
V
v
But in our case the shell consists of 3 sections, where the middle section can be
considered, for example, a frame expressed in terms of its differential equations.
Then instead of vectors Y ( x0 ) , Y ( x1 ) , Y ( x2 ) , Y ( x3 ) we should consider vectors:
Y01 ( x0 ), Y01 ( x1 ), Y12 ( x1 ), Y12 ( x2 ), Y23 ( x2 ), Y23 ( x3 ) .
46
Then the matrix equations
UY ( x0 ) u,
VY ( x3 ) v.
EY ( x1 ) K ( x1 x0 )Y ( x0 ) Y ( x1 x0 ) ,
EY ( x2 ) K ( x2 x1 )Y ( x1 ) Y ( x2 x1 ) ,
EY ( x3 ) K ( x3 x 2 )Y ( x 2 ) Y ( x3 x 2 )
will take the form:
UY01 ( x0 ) u,
VY23 ( x3 ) v.
EY01 ( x1 ) K 01 ( x1 x0 )Y01 ( x0 ) Y01* ( x1 x0 ) ,
M 01Y01 ( x1 ) P0112 L0112 M 12Y12 ( x1 ) ,
EY12 ( x2 ) K12 ( x2 x1 )Y12 ( x1 ) Y12* ( x2 x1 ) ,
M 12Y12 ( x2 ) P1223 L1223 M 23Y23 ( x2 ) ,
EY23 ( x3 ) K 23 ( x3 x 2 )Y23 ( x 2 ) Y23* ( x3 x 2 ) .
After rearranging the summands, we get:
UY01 ( x0 ) u,
VY23 ( x3 ) v.
K 01 ( x1 x0 )Y01 ( x0 ) EY01 ( x1 ) Y01* ( x1 x0 ) ,
M 01Y01 ( x1 ) L0112 M 12Y12 ( x1 ) P0112 ,
K12 ( x2 x1 )Y12 ( x1 ) EY12 ( x2 ) Y12* ( x2 x1 ) ,
M 12Y12 ( x2 ) L1223 M 23Y23 ( x2 ) P1223 ,
47
K 23 ( x3 x 2 )Y23 ( x 2 ) EY23 ( x3 ) Y23* ( x3 x 2 ) .
As a result, we can write down the final system of linear algebraic equations:
U
0
K 01 ( x1 x0 ) E
0
M 01
0
0
0
0
0
0
0
0
0
0
0
0
L0112 M 12
0
K 12 ( x 2 x1 ) E
0
M 12
0
0
0
0
0
0
0
0
0
u
Y01 ( x0 )
*
0
Y01 ( x1 x0 )
Y01 ( x1 )
0
P0112
Y12 ( x1 )
*
Y12 ( x 2 x1 )
0
Y12 ( x 2 )
L12 23 M 23
0
P12 23
Y23 ( x 2 )
*
K 23 ( x3 x 2 ) E
Y23 ( x3 x 2 )
Y23 ( x3 )
0
V
v
This system is solved by the Gauss method with the separation of the main element.
At points located between nodes, the solution is to be solved by solving Cauchy’s
problems with the initial conditions in the i-th node:
Y ( x) K ( x xi )Y ( xi ) Y ( x xi ) .
It is not necessary to apply orthonormalization for boundary value problems for stiff
ordinary differential equations.
8.3. Frame, expressed not by differential, but algebraic equations.
Let us consider the case when the frame (at a point x1 ) is expressed not in terms of
differential equations, but in terms of algebraic equations.
Above we wrote down that:
P01 ( x1 ) P0112 L0112 P12 ( x1 )
We can represent the vector P01 ( x1 ) of force factors and displacements in the form:
P01 ( x1 )
R01 ( x1 )
,
S01 ( x1 )
48
where R01 ( x1 ) is the displacement vector, S01 ( x1 ) is the vector of forces and moments.
Algebraic equation for the frame:
GR S ,
where G is the matrix of the rigidity of the frame, R is the vector of the frame movements, S is
the vector of force factors that act on the frame.
At the point of the frame we have:
R 0, S GR ,
that is, there is no discontinuity in the movements R 0 , but there is a resultant vector of force
factors S GR , which consists of forces and moments on the left plus forces and moments to
the right of the point of the frame.
P01 ( x1 )
R
L0112 P12 ( x1 ) ,
S
P01 ( x1 )
0
L0112 P12 ( x1 ) ,
GR
R01 ( x1 )
R (x )
0
L0112 12 1 ,
S 01 ( x1 ) GR
S12 ( x1 )
M 01Y01 ( x1 )
0
L0112 M 12Y12 ( x1 ) ,
GR
M 01Y01 ( x1 ) g * L0112 M 12Y12 ( x1 ) , где g *
0
,
GR
which is true if we do not forget that in this case we have:
P01 ( x1 )
R01 ( x1 )
,
S 01 ( x1 )
that is, the vector of displacements and force factors is first compiled from displacements
(above) R01 ( x1 ) , and then from force factors (below) S01 ( x1 ) .
49
Here it is necessary to remember that the displacement vector R01 ( x1 ) is expressed in
terms of the required state vector Y01 ( x1 ) :
g*
0
0
,
GR01 ( x1 )
GR
R01 ( x1 )
M 11p
p
P01 ( x1 )
M 01Y01 ( x1 ) M Y01 ( x1 )
S 01 ( x1 )
M 21p
M 12p
Y01 ( x1 ) ,
M 22p
where for convenience was introduced re-designation M 01 M p .
Then we can write:
R01 ( x1 ) M 11p
g*
0
GR01 ( x1 )
G M 11p
M 12p Y01 ( x1 ) ,
0
0...0
0
Y
(
x
)
Y01 ( x1 )
01
1
M 12p Y01 ( x1 )
G M 11p M 12p
G M 11p M 12p
We write the matrix equations for this case:
UY01 ( x0 ) u,
VY12 ( x 2 ) v.
EY01 ( x1 ) K 01 ( x1 x0 )Y01 ( x0 ) Y01* ( x1 x0 ) ,
M 01Y01 ( x1 ) g * L0112 M 12Y12 ( x1 ) ,
EY12 ( x2 ) K12 ( x2 x1 )Y12 ( x1 ) Y12* ( x2 x1 ) .
Let us write the vector g * in the equation:
M 01Y01 ( x1 )
( M 01
0
G M 11p
M 12p
0
GM
p
11
M 12p
Y01 ( x1 ) L0112 M 12Y12 ( x1 ) ,
) Y01 ( x1 ) L0112 M 12Y12 ( x1 ) .
50
To ensure the non-cumbersomeness, we introduce the notation:
( M 01
0
GM
p
11
M
p
12
) M *.
Then equation
M 01Y01 ( x1 ) g * L0112 M 12Y12 ( x1 )
will take the form:
M *Y01 ( x1 ) L0112 M 12Y12 ( x1 ) .
For convenience, we rearrange the terms in the matrix equations so that the resulting
system of linear algebraic equations is written clearly:
UY01 ( x0 ) u,
VY12 ( x 2 ) v.
K 01 ( x1 x0 )Y01 ( x0 ) EY01 ( x1 ) Y01* ( x1 x0 ) ,
M *Y01 ( x1 ) L0112 M 12Y12 ( x1 ) 0 ,
K12 ( x2 x1 )Y12 ( x1 ) EY12 ( x2 ) Y12* ( x2 x1 ) .
Thus, we obtain the resulting system of linear algebraic equations:
U
0
K 01 ( x1 x0 ) E
0
M*
0
0
0
0
0
0
0
u
Y01 ( x0 )
*
0
Y01 ( x1 x0 )
Y01 ( x1 )
.
L0112 M 12
0
0
Y12 ( x1 )
*
K12 ( x 2 x1 ) E
Y12 ( x 2 x1 )
Y12 ( x 2 )
0
V
v
If an external force-moment action g p is applied to the frame, then g *
should be rewritten in the form g *
g*
GM
p
11
0
GR
0
0
, then:
p
GR01 ( x1 ) g p
GR g
0
0...0
0
.
p
p
p Y01 ( x1 )
M Y01 ( x1 ) g
G M 11 M 12
gp
p
12
51
Then the matrix equation
M 01Y01 ( x1 )
0
GM
p
11
M 12p
Y01 ( x1 ) L0112 M 12Y12 ( x1 )
will take the form:
M 01Y01 ( x1 )
0
GM
p
11
M
Y01 ( x1 )
p
12
0
L0112 M 12Y12 ( x1 ) ,
gp
M *Y01 ( x1 ) L0112 M 12Y12 ( x1 )
0
.
gp
The resulting system of linear algebraic equations takes the form:
U
0
K 01 ( x1 x0 ) E
0
M*
0
0
0
0
u
0
Y01 ( x0 )
Y01* ( x1 x0 )
0
0
Y01 ( x1 )
p
.
L0112 M 12
0
g
Y12 ( x1 )
K12 ( x 2 x1 ) E
Y12* ( x 2 x1 )
Y12 ( x 2 )
0
V
v
0
0
8.4. The case where the equations (of shells and frames) are expressed not with abstract
vectors, but with vectors, consisting of specific physical parameters.
Let us consider the case when the parts of the shell structure and the frame are expressed
with state vectors (such as Y12 ( x) ), which (in a particular case) coincide with vectors of physical
parameters (such as P12 ( x) - displacements, angles, forces, moment).
Then the matrices of
type M12 are unit: M 12 E . And let the positive directions of the physical parameters be the
same for all parts of the shell and the frame ( L0112 E ).
Then we have the equations:
P12 ( x) M12Y12 ( x) ,
P01 ( x1 ) P0112 L0112 P12 ( x1 ) ,
52
M 01Y01 ( x1 ) P0112 L0112 M 12Y12 ( x1 ) ,
as:
P12 ( x) EY12 ( x) ,
P01 ( x1 ) P0112 EP12 ( x1 ) ,
EY01 ( x1 ) P0112 EY12 ( x1 ) ,
where E is the identity matrix.
Equations
M 01Y01 ( x1 ) g * L0112 M 12Y12 ( x1 ) ,
P01 ( x1 )
P01 ( x1 )
0
L0112 P12 ( x1 ) ,
GR01 ( x1 ) g p
R01 ( x1 )
Mp
M 01Y01 ( x1 ) M pY01 ( x1 ) 11p
S 01 ( x1 )
M 21
R01 ( x1 ) M 11p
g*
GM
p
11
M 12p
Y01 ( x1 ) ,
M 22p
M 12p Y01 ( x1 ) ,
0
0
0
0
p
Y01 ( x1 ) p ,
p
p
p
M 12 Y01 ( x1 )
G M 11 M 12
g
g
will take the form:
EY01 ( x1 ) g * EY12 ( x1 ) ,
EY01 ( x1 )
P01 ( x1 )
0
EY12 ( x1 ) ,
GR01 ( x1 ) g p
R01 ( x1 )
E 0
EY01 ( x1 ) M pY01 ( x1 )
Y01 ( x1 ) ,
S 01 ( x1 )
0 E
R01 ( x1 ) E 0 Y01 ( x1 ) ,
g*
0
GR01 ( x1 )
0
0
0
0
0
Y01 ( x1 ) p .
p
p
G E 0 Y01 ( x1 )
GE 0
g
g
g
53
Equations
M 01Y01 ( x1 )
0
GM
p
11
M
Y01 ( x1 )
p
12
0
L0112 M 12Y12 ( x1 ) ,
gp
M *Y01 ( x1 ) L0112 M 12Y12 ( x1 )
0
.
gp
will take the form:
EY01 ( x1 )
0
0
Y01 ( x1 ) p EY12 ( x1 ) ,
GE 0
g
0
0
,
где
(
E
) M*
p
GE 0
g
M *Y01 ( x1 ) EY12 ( x1 )
The resulting system of linear algebraic equations
U
0
K 01 ( x1 x0 ) E
0
M*
0
0
0
0
u
0
*
Y01 ( x0 )
Y01 ( x1 x0 )
0
0
Y01 ( x1 )
p
L0112 M 12
0
g
Y12 ( x1 )
K12 ( x 2 x1 ) E
*
Y12 ( x 2 x1 )
Y12 ( x 2 )
0
V
v
0
0
will take the form:
U
0
K 01 ( x1 x0 ) E
0
M*
0
0
0
0
where M ( E
*
8 x8
0
4 x8
G E
4x4 4x4
0
u
0
0
*
Y01 ( x0 )
Y01 ( x1 x0 )
0
0
0
Y01 ( x1 )
p
,
E
0
g
Y12 ( x1 )
K12 ( x 2 x1 ) E
Y12* ( x 2 x1 )
Y12 ( x 2 )
0
V
v
) (E
4x4
This means that equation
8 x8
0
4 x8
G
4x4
0
4x4
)
E
4x4
G
4x4
0
4x4
E
4x4
.
54
M *Y01 ( x1 ) EY12 ( x1 )
0
gp
will take the form:
E
4x4
G
4x4
E
4x4
G
4x4
0
4x4
E
Y01 ( x1 ) EY12 ( x1 )
4x4
0 R01 ( x1 )
E
4x4
E S 01 ( x1 )
0
4x4
4x4
4x4
0
gp
0 R12 ( x1 )
0
p
E S12 ( x1 )
g
4x4
4x4
R01 ( x1 ) R12 ( x1 ) 0 , (there is no jump in the displacements and the angle) and
GR01 ( x1 ) S 01 ( x1 ) S12 ( x1 ) g p - base balance,
i.e:
R01 ( x1 ) R12 ( x1 ) (displacement and angle: no rupture)
S 01 ( x1 ) S g p S12 ( x1 ) , где S GR01 ( x1 ) (strength and momentum: balance).
55
Appendix. Computational experiments (a C++ program).
Computational experiments have been carried out in comparison with the method of the
boundary conditions transfer of Alexei Vinogradov. Line-by-line orthonormalization is used in
this method.
Without use of orthonormalization it is possible by means of the boundary conditions
transfer method to solve a problem of cylindrical shell loading that is fixed in cantilever fashion
at the right boundary and loaded at the left boundary by the force distributed uniformly by the
circular arch with the ratio of the length to the radius L/R=2 and with ratio of the radius to the
thickness R/h=100. In case of ratio R/h=200 the problem by means of the method of the
boundary conditions transfer without orthonormalozation cannot be solved by this time because
there are mistakes resulted in due to counting instability. However, in case of use of
orthonormalization in the method of the boundary conditions transfer the problems related to the
parameters R/h=300, R/h=500, R/h=1000 can be solved.
A new method proposed in this paper allows solving of all above mentioned test
problems without use of orthonormalization operation that results in significant simplification of
its programming.
In case of the test computations of the problems characterized by the above mentioned
parameters by means of this new proposed method the integration interval is divided into 10
segments while between the nodes as aforesaid the solution was found as a solution of Cauchy’s
problem. 50 harmonics of Fourier’s series were used for solving the problem since the result in
case of usage of 50 harmonics didn’t differ from the case when 100 harmonics were used.
Test problem computing speed by means of the proposed method is not less compared to
the boundary conditions transfer method because both methods when used for test problems
while using 50 harmonics of the Fourier series produced a final solution instantaneously after
launching a program (notebook ASUS M51V CPU Duo T5800). At the same time programming
of this newly proposed method is significantly simpler because there is no need in
orthonormalization procedure programming.
C++ program
// sopryazhenie.cpp: main file of the project.
//Solution of the boundary value problem – a cylindrical shell problem.
//Integration interval is divided into 10 matching segments: left boundary – point 0 and
right boundary – point 10.
//WITHOUT ORTHONORMALIZATION
#include "stdafx.h"
#include <iostream>
56
#include <conio.h>
using namespace std;
//Multiplication of A matrix by b vector and obtain rezult.
void mat_on_vect(double A[8][8], double b[8], double rezult[8]){
for(int i=0;i<8;i++){
rezult[i]=0.0;
for(int k=0;k<8;k++){
rezult[i]+=A[i][k]*b[k];
}
}
}
//Computation of the matrix exponent EXP=exp(A*delta_x)
void exponent(double A[8][8], double delta_x, double EXP[8][8]) {
//n – number of the terms of the series in the exponent, m – a counter of the
number of the terms of the series (m<=n)
int n=100, m;
double E[8][8]={0}, TMP1[8][8], TMP2[8][8];
int i,j,k;
//E – unit matrix – the first term of the series of the exponent
E[0][0]=1.0; E[1][1]=1.0; E[2][2]=1.0; E[3][3]=1.0;
E[4][4]=1.0; E[5][5]=1.0; E[6][6]=1.0; E[7][7]=1.0;
//initial filling-in of the auxiliary array TMP1 – the previous term of the series
for follow-up multiplication
//and initial filling-in of the exponent by the first term of the series
for(i=0;i<8;i++) {
for(j=0;j<8;j++) {
TMP1[i][j]=E[i][j];
EXP[i][j]=E[i][j];
}
}
//series of EXP exponent computation starting from the 2nd term of the series
(m=2;m<=n)
for(m=2;m<=n;m++) {
for(i=0;i<8;i++) {
for(j=0;j<8;j++) {
TMP2[i][j]=0;
for(k=0;k<8;k++) {
//TMP2[i][j]+=TMP1[i][k]*A[k][j]*delta_x/(m-1);
TMP2[i][j]+=TMP1[i][k]*A[k][j];
}
TMP2[i][j]*=delta_x;//taken out beyond the cycle of
multiplication of the row by the column
TMP2[i][j]/=(m-1);// taken out beyond the cycle of
multiplication of the row by the column
EXP[i][j]+=TMP2[i][j];
}
}
//filling-in of the auxiliary array TMP1 for computing the next term of the
series TMP2 in the next step of the cycle by m
if (m<n) {
for(i=0;i<8;i++) {
for(j=0;j<8;j++) {
TMP1[i][j]=TMP2[i][j];
}
57
}
}
}
}
//computation of the matrix MAT_ROW in the form of the matrix series for follow-up use
//when computing a vector - partial vector – a vector of the partial solution of the
heterogeneous system of the ordinary differential equations at the step delta x
void mat_row_for_partial_vector(double A[8][8], double delta_x, double MAT_ROW[8][8]) {
//n – number of the terms of the series in MAT_ROW, m – a counter of the number of
the terms of the series (m<=n)
int n=100, m;
double E[8][8]={0}, TMP1[8][8], TMP2[8][8];
int i,j,k;
//E – unit matrix – the first term of the series MAT_ROW
E[0][0]=1.0; E[1][1]=1.0; E[2][2]=1.0; E[3][3]=1.0;
E[4][4]=1.0; E[5][5]=1.0; E[6][6]=1.0; E[7][7]=1.0;
//initial filling-in of the auxiliary array TMP1 – the previous term of the series
for following multiplication
//and initial filling-in of MAT_ROW by the first term of the series
for(i=0;i<8;i++) {
for(j=0;j<8;j++) {
TMP1[i][j]=E[i][j];
MAT_ROW[i][j]=E[i][j];
}
}
//a series of computation of MAT_ROW starting from the second term of the series
(m=2;m<=n)
for(m=2;m<=n;m++) {
for(i=0;i<8;i++) {
for(j=0;j<8;j++) {
TMP2[i][j]=0;
for(k=0;k<8;k++) {
TMP2[i][j]+=TMP1[i][k]*A[k][j];
}
TMP2[i][j]*=delta_x;
TMP2[i][j]/=m;
MAT_ROW[i][j]+=TMP2[i][j];
}
}
//filling-in of the auxiliary array TMP1 for computing the next term of the
series – TMP2 in the next step of the cycle by m
if (m<n) {
for(i=0;i<8;i++) {
for(j=0;j<8;j++) {
TMP1[i][j]=TMP2[i][j];
}
}
}
}
}
//specifying the external influence vector in the system of ordinary differential
equations – POWER vector: Y'(x)=A*Y(x)+POWER(x):
void power_vector_for_partial_vector(double x, double POWER[8]){
POWER[0]=0.0;
58
POWER[1]=0.0;
POWER[2]=0.0;
POWER[3]=0.0;
POWER[4]=0.0;
POWER[5]=0.0;
POWER[6]=0.0;
POWER[7]=0.0;
}
//computation of the vector – ZERO (particular case) vector of the partial solution
//heterogeneous system of differential equations in the segment of interest:
void partial_vector(double vector[8]){
for(int i=0;i<8;i++){
vector[i]=0.0;
}
}
//computation of the vector – the vector of the partial solution of the heterogeneous
system of differential equations in the segment of interest delta x:
void partial_vector_real(double expo_[8][8], double mat_row[8][8], double
x_, double delta_x, double vector[8]){
double POWER_[8]={0};//External influence vector on the shell
double REZ[8]={0};
double REZ_2[8]={0};
power_vector_for_partial_vector(x_, POWER_);//Computing POWER_ at coordinate x_
mat_on_vect(mat_row, POWER_, REZ);//Multiplication of the matrix mat_row by POWER
vector and obtain REZ vector
mat_on_vect(expo_, REZ, REZ_2);// Multiplication of matrix expo_ by vector REZ and
obtain vector REZ_2
for(int i=0;i<8;i++){
vector[i]=REZ_2[i]*delta_x;
}
}
//Solution of SLAE of 88 dimensionality by the Gauss method with discrimination of the
basic element
int GAUSS(double AA[8*11][8*11], double bb[8*11], double x[8*11]){
double A[8*11][8*11];
double b[8*11];
for(int i=0;i<(8*11);i++){
b[i]=bb[i];//we will operate with the vector of the и right parts to
provide that initial vector bb would not change when exiting the subprogram
for(int j=0;j<(8*11);j++){
A[i][j]=AA[i][j];//we will operate with A matrix to provide that
initial AA matrix would not change when exiting the subprogram
}
}
int e;//number of the row where main (maximal) coefficient in the column jj is
found
double s, t, main;//Ancillary quantity
for(int jj=0;jj<((8*11)-1);jj++){//Cycle by columns jj of transformation of A
matrix into upper-triangle one
e=-1; s=0.0; main=A[jj][jj];
for(int i=jj;i<(8*11);i++){//there is a number of у row where main
(maximal) element is placed in the column jj and row interchanging is made
if ((A[i][jj]*A[i][jj])>s) {//Instead of multiplication (potential
minus sign is deleted) it could be possible to use a function by abs() module
e=i; s=A[i][jj]*A[i][jj];
59
}
}
if (e<0) {
cout<<"Mistake "<<jj<<"\n"; return 0;
}
if (e>jj) {//If the main element isn’t placed in the row with jj number but
is placed in the row with у number
main=A[e][jj];
for(int j=0;j<(8*11);j++){//interchanging of two rows with e and jj
numbers
t=A[jj][j]; A[jj][j]=A[e][j]; A[e][j]=t;
}
t=b[jj]; b[jj]=b[e]; b[e]=t;
}
for(int i=(jj+1);i<(8*11);i++){//reduction to the upper-triangle matrix
for(int j=(jj+1);j<(8*11);j++){
A[i][j]=A[i][j]-(1/main)*A[jj][j]*A[i][jj];//re-calculation of
the coefficients of the row i>(jj+1)
}
b[i]=b[i]-(1/main)*b[jj]*A[i][jj];
A[i][jj]=0.0;//nullified elements of the row under diagonal element
of A matrix
}
}//Cycle by jj columns of transformation of A matrix into upper-triangle one
x[(8*11)-1]=b[(8*11)-1]/A[(8*11)-1][(8*11)-1];//initial determination of the last
element of the desired solution x (87th)
for(int i=((8*11)-2);i>=0;i--){//Computation of the elements of the solution x[i]
from 86th to 0th
t=0;
for(int j=1;j<((8*11)-i);j++){
t=t+A[i][i+j]*x[i+j];
}
x[i]=(1/A[i][i])*(b[i]-t);
}
return 0;
}
int main()
{
int nn;//Number of the harmonic starting from the 1st (without zero one)
int nn_last=50;//Number of the last harmonic
double Moment[100+1]={0};//An array of the physical parameter (momentum) that is
computed in each point between the boundaries
double step=0.05; //step=(L/R)/100 – step size of shell computation – a step of
integration interval (it should be over zero, i.e. be positive)
double h_div_R;//Value of h/R
h_div_R=1.0/100;
double c2;
c2=h_div_R*h_div_R/12;//Value of h*h/R/R/12
double nju;
nju=0.3;
double gamma;
gamma=3.14159265359/4;//The force distribution angle by the left boundary
//printing to files:
60
FILE *fp;
// Open for write
if( (fp = fopen( "C:/test.txt", "w" )) == NULL ) // C4996
printf( "The file 'C:/test.txt' was not opened\n" );
else
printf( "The file 'C:/test.txt' was opened\n" );
for(nn=1;nn<=nn_last;nn++){ //A CYCLE BY HARMONICS STARTING FROM THE 1st HARMONIC
(EXCEPT ZERO ONE)
double x=0.0;//A coordinate from the left boundary – it is needed in case of
heterogeneous system of the ODE for computing the particular vector FF
double expo_from_minus_step[8][8]={0};//The matrix for placement of the exponent
in it at the step of (0-x1) type
double expo_from_plus_step[8][8]={0};// The matrix for placement of the exponent
in it in the step of (x1-0) type
double mat_row_for_minus_expo[8][8]={0};//the auxiliary matrix for particular
vector computing when moving at step of (0-x1) type
double mat_row_for_plus_expo[8][8]={0};// the auxiliary matrix for particular
vector computing when moving at step of (x1-0) type
double U[4][8]={0};//The matrix of the boundary conditions of the left boundary of
the dimensionality 4x8
double u_[4]={0};//Dimensionality 4 vector
of the external influence for the
boundary conditions of the left boundary
double V[4][8]={0};//The boundary conditions matrix of the right boundary of the
dimensionality 4x8
double v_[4]={0};// The dimensionality 4 vector of the external influence for the
boundary conditions for the right boundary
double Y[100+1][8]={0};//The array of the vector-solutions of the corresponding
linear algebraic equations system (in each point of the interval between the boundaries):
MATRIXS*Y=VECTORS
double A[8][8]={0};//Matrix of the coefficients of the system of ODE
double FF[8]={0};//Vector of the particular solution of the heterogeneous ODE at
the integration interval sector
double Y_many[8*11]={0};// a composite vector consisting of the vectors Y(xi) in
11 points from point 0 (left boundary Y(0) to the point 10 (right boundary Y(x10))
double MATRIX_many[8*11][8*11]={0};//The matrix of the system of the linear
algebraic equations
double B_many[8*11]={0};// a vector of the right parts of the SLAE:
MATRIX_many*Y_many=B_many
double Y_vspom[8]={0};//an auxiliary vector
double Y_rezult[8]={0};//an auxiliary vector
double nn2,nn3,nn4,nn5,nn6,nn7,nn8;//Number of the
corresponding powers
nn2=nn*nn;
nn3=nn2*nn;
nn4=nn2*nn2;
nn5=nn4*nn;
nn
harmonic
nn6=nn4*nn2;
raised
to
nn7=nn6*nn;
nn8=nn4*nn4;
//Filling-in of non-zero elements of the A matrix of the coefficients of ODE
system
A[0][1]=1.0;
A[1][0]=(1-nju)/2*nn2; A[1][3]=-(1+nju)/2*nn; A[1][5]=-nju;
A[2][3]=1.0;
61
A[3][1]=(1+nju)/(1-nju)*nn; A[3][2]=2*nn2/(1-nju); A[3][4]=2*nn/(1-nju);
A[4][5]=1.0;
A[5][6]=1.0;
A[6][7]=1.0;
A[7][1]=-nju/c2; A[7][2]=-nn/c2; A[7][4]=-(nn4+1/c2); A[7][6]=2*nn2;
//Here firstly it is necessary to make filling-in by non-zero values of the matrix
and boundary conditions vector U*Y[0]=u_ (on the left) and V*Y[100]=v_ (on the right):
U[0][1]=1.0; U[0][2]=nn*nju; U[0][4]=nju; u_[0]=0.0;//Force T1 at the left
boundary is equal to zero
U[1][0]=-(1-nju)/2*nn; U[1][3]=(1-nju)/2; U[1][5]=(1-nju)*nn*c2; u_[1]=0.0;//Force
S* at the left boundary is equal to zero
U[2][4]=-nju*nn2; U[2][6]=1.0; u_[2]=0;//Momentum M1 at the left boundary is equal
to zero
U[3][5]=(2-nju)*nn2; U[3][7]=-1.0;
u_[3]=-sin(nn*gamma)/(nn*gamma);//Force Q1* at the left boundary is distributed at
the angle -gamma +gamma
V[0][0]=1.0; v_[0]=0.0;// The right boundary displacement u is equal to zero
V[1][2]=1.0; v_[1]=0.0;//The right boundary displacement v is equal to zero
V[2][4]=1.0; v_[2]=0.0;//The right boundary displacement w is equal to zero
V[3][5]=1.0; v_[3]=0.0;//The right boundary rotation angle is equal to zero
//Here initial filling-in of U*Y[0]=u_ и V*Y[100]=v_ terminates
exponent(A,(-step*10),expo_from_minus_step);//A negative step (step value is less
than zero due to direction of matrix exponent computation)
//x=0.0;//the initial value of coordinate – for partial vector computation
//mat_row_for_partial_vector(A, step, mat_row_for_minus_expo);
//Filling-in of the SLAE coefficients matrix MATRIX_many
for(int i=0;i<4;i++){
for(int j=0;j<8;j++){
MATRIX_many[i][j]=U[i][j];
MATRIX_many[8*11-4+i][8*11-8+j]=V[i][j];
}
B_many[i]=u_[i];
B_many[8*11-4+i]=v_[i];
}
for(int kk=0;kk<(11-1);kk++){//(11-1) unit matrix and EXPO matrix should be
written into MATRIX_many
for(int i=0;i<8;i++){
MATRIX_many[i+4+kk*8][i+kk*8]=1.0;//filling-in by unit matrixes
for(int j=0;j<8;j++){
MATRIX_many[i+4+kk*8][j+8+kk*8]=expo_from_minus_step[i][j];//filling-in by matrix exponents
}
}
}
//Solution of the system of linear algebraic equations
GAUSS(MATRIX_many,B_many,Y_many);
//Computation of the state vectors in 101 points – left point 0 and right point
100
exponent(A,step,expo_from_plus_step);
for(int i=0;i<11;i++){//passing points filling-in in all 10 intervals (we will
obtain points from 0 to 100) between 11 nodes
for(int j=0;j<8;j++){
62
Y[0+i*10][j]=Y_many[j+i*8];//in 11 nodes the vectors are taken from
SLAE solutions – from Y many
}
}
for(int i=0;i<10;i++){//passing points filling-in in 10 intervals
for(int j=0;j<8;j++){
Y_vspom[j]=Y[0+i*10][j];//the initial vector for ith segment, zero
point, starting point of the ith segment
}
mat_on_vect(expo_from_plus_step, Y_vspom, Y_rezult);
for(int j=0;j<8;j++){
Y[0+i*10+1][j]=Y_rezult[j];//the first point of the interval fillingin
Y_vspom[j]=Y_rezult[j];//for the next step
}
mat_on_vect(expo_from_plus_step, Y_vspom, Y_rezult);
for(int j=0;j<8;j++){
Y[0+i*10+2][j]=Y_rezult[j];//filling-in of the 2nd point of the
interval
Y_vspom[j]=Y_rezult[j];//for the next step
}
mat_on_vect(expo_from_plus_step, Y_vspom, Y_rezult);
for(int j=0;j<8;j++){
Y[0+i*10+3][j]=Y_rezult[j];//filling-in of the 3rd point of the
interval
Y_vspom[j]=Y_rezult[j];//for the next step
}
mat_on_vect(expo_from_plus_step, Y_vspom, Y_rezult);
for(int j=0;j<8;j++){
Y[0+i*10+4][j]=Y_rezult[j];//filing-in of the 4th point of the
interval
Y_vspom[j]=Y_rezult[j];//for the next step
}
mat_on_vect(expo_from_plus_step, Y_vspom, Y_rezult);
for(int j=0;j<8;j++){
Y[0+i*10+5][j]=Y_rezult[j];// filing-in of the 5th point of the
interval
Y_vspom[j]=Y_rezult[j];// for the next step
}
mat_on_vect(expo_from_plus_step, Y_vspom, Y_rezult);
for(int j=0;j<8;j++){
Y[0+i*10+6][j]=Y_rezult[j];// filing-in of the 6th point of the
interval
Y_vspom[j]=Y_rezult[j];// for the next step
}
mat_on_vect(expo_from_plus_step, Y_vspom, Y_rezult);
for(int j=0;j<8;j++){
Y[0+i*10+7][j]=Y_rezult[j];// filing-in of the 7th point of the
interval
Y_vspom[j]=Y_rezult[j];// for the next step
}
mat_on_vect(expo_from_plus_step, Y_vspom, Y_rezult);
for(int j=0;j<8;j++){
63
Y[0+i*10+8][j]=Y_rezult[j];// filing-in of the 8th point of the
interval
Y_vspom[j]=Y_rezult[j];// for the next step
}
mat_on_vect(expo_from_plus_step, Y_vspom, Y_rezult);
for(int j=0;j<8;j++){
Y[0+i*10+9][j]=Y_rezult[j];// filing-in of the 9th point of the
interval
Y_vspom[j]=Y_rezult[j];// for the next step
}
}
//Computation of the momentum in all points between the boundaries
for(int ii=0;ii<=100;ii++){
Moment[ii]+=Y[ii][4]*(-nju*nn2)+Y[ii][6]*1.0;//Momentum M1 in the point
[ii]
//U[2][4]=-nju*nn2; U[2][6]=1.0; u_[2]=0;//Momentum M1
}
}//CYCLE BY HARMONICS TERMINATES HERE
for(int ii=0;ii<=100;ii++){
fprintf(fp,"%f\n",Moment[ii]);
}
fclose(fp);
printf( "PRESS any key to continue...\n" );
_getch();
return 0;
}
64
List of published works.
1. Виноградов А.Ю. Вычисление начальных векторов для численного решения краевых
задач // - Деп. в ВИНИТИ, 1994. -N2073- В94. -15 с.
2. Виноградов А.Ю., Виноградов Ю.И. Совершенствование метода прогонки Годунова
для задач строительной механики // Изв. РАН Механика твердого тела, 1994. -№4. -С. 187191.
3. Виноградов А.Ю. Вычисление начальных векторов для численного решения краевых
задач // Журнал вычислительной математики и математической физики, 1995. -Т.З5. -№1. С. 156-159.
4. Виноградов А.Ю. Численное моделирование произвольных краевых условий для
задач строительной механики тонкостенных конструкций // Тез. докладов Белорусского
Конгресса по теоретической и прикладной механике "Механика-95". Минск, 6-11 февраля
1995 г., Гомель: Изд- во ИММС АНБ, 1995. -С63-64.
5. Vinogradov A.Yu. Numerical modeling of boundary conditions in deformation problems of
structured material in thin wall constructions // International Symposium "Advances in
Structured and Heterogeneous Continua II". Book Absfracts. August, 14-16, 1995, Moscow,
Russia. -P.51.
6. Виноградов А.Ю. Численное моделирование краевых условий в задачах
деформирования тонкостенных конструкций из композиционных материалов // Механика
Композиционных Материалов и Конструкций, 1995. -T.I. -N2. - С. 139-143.
7. Виноградов А.Ю. Приведение краевых задач механики элементов приборных
устройств к задачам Коши для выбранной точки // Прикладная механика в приборных
устройствах. Меж вуз. сб. научных трудов. -Москва: МИРЭА, 1996.
8. Виноградов А.Ю. Модификация метода Годунова// Труды Международной научнотехнической конференции "Современные проблемы машиноведения", Гомель: ГПИ им.
П.О. Сухого, 1996. - С.39-41.
9. Виноградов А.Ю., Виноградов Ю.И. Методы переноса сложных краевых условий для
жестких дифференциальных уравнений строительной механики //
Труды
Международной конференции «Ракетно-космическая техника: фундаментальные
проблемы механики и теплообмена», Москва, 1998.
10. Виноградов А.Ю. Метод решения краевых задач путем переноса условий с краев
интервала интегрирования в произвольную точку // Тез. докладов Международной
конференции "Актуальные проблемы механики оболочек", Казань, 2000.-С. 176.
11. Виноградов Ю.И., Виноградов А.Ю., Гусев Ю.А. Метод переноса краевых
условий для дифференциальных уравнений
теории
оболочек //
Труды
Международной
конференции "Актуальные проблемы механики оболочек", Казань,
2000. -С. 128-132.
12. Виноградов А.Ю., Виноградов Ю.И. Метод переноса краевых условий функциями
Коши-Крылова для жестких линейных обыкновенных дифференциальных уравнений. //
ДАН РФ, – М.: 2000, т. 373, №4, с. 474-476.
13. Виноградов А.Ю., Виноградов Ю.И. Функции Коши-Крылова и алгоритмы решения
краевых задач теории оболочек // ДАН РФ, - М.: 2000. -Т.375.-№3.-С. 331-333.
14. Виноградов А.Ю. Численные методы переноса краевых условий // Журнал
"Математическое
моделирование",
изд-во
РАН,
Институт
математического
моделирования, - М.: 2000, Т. 12 , № 7, с. 3-6.
15. Виноградов А.Ю., Гусев Ю.А. Перенос краевых условий функциями Коши-Крылова в
задачах строительной механики // Тез. докладов VIII Всероссийского съезда по
теоретической и прикладной механике, Пермь, 2001. -С. 109-110.
65
16. Виноградов Ю.И., Виноградов А.Ю. Перенос краевых условий функциями КошиКрылова // Тез. докладов Международной конференции "Dynamical System Modeling and
Stability Investigation"- "DSMSI-2001", Киев, 2001.
17. Виноградов А.Ю., Виноградов Ю.И., Гусев Ю.А, Клюев Ю.И. Перенос краевых
условий функциями Коши-Крылова и его свойства // Изв. РАН МТТ, №2. 2001. -С.155161.
18. Виноградов Ю.И., Виноградов А.Ю., Гусев Ю.А.Численный метод переноса краевых
условий для жестких дифференциальных уравнений строительной механики // Журнал
"Математическое
моделирование",
изд-во
РАН,
Институт
математического
моделирования, - М.: 2002, Т. 14, №9, с.3-8.
19. Виноградов Ю.И., Виноградов А.Ю. Простейший метод решения жестких краевых
задач // Фундаментальные исследования. – 2014. – № 12–12. – С. 2569-2574;
20. Виноградов Ю.И., Виноградов А.Ю. Решение жестких краевых задач строительной
механики (расчет оболочек составных и со шпангоутами) методом Виноградовых (без
ортонормирования) // Современные проблемы науки и образования. – 2015. - №1;
21. Виноградов Ю.И., Виноградов А.Ю. Уточненный метод Виноградовых переноса
краевых условий в произвольную точку интервала интегрирования для решения жестких
краевых задач // Современные наукоемкие технологии. – 2016. – № 11-2. – С. 226-232;
22. Виноградов А.Ю. Методы решения жестких и нежестких краевых задач: монография //
Волгоград: Изд-во ВолГУ, 2016. – 128 с.
23. Виноградов А.Ю. Программа (код) на С++ решения жесткой краевой задачи методом
А.Ю. Виноградова // exponenta.ru/educat/systemat/vinogradov/index.asp
24. Виноградов А.Ю. Методы А.Ю.Виноградова решения краевых задач, в том числе
жестких краевых задач // exponenta.ru/educat/systemat/vinogradov/index2.asp
25. Виноградов А.Ю. Расчет составных оболочек и оболочек со шпангоутами простейшим
методом решения жестких краевых задач (без ортонормирования) //
exponenta.ru/educat/systemat/vinogradov/index3.asp
26. Виноградов А.Ю. Методы решения жестких и нежестких краевых задач – народная
докторская // exponenta.ru/educat/systemat/vinogradov/index5.asp
27. Виноградов А.Ю. Методы решения жестких и нежестких краевых задач (9 методов) //
exponenta.ru/educat/systemat/vinogradov/index6.asp
Отзывы:
Авторизуйтесь, чтобы оставить отзыв