The reverse procedure of converting a system of linear ODEs into a single equation is discussed in section viii of Part 2. Not every system of n linear first order differential equations can be transfered to a single n-th order differential equation.

Given a single ordinary differential equation, one method of finding numerical solutions entails transferring it into an equivalent system of differential equations of the first order. When a single differential equation has an isolated highest derivative, it is always possible to transfer the differential equation into an equivalent system of differential equations of the first order. This can be achieved by denoting sequential derivatives of the unknown variable with a new dependent variable with the exception of the last derivative, which is used to incorporate the given single differential equation. Other options are also possible for such a conversion, which will be clear from the examples.

The word "equivalent" indicates that both a single differential equation and the corresponding system of equations have the same set of solutions (so no solution is lost or added). The opposite is not always true, and some systems of n ordinary differential equations cannot be reduced to an equivalent single n-th order differential equation. The general solution of the single n-th order differential equation depends on n arbitrary constants and it belongs to the space of functions having n continuous derivatives, which is usually denoted by Cⁿ. A system of n differential equations has the same properties and its general solution also depends on n arbitrary constants. Therefore, a single n-th order ODE can be equivalent to the system of (n+1) first order ODEs only when the latter does not require extra differentiation of the unknown function. However, a single n-th order ODE cannot be equivalent to the (n−1) system of first order ODEs because the latter contains (n−1) arbitrary constants and requires utilization of C^n-1. A key illustration highlighting this is the pendulum equation, presented in Example 5, which does not use the derivative operator, and where a single second order equation is transformed into four first order equations representing an equivalent system. On the other hand, the derivative operator is used in Example 1, showing that the resulting system of equations is not equivalent to the original equation. Therefore, the key takeaway is to avoid using the derivative operator because it is an unbounded operation.

There are a couple of reasons why it is convenient to transfer a single differential equation of an order higher than one to a system of differential equations of first order. First of all, it is simpler to analyze theoretically the first order vector differential equation than deal with a higher order differential equation because the latter can be naturally included into a broader topic: systems of differential equations. Second, it is easier to use and implement a universal numerical algorithm for a first order derivative operator instead of utilizing special numerical procedures for higher order derivatives. It is more accurate to monitor errors of numerical calculations for a system of first order differential equations than for a single equation. Also, there are many situations in which we not only want to know the solution to an ODE (acronym for Ordinary Differential Equations), but also the derivative or acceleration of the solution. This information can be naturally extracted from the solution to the system of first order differential equations (with no extra work).

When treating the independent variable as a time, it is common to use Newton's notation \( \dot{y} \) for the derivative instead of Lagrange notation y' or Leibniz's notation \( \displaystyle {\text d}y/{\text d}t . \) Sometimes, we also utilize Euler's notation for the derivative operator: \( \displaystyle \texttt{D} = {\text d}/{\text d}t . \)

Example 1: Linear second order differential equation

Example 1: Let us start with a simple case of a constant coefficient linear differential equation when all manipulations become transparent. So we consider the following initial value problem for a linear second order differential equation

\begin{equation*} %\label{EqConv.1} x''(t) + 2x'(t) - 3x(t) = 0, \qquad x(0) =3, \quad x' (0) =-1. \tag{1.1} \end{equation*}

This is a homogeneous, linear, second-order ODE with constant coefficients. As a brief refresher from the previous course, recall that its solution can be computed by finding the roots of the corresponding characteristic equation:

\[ \lambda^2 + 2\,\lambda - 3 = 0. \tag{1.2} \]

There are two distinct, real roots to the quadratic (characteristic) equation: λ₁ = 1 and λ₂ = -3, and the general solution is a linear combination of two exponent functions:

\[ x(t) = c_{1}e^{\lambda_{1}t} + c_{2}e^{\lambda_{2}t} = c_{1}e^{t} + c_{2}e^{-3t} , \]

where c₁ and c₂ are arbitrary constants. To satisfy the initial conditions, we have to choose these constants so that x(0) = c₁ + c₂ = 3 and x'(0) = c₁ -3 c₂ = -1. This yields c₁ = 2 and c₂ = 1; and the solution becomes

\begin{equation*} %\label{EqConv.2} x(t) = 2\, e^t + e^{-3t} \qquad \Longrightarrow \qquad x'(t) = 2\, e^t -3\, e^{-3t} . \tag{1.3} \end{equation*}

Click to see matlab code

syms x(t) t

% derivative
Dx = diff(x);

eqn = diff(x, t, 2) == - 2*Dx + 3*x;

cond1 = x(0) == 3;
cond2 = Dx(0) == -1;
conds = [cond1 cond2];
 
soln(t) = dsolve(eqn, conds)
soln = simplify(soln)
end

Copy

Now we convert the given problem to a vector equation by letting x₁ = x and its first derivative we denote by x₂. Note that second derivative is known from the given equation \( x'' = 3\,x(t) - 2\, x' (t) = 3\, x_1 - 2\, x_2 . \) This leads to the system of first order differential equations

\[ \begin{split} x'_1 &= x_2 = x', \\ x'_2 &= x'' = -2\, x_2 + 3\,x_1 . \end{split} \]

Upon introducing a column vector x = [ x₁, x₂]^T (here "T" stands for transposition; see first chapter of this tutorial), we can rewrite the system of ODEs in more appropriate vector form:

\begin{equation*} %\label{EqConv.3} \frac{\text d}{{\text d}t} \begin{bmatrix} x_1 (t) \\ x_2 (t) \end{bmatrix} = \begin{bmatrix} 0& \phantom{-}1 \\ 3 & -2 \end{bmatrix} \begin{bmatrix} x_1 (t) \\ x_2 (t) \end{bmatrix} , \qquad \begin{bmatrix} x_1 (0) \\ x_2 (0) \end{bmatrix} = \begin{bmatrix} \phantom{-}3 \\ -1 \end{bmatrix} . \tag{1.4} \end{equation*}

Click to see matlab code

syms x1(t) x2(t) 
 
eqn1 = diff(x1, t) == x2;
eqn2 = diff(x2, t) == -2*x2 + 3*x1;
eqns = [eqn1, eqn2];

 
cond1 = x1(0) == 3;
cond2 = x2(0) == -1;
conds = [cond1 cond2];

 
[soln1(t), soln2(t)] = dsolve(eqns, conds);
soln = [soln1(t), soln2(t)];
soln = simplify(soln)

Copy

With matrix notation,

\[ {\bf x} (t) = \begin{bmatrix} x_1 (t) \\ x_2 (t) \end{bmatrix} , \qquad {\bf A} = \begin{bmatrix} 0& \phantom{-}1 \\ 3 & -2 \end{bmatrix} , \qquad {\bf x}_0 = \begin{bmatrix} \phantom{-}3 \\ -1 \end{bmatrix} , \]

we rewrite the system of ODEs and corresponding initial conditions in a more appropriate vector form:

\[ \frac{\text d}{{\text d}t}\,{\bf x}(t) = {\bf A}\,{\bf x}(t) , \qquad {\bf x}(0) = {\bf x}_0 . \]

In order to find the solution to this vector differential equation, we first need to find the eigenvalues of the corresponding matrix (which is called the companion matrix):

\[ \det \left( \lambda {\bf I} - {\bf A} \right) =0 \qquad \Longleftrightarrow \qquad \det \left( \lambda \begin{bmatrix} 1&0 \\ 0&1 \end{bmatrix} - \begin{bmatrix} 0&\phantom{-}1 \\ 3&-2 \end{bmatrix} \right) = \det \begin{bmatrix} \lambda & -1 \\ -3& \lambda +2 \end{bmatrix} = \lambda^2 + 2\,\lambda - 3 =0 , \]

where I is the identity matrix. From observation, the eigenvalues of matrix A are the same as the roots of the characteristic equation λ² + 2 λ −3 = 0 for the given second order differential equation.

Sometimes, we need to monitor the residual of numerical calculations by introducing an additional auxiliary variable y₃ = x'', the second derivative of the unknown solution x(t), and keep other variables: y₁ = x and y₂ = x'. Then its derivative, according to the given differential equation, will be:

\[ y'_3 = x''' = \frac{\text d}{{\text d}t} \left( x'' \right) = \frac{\text d}{{\text d}t} \left( -2\,y_2 + 3\, y_1 \right) = -2\,y'_2 + 3\, y'_1 = -2 \left( -2\, y_2 + 3\, y_1 \right) + 3\, y_2 = 7\, y_2 -6\, y_1 . \]

This allows us to compose the system of three equations. There are two options to achieve it. First, we recall that the derivative of y₂ is our new dependent variable y₃. So we get

\begin{equation*} %\begin{split} \begin{cases} y'_1 &= y_2 = x' , \\ y'_2 &= y_3 = x'' , \\ y'_3 &= -2\, y_3 +3 \, y_2 = x''' ; \end{cases} %\end{split} \qquad\quad %\begin{split} \begin{cases} y_1 (0) &= 3 \\ y_2 (0) &= -1 \\ y_3 (0) &= 11 . \end{cases} %\end{split} \tag{1.5} \end{equation*}

Therefore, we obtain a vector differential equation \( \dot{\bf y} = {\bf C}\,{\bf y}(t) \) that is equivalent (we hope?) to the the original IVP (1.1)

\[ \frac{\text d}{{\text d}t} \begin{bmatrix} y_1 (t) \\ y_2 (t) \\ y_3 (t) \end{bmatrix} = \begin{bmatrix} 0&1&\phantom{-}0 \\ 0&0&\phantom{-}1 \\ 0&3&-2 \end{bmatrix} \begin{bmatrix} y_1 (t) \\ y_2 (t) \\ y_3 (t) \end{bmatrix} , \qquad \begin{bmatrix} y_1 (0) \\ y_2 (0) \\ y_3 (0) \end{bmatrix} =\begin{bmatrix} \phantom{1}3 \\ -1 \\ 11 \end{bmatrix} , \qquad\mbox{with}\quad {\bf C} = \begin{bmatrix} 0&1&\phantom{-}0 \\ 0&0&\phantom{-}1 \\ 0&3&-2 \end{bmatrix} . \]

The initial condition for this function also follows from the given differential equation (1.1):

\[ y_3 = x'' = 3\, x(t) - 2\, x' (t) \qquad \Longrightarrow \qquad y_3 (0) = 3\,x(0) -2\, x' (0) = 3\cdot 3 -2 \cdot (-1) = 11 . \]

Note that the derivative of y₃ is expressed through the original function x(t) as \( y'_3 = 7\, y_2 -6\, y_1 . \) So another option is to rewrite the given system of ODEs as

\begin{equation*} %\label{EqConv.4} %\begin{split} \begin{cases} y'_1 &= y_2 = x' , \\ y'_2 &= y_3 = x'' , \\ y'_3 &= -6\, y_1 +7\, y_2 = x''' ; \end{cases} %\end{split} \qquad\quad %\begin{split} \begin{cases} y_1 (0) &= 3 \\ y_2 (0) &= -1 \\ y_3 (0) &= 11 . \end{cases} %\end{split} \tag{1.6} \end{equation*}

In vector form, it becomes:

\[ \frac{\text d}{{\text d}t} \begin{bmatrix} y_1 (t) \\ y_2 (t) \\ y_3 (t) \end{bmatrix} = \begin{bmatrix} \phantom{-}0&1&0 \\ \phantom{-}0&0&1 \\ -6&7&0 \end{bmatrix} \begin{bmatrix} y_1 (t) \\ y_2 (t) \\ y_3 (t) \end{bmatrix} , \qquad \begin{bmatrix} y_1 (0) \\ y_2 (0) \\ y_3 (0) \end{bmatrix} =\begin{bmatrix} \phantom{1}3 \\ -1 \\ 11 \end{bmatrix} , \qquad\mbox{with}\quad {\bf B} = \begin{bmatrix} \phantom{-}0&1&0 \\ \phantom{-}0&0&1 \\ -6&7&0 \end{bmatrix} . \]

DSolve[{y1'[t] == y2[t], y2'[t] == -2*y2[t] + 3*y1[t], y3'[t] == 7*y2[t] - 6*y1[t], y1[0] == 3, y2[0] == -1, y3[0] == 11}, {y1[t], y2[t], y3[t]}, t]

{{y1[t] -> E^(-3 t) (1 + 2 E^(4 t)), y2[t] -> E^(-3 t) (-3 + 2 E^(4 t)), y3[t] -> E^(-3 t) (9 + 2 E^(4 t))}}

Click to see matlab code

syms y1(t) y2(t) y3(t)

eqn1 = diff(y1, t) == y2;
eqn2 = diff(y2, t) == -2*y2 + 3*y1;
eqn3 = diff(y3, t) == 7*y2 - 6*y1;
eqns = [eqn1, eqn2, eqn3];

cond1 = y1(0) == 3;
cond2 = y2(0) == -1;
cond3 = y3(0) == 11;
conds = [cond1 cond2 cond3];
 
[soln1(t), soln2(t), soln3(t)] = dsolve(eqns, conds);
    soln = [soln1(t), soln2(t), soln3(t)];
    soln = simplify(soln)

Copy

Despite that the first two components (y₁ and y₂) of the solution to the initial value problem (1.6) are the same as (1.3),

\[ \begin{bmatrix} y_1 (t) \\ y_2 (t) \\ y_3 (t) \end{bmatrix} = \begin{bmatrix} e^{-3t} + 2\,e^{t} \\ -3\,e^{-3t} + 2\,e^{t} \\ 9\,e^{-3t} + 2\,e^{t} \end{bmatrix} , \]

the system of differential equations (1.6) is not equivalent to (1.1) because the corresponding matrix B has three eigenvalues {−3, 2, 1}. Recall that matrix A has two eigenvalues {−3, 1}. It becomes clear when we rewrite the system of three equations with matrix B in the following form with Euler's notation \( \texttt{D} = {\text d}/{\text d}t \) for derivative operator:

\[ \begin{split} \texttt{D} y_1 &= y_2 , \\ \texttt{D} y_2 &= y_3 , \\ \texttt{D} y_3 &= -6\,y_1 + 7\,y_2 . \end{split} \]

The system above contains three unbounded derivative operators, but the original differential equation (1.1) is of order two. So its general solution is

\[ y_1 (t) = x(t) = C_1 e^t + C_2 e^{-3t} + C_3 e^{2t} , \]

with some arbitrary constants C₁, C₂, and C₃. Upon using the appropriate initial conditions we eliminated unwanted additional constant C₃ from its general solution. So we obtain the same solution from the system of three equations as the original IVP, at least from theoretical point of view. However, matrix B has a dominant eigenvalue λ = 2>1 and the corresponding numerical algorithm should avoid it.

Then the difference of y₃ and \( 3\, x(t) - 2\,x' (t) = 3\,y_1 - 2\, y_2 \) gives the residual (error) of calculations.

We use the standard Mathematica command NDSolve to solve numerically the system of three equations (1.6) and plot the residual of calculations for the equation. Blue curve gives the error with the true second derivative, and the red curve gives the residual based on numerical calculation.

sol3 = NDSolve[{y1'[t] == y2[t], y2'[t] == y3[t], y3'[t] == -6*y1[t] + 7*y2[t], y1[0] == 3, y2[0] == -1, y3[0] == 11}, {y1, y2, y3}, {t, 0, 3}];;
sol2 = NDSolve[{x1'[t] == x2[t], x2'[t] == 3*x1[t] - 2*x2[t], x1[0] == 3, x2[0] == -1}, {x1, x2}, {t, 0, 3}];
a1 = Plot[ Evaluate[y3[t] /. sol3] + Evaluate[(2*x2[t] - 3*x1[t]) /. sol2], {t, 0, 1}, PlotTheme -> "Web"]
a2 = Plot[9 E^(-3 t) + 2 E^t - Evaluate[(y3[t]) /. sol3], {t, 0, 1}, PlotTheme -> "Business"]

Figure 1. Residuals of the second derivative: the error between the second derivative y3(t) and the true value in blue.		Figure 2. The residual of the second derivative y3(t) = d²x/dt² and the numerical solution 3x2(t) - 2x1(t).

Adjust coefficients:

odefcn2 = @(t,y)[y(2);y(3);((31/4)*y(2)) - ((15/4)*y(1))];
y0 = [ 3; -2 ; 9.5]; 
t = [0,3] ; 
sol = ode45(odefcn2, t , y0); 
plot(sol.x,((2*sol.y(3,:))-(5*sol.y(2,:))+(3*sol.y(1,:))))

Note that introducing the system of three equations with artificial differentiation does not provide an equivalent system for the given equation because the derivative operator is unbounded.

The transformation to a system of first order ODEs is straight forward, but it is not unique. Let us try another set of new dependent variables:

\[ z_1 = x, \quad z_2 = \dot{x} - x \qquad\Longrightarrow \qquad \dot{z}_2 = \ddot{x} - \dot{x} = -3 \left( \dot{x} - x \right) . \]

This leads to the system of first order differential equations:

\[ \begin{split} \dot{z}_1 &= z_1 + z_2 \qquad \Longleftrightarrow \qquad \dot{x} = \dot{x} , \\ \dot{z}_2 &= -3\, z_2 . \end{split} \]

We can also rewrite the above system in vector form

\[ \frac{\text d}{{\text d}t} \begin{bmatrix} z_1 (t) \\ z_2 (t) \end{bmatrix} = \begin{bmatrix} 1& 1 \\ 0 & -3 \end{bmatrix} \begin{bmatrix} z_1 (t) \\ z_2 (t) \end{bmatrix} , \qquad \begin{bmatrix} z_1 (0) \\ z_2 (0) \end{bmatrix} = \begin{bmatrix} 3 \\ -4 \end{bmatrix} , \]

because \( z_2 (0) = \dot{x}(0) - x(0) = -1 -3 =-4 . \) The characteristic polynomial for the corresponding matrix is

\[ \det \left( \lambda \begin{bmatrix} 1&0 \\ 0&1 \end{bmatrix} - \begin{bmatrix} 1&1 \\ 0&-3 \end{bmatrix} \right) = \det \begin{bmatrix} \lambda -1& -1 \\ 0& \lambda +3 \end{bmatrix} = \lambda^2 + 2\,\lambda - 3 =0 . \]

   ■

function Converting_1A
 syms A B
 eqn1 = -4*A + 2*B == 0;
 eqn2 = -4*B - 2*A == 10;
 eqns = [eqn1, eqn2];
 
 cond1 = A;
 cond2 = B;
 conds = [cond1 cond2];
 
 [soln1, soln2] = solve(eqns, conds);
 % A = soln1, and B = soln2
 soln = [soln1, soln2];
 soln = simplify(soln)
end

 syms c1 c2
 eqn1 = c1 + c2 == 4;
 eqn2 = c1 - 3*c2 == 1;
 eqns = [eqn1, eqn2];
 
 cond1 = c1;
 cond2 = c2;
 conds = [cond1 cond2];
 [soln1, soln2] = solve(eqns, conds);
 % c1 = soln1, c2 = soln2
 soln = [soln1, soln2];
 soln = simplify(soln)

syms x(t) t

% derivative
Dx = diff(x);

eqn = diff(x, t, 2) == 10*sin(t) + 3*x - 2*Dx;

cond1 = x(0) == 3;
cond2 = Dx(0) == -1;
conds = [cond1 cond2];

soln(t) = dsolve(eqn, conds);
soln = simplify(soln)
 
end

Example 2: Mathieu equation

Example 2: With our next example, we show that variable coefficients do not prevent conversion of a single equation of the second order into a system of two equations of the first order. The differential equation for modeling non-linear electrical oscillations has the form

\[ \texttt{D}^2 x(t) + q(t)\, x(t) = 0 , \tag{2.1} \]

where \( \ddot{x} = \frac{{\text d}^2 x}{{\text d} t^2} = \texttt{D}^2 x(t) \) stands for the second derivative with respect to time variable t, and q(t) is a single-valued periodic function. This linear differential equation has practical importance in the analysis of time-variable circuits. If function q(t) is a periodic function, then Eq.(2.1) is called Hill's equation. It was discovered by an American astronomer and mathematician from Rutgers University George William Hill (1838--1914) in 1886. Hill's differential equation is widely applicable to many problems involving sawtooth variations, rectangular ripple analysis, exponential variations, etc.

A particular case of Hill's differential equation when the function q(t) has the form

\[ q(t) = a_0 - a_1 \cos \omega t \tag{2.2} \]

is known as Mathieu's equation. The problem of vibrations of an elliptical membrane was first seriously studied by the French mathematician Émile Léonard Mathieu (1835--1890) in 1868, in his celebrated paper, “Mémoire sur le mouvement vibratoire d’une membrane de forme elliptique”, J. Math. Pures Appl. 13, 137 (1868). (http://eudml.org/doc/234720). Mathieu proved that the vertical displacement produced from oscillations on an elliptical drumhead is governed, in general, by the 2D Helmholtz equation. Furthermore, Mathieu showed that when the 2D Helmholtz equation is transformed into elliptical coordinates, the differential equations, respectively, for the radial and angular coordinates of the elliptical system become separable. The result is two almost identical linear second-order differential equations, each of which is aptly named the “Mathieu equation”. So, to this end, Mathieu’s equation has the canonical form

\[ \texttt{D}^2 x(t) + \left( \lambda - 2q\,\cos 2t \right) x(t) = 0 , \tag{2.3} \]

where the separation constant λ and the physical parameter q are both real.

Consider the initial value problem for the linear homogeneous equation

\[ \ddot{x} + \left( A - \sin 2t \right) x =0 , \qquad x(0) =1, \quad \dot{x} (0) =0 , \]

where A is a parameter. Here double superdot \( \ddot{x} = \frac{{\text d}^2 x}{{\text d} t^2} \) stands for the second derivative with respect to time variable t. Upon converting the given differential equation to the system of first order equations

\[ \begin{cases} \dot{x} =y , \\ \dot{y} = - \left( A - \sin 2t \right) x , \end{cases} \qquad x(0) =1, \quad y(0) =0 , \]

 

So we see the behavior of solutions depends on the value of parameter A: when it is greater than 1, solutions have different properties. Therefore, such value of A = 1 is called the bifurcation point (or node). Mathematica has a dedicated command NDSolve to find an a approximate solution to the initial value problem. Correspondingly, we use it to plot solutions of the given initial value problem for different values of parameter A. Click to see matlab code

A = 0.5;
tspan = [0 30];
y0 = [1 0];
[t,y] = ode45(@(t,y) odefcn(t,y,A), tspan, y0);
figure
plot(t,y(:,1),'-o')

A=1;
tspan = [0 30];
y0 = [1 0];
[t,y] = ode45(@(t,y) odefcn(t,y,A), tspan, y0);
figure
plot(t,y(:,1),'-o')

A=2;
tspan = [0 30];
y0 = [1 0];
[t,y] = ode45(@(t,y) odefcn(t,y,A), tspan, y0);
figure
plot(t,y(:,1),'-o')

function dydt = odefcn(t,y,A)
  dydt = zeros(2,1);
  dydt(1) = y(2);
  dydt(2) = -1*(A-sin(2*t))*y(1);
end

Copy

Solution for A = 1/2

Solution for A = 1

Solution for A = 2

All three solutions plotted together

When you give a setting for WorkingPrecision (which is an option for various numerical operations that specifies how many digits of precision should be maintained in internal computations), this typically defines an upper limit on the precision of the results from a computation. But within this constraint you can tell the Wolfram Language how much precision and accuracy you want to get. You should realize that for many kinds of numerical operations, increasing precision and accuracy goals by only a few digits can greatly increase the computation time required. Given a particular setting for WorkingPrecision, each of the functions for numerical operations in the Wolfram Language uses certain default settings for PrecisionGoal (which is an option for various numerical operations that specifies how many effective digits of precision should be sought in the final result) and AccuracyGoal (which is an option for various numerical operations which specifies how many effective digits of accuracy should be sought in the final result). Typical is the case of NDSolve, in which these default settings are equal to half the settings given for WorkingPrecision.

One can try to control the accuracy of numerical calculations by introducing the residual variable: \( z (t) = \ddot{x} = - \left( A - \sin 2t \right) x = \left( \sin 2t - A \right) x . \) Upon its differentiation, we get another system of differential equations with three unknowns:

\[ \begin{cases} \dot{x} =y , \\ \dot{y} = z , \\ \dot{z} = 2x\,\cos 2t + \left( \sin 2t - A \right) y , \end{cases} \qquad x(0) =1, \quad y(0) (0) =0 , \quad z(0) = -A. \]

This system of 3 differential equations is not equivalent to the original Mathieu equation because its general solution depends on 3 arbitrary constants and its existence requires a solution to be from the space C³. To eliminate extra constant, we need to impose the third initial condition, so the initial value problem will have the same solution as the original IVP for the Mathieu equation.

⁎ ✱ ✲ ✳ ✺ ✻ ✼ ✽ ❋
◾

   

Now we generalize the previous examples by considering the general linear differential equation. Suppose we are given an n-th order linear differential equation

\begin{equation} \label{EqConv.5} x^{(n)} (t) + a_{n-1} x^{(n-1)} (t) + \cdots + a_1 x' (t) + a_0 x(t) = f(t) , \end{equation}

where all coefficients 𝑎_i and function f are known. Upon introducing n new dependent variables

\[ y_1 (t) = x(t), \quad y_2 (t) = y'_1 (t) = x' (t), \quad y_3 (t) = y'_2 (t) = x'' (t), \quad \ldots , \quad y_n (t) = y'_{n-1} (t) = x^{(n-1)} (t) , \]

we reduce the given n-th order linear differential equation \eqref{EqConv.5} to an equivalent system of the first order equations.

We can assign to equation \eqref{EqConv.5} the differential operator \( L\left[ t, \texttt{D} \right] = \texttt{D}^n + a_{n-1} (t)\, \texttt{D}^{n-1} + \cdots + a_1 (t)\,\texttt{D} + a_0 (x)\,\texttt{I} , \) where \( \texttt{D} = {\text d}/{\text d}t \) and \( \texttt{D}^0 = \texttt{I} \) are the derivative operator and the identical operator, respectively. Then we rewrite the differential equation in a more concise form:

\[ L\left[ t, \texttt{D} \right] x(t) = f(t) . \]

When coefficients in this differential operator of order n do not depend on t, we obtain a constant coefficient differential operator \( L\left[ \texttt{D} \right] = \texttt{D}^n + a_{n-1} \, \texttt{D}^{n-1} + \cdots + a_1 \,\texttt{D} + a_0 \,\texttt{I}, \) to which we can assign a polynomial

\[ p(\lambda ) = \lambda^n + a_{n-1} \lambda^{n-1} + \cdots + a_1 \lambda + a_0 , \]

called the characteristic polynomial. Upon introducing n new dependent variables

\[ y_1 (t) = x(t), \quad y_2 = \dot{y}_1 , \quad y_3 = \ddot{y}_1 = \dot{y}_2 , \ldots , \quad y_{n} = \dot{y}_{n-1} , \]

we convert the given differential equation into the system of first order equations:

\[ \frac{\text d}{{\text d}t} \begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{bmatrix} = {\bf C}_{p} \begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ \vdots \\ f(t) \end{bmatrix} . \]

where matrix

\begin{equation} \label{EqConv.6} {\bf C}_{p} = \begin{bmatrix} 0&1&0& \cdots &0&0 \\ 0&0&1& \cdots &0&0 \\ 0&0&0 \cdots &0&0 \\ \vdots&\vdots&\vdots& \ddots&\vdots&\vdots \\ 0&0&0& \cdots &0&1 \\ - a_0 & -a_1 & -a_2 & \cdots & -a_{n-2}& -a_{n-1} \end{bmatrix} \end{equation}

is called companion matrix for polynomial p(λ).

Theorem: For any polynomial p(λ), its companion matrix has the characteristic polynomial

\begin{equation} \label{EqConv.7} \chi (\lambda ) = \det \left( \lambda {\bf I} - {\bf C}_{p} \right) = p(\lambda ) \end{equation}

equal to the given polynomial.    ▣

Example 3: Inverse of a companion matrix

Example 3: Consider a polynomial of the fourth degree \( p(\lambda ) = \lambda^4 + a\,\lambda^3 + b\, \lambda^2 + c\,\lambda + d , \) where d ≠ 0, and its companion matrix

\[ {\bf C}_4 = \begin{bmatrix} 0&1&0&0 \\ 0&0&1&0 \\ 0&0&0&1 \\ -d& -c& -b& -a \end{bmatrix} \qquad \Longrightarrow \qquad {\bf C}_4^{-1} = \begin{bmatrix} -c/d & -b/d & -a/d & -1/d \\ 1&0&0&0 \\ 0&1&0&0 \\ 0&0&1&0 \end{bmatrix} . \]

Click to see matlab code

function Converting_3
%% Mathematica code:
% Inverse[{{0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}, {-d, -c, -b, -a}}]
syms a b c d
X = [0 1 0 0; 0 0 1 0; 0 0 0 1; -d -c -b -a]
Y = inv(X)
end

Copy

Thus, the inverse of this companion matrix corresponds to the polynomial:

\[ \lambda^4 + \frac{c}{d}\, \lambda^3 + \frac{b}{d}\, \lambda^2 + \frac{a}{d}\, \lambda + \frac{1}{d} . \qquad \blacksquare \]

The same approach can be used to convert a nonlinear differential equation to an equivalent system of first order equations. It should be noted that nonlinear systems of differential equations can be solved only numerically (except for some rare cases). When a numerical solver is applied, it is common to introduce an auxiliary variable---the residual---in order to control accuracy. This approach is demonstrated in the following series of examples.

Example 4: Pendulum equation

Example 4: Consider the initial value problem for the ideal pendulum equation

\[ \ddot{\theta} + \omega^2 \sin \theta =0 , \qquad \theta (0) = a , \quad \dot{\theta} (0) = v , \tag{4.1} \]

where initial displacement 𝑎 and velocity v are specified, and ω² = (g/ℓ), with g being acceleration due to gravity and ℓ being length of the rod containing as point mass. This equation is called the ideal pendulum equation because friction is not taken into account. When resistance is presented and assumed to be proportional to the velocity of the bob, the equation becomes

\[ \ddot{\theta} + k\,\dot{\theta} + \omega^2 \sin \theta =0 , \qquad \theta (0) = a , \quad \dot{\theta} (0) = v , \tag{4.2} \]

where k is a positive coefficient.

These problems (for ideal pendulum and pendulum with resistance) can be converted to equivalent systems of differential equations of the first order:

\[ \begin{cases} \dot{x} = y , \\ \dot{y} = - \omega^2 \sin x , \end{cases} \qquad x(0) =a, \quad y(0) = v ; \tag{4.3} \]

and

\[ \begin{cases} \dot{x} = y , \\ \dot{y} = - \omega^2 \sin x - k\,y , \end{cases} \qquad x(0) =a, \quad y(0) = v , \tag{4.4} \]

respectively. We plot its phase portraits using matlab:

 

The ideal pendulum equation can be converted to a system of three equations upon introducing the residual variable:

\[ \begin{cases} \dot{x} = y , \\ \dot{y} = z , \\ \dot{z} = - \omega^2 y\,\cos x , \end{cases} \qquad x(0) =a, \quad y(0) = v , \quad z(0) = - \omega^2 \sin a . \]

So z is actually the second derivative of the original unknown function x(t). Therefore, the last equation is obtained by differentiation of the original pendulum equation with respect to time variable t. WE control approximation of the pendulum differential equation with the following Mathematica code:

		Since the function z(t) is supposed to be the second derivative, we can control the accuracy of matlab approximation of the pendulum equation by evaluating \( \ddot{\theta} + \sin \theta . \) Click to see matlab code % Ideal pendulum equation converted to a system of three equations % Define pendulum function to be solved pf = @(t,y) [y(2);y(3); -y(2)*cos(y(1))]; % Solve pendulum function using ODE45 solver and initial conditions [t,y] = ode45(pf,[0 20] , [pi/3,-1/10,-sin(pi/3)]); % Create figure and plot solution figure hold on plot(t, y(:,3)+sin(y(:,1)),'DisplayName','l=1','LineWidth',1.5) % Add x and y grid lines ax = gca; ay = gca; ax.YGrid = 'on'; ay.XGrid = 'on'; Copy
Plot of \( \ddot{\theta} + \sin \theta . \)		matlab code

Example 5: Pendulum equation again

Example 5: In the previous example, we converted the ideal pendulum equation \( \ddot{\theta} + \sin \theta =0 \) into the system of first order equations upon substituting y₁ = θ and let y₂ to be its derivative. Now we may define a new variable y₃ = sinθ. According to the chain rule, its derivative is

\[ \frac{{\text d}y_3}{{\text d}t} = \left(\frac{{\text d}y_1}{{\text d}t} \right) \cos y_1 = y_2 \cos \theta . \]

Then we introduce a fourth dependent variable y₄ = cosθ and obtain for it

\[ \frac{{\text d}y_4}{{\text d}t} = - \left(\frac{{\text d}y_1}{{\text d}t} \right) \sin y_1 = - y_2 \,y_3 . \]

Finally, we combine all equations together:

\begin{align*} \dot{y}_1 &= y_2 , \\ \dot{y}_2 &= -y_3 , \\ \dot{y}_3 &= y_2\, y_4 , \tag{5.1} \\ \dot{y}_4 &= -y_2\, y_3 . \end{align*}

The initial conditions for y₃ and y₄ will be related to those imposed on y₁ and y₂:

\[ y_1 (0) = a, \quad y_2 (0) = v, \quad y_3 (0) = \sin a, \quad y_4 (0) = \cos a. \]

This polynomial system of equations is suitable for Picard's iterations (see section in this tutorial):

\[ {\bf \phi}_{m+1} (t) = \begin{bmatrix} y_1 (t) \\ y_2 (t) \\ y_3 (t) \\ y_4 (t) \end{bmatrix}_{(m+1)} = \begin{bmatrix} a \\ v \\ \sin a \\ \cos a \end{bmatrix} + \int_0^t {\text d}s \begin{bmatrix} y_2 (s) \\ -y_3 (s) \\ y_2 (s)\,y_4 (s) \\ -y_2 (s)\,y_3 (s) \end{bmatrix}_{(m)} , \qquad m= 0,1,2,\ldots . \tag{5.2} \]

Performing a couple of iterations, we obtain

\[ \theta_{(3)} (t) = a + tv - \frac{t^2}{2}\, \sin a - \frac{t^3}{6}\,v\,\cos a . \]

phi1[t_] = {a, v, Sin[a], Cos[a]} + Integrate[{v, -Sin[a], v*Cos[a], -v*Sin[a]}, {s, 0, t}];
phi2[t_] = {a, v, Sin[a], Cos[a]} + Integrate[{phi1[s][[2]], -phi1[s][[3]], phi1[s][[2]]*phi1[s][[4]], -phi1[s][[2]]*phi1[s][[3]]}, {s, 0, t}];
phi3[t_] = {a, v, Sin[a], Cos[a]} + Integrate[{phi2[s][[2]], -phi2[s][[3]], phi2[s][[2]]*phi2[s][[4]], -phi2[s][[2]]*phi2[s][[3]]}, {s, 0, t}];
phi4[t_] = {a, v, Sin[a], Cos[a]} + Integrate[{phi3[s][[2]], -phi3[s][[3]], phi3[s][[2]]*phi3[s][[4]], -phi3[s][[2]]*phi3[s][[3]]}, {s, 0, t}];
phi3[t][[1]]

a + t v - 1/6 t^3 v Cos[a] - 1/2 t^2 Sin[a]

pphi4[t_] = phi4[t][[1]]

Now we plot the third approximation along with the true solution:

plot = Plot[{Evaluate[y[t] /. s], pphi3[t]}, {t, 0, 8}, PlotRange -> All, PlotStyle -> Thick, Epilog -> {Text[Style["approximation", 25, Orange], {5.2, 6}]}]

   ◾

Example 6: Forced Duffing's equation

Example 6: Consider the initial value problem for Duffing's equation

\[ \ddot{x} + k\,\dot{x} + x - \varepsilon\,x^3 = f(t) , \qquad x(0) =a, \quad \dot{x} = v , \tag{6.1} \]

where f(t) is a periodic forcing function, which can be chosen as \( f(t) = A\,\cos \omega t \) or \( f(t) = A\,\sin \omega t \) or their combination \( f(t) = A\,\sin \left( \omega t + \phi \right) . \) Here A is the amplitude, k is the dumping coefficient, and ε is a real number.

We convert it to the system of two differential equations

\[ \begin{cases} \dot{x} = y , \\ \dot{y} = - k\,y -x + \varepsilon\,x^3 + f(t) , \end{cases} \qquad x(0) =a, \quad y(0) = v . \tag{6.2} \]

Its direction field for undriven case is plotted with the following commands:

StreamPlot[{y, -y - x - x^3}, {x, -6, 6}, {y, -2, 2}, StreamScale -> Medium]

The residual system of equations becomes

\[ \begin{cases} \dot{x} = y , \\ \dot{y} = z , \\ \dot{z} = - k\,z - y - 3\varepsilon \, x^2 y - \dot{f} , \end{cases} \qquad x(0) =a, \quad y(0) = v , \quad z(0) = - k\,v - a - \varepsilon \, a^3 + f(0). \]

   ◾

Example 7: Single linear differential equation of third order

Example 7: Consider a third order constant coefficient differential equation

\[ y''' +2\, y'' - y' -2 \,y = \sin (2t) . \tag{7.1} \]

Upon introducing three new dependent variables

\[ y_1 = y(t), \quad y_2 = \dot{y} = \dot{y}_1 , \quad y_3 = \ddot{y} = \dot{y}_2 , \tag{7.2} \]

we transfer the given single third order differential equation to the system of three equations of the first order

\begin{align*} \dot{y}_1 &= y_2 , \\ \dot{y}_2 &= y_3 , \tag{7.3} \\ \dot{y}_3 &= - 2\,y_3 + y_2 + 2\, y_1 + \sin (2t) . \end{align*}

We rewrite this system of equations in compact vector form:

\[ \frac{\text d}{{\text d}t} \begin{bmatrix} y_1 (t) \\ y_2 (t) \\ y_3 (t) \end{bmatrix} = \begin{bmatrix} 0&1&0 \\ 0&0&1 \\ -2& 1& 2 \end{bmatrix} \begin{bmatrix} y_1 (t) \\ y_2 (t) \\ y_3 (t) \end{bmatrix} + \begin{bmatrix} 0 \\ 0 \\ \sin (2t) \end{bmatrix} \tag{7.4} \]

The characteristic polynomial of the companion matrix

\[ \det \left( \lambda \begin{bmatrix} 1&0&0 \\ 0&1&0 \\ 0&0&1 \end{bmatrix} - \begin{bmatrix} 0&1&0 \\ 0&0&1 \\ 2&1&-2 \end{bmatrix} \right) = \lambda^3 + 2\,\lambda^2 - \lambda - 2\,\lambda = \left( \lambda^2 -1 \right) \left( \lambda +2 \right) \]

is in agreement with the characteristic polynomial for the given single differential equation.    ◾

Example 8: Blasius equation

Example 8: The Blasius equation for laminar flow over a flat plate can be reduced to the following differential equation of the third order:

\[ f''' + \frac{1}{2}\, f \, f'' =0 . \tag{8.1} \]

Upon introducing new functions:

\begin{eqnarray*} u_1 &=& f , \\ u_2 &=& f' = u'_1 , \\ u_3 &=& f'' = u'_2 , \end{eqnarray*}

we reduce the Blasius equation to a system of first order equations:

\[ \frac{{\text d}}{{\text d}y} \begin{bmatrix} u_1 \\ u_2 \\ u_3 \end{bmatrix} = \begin{bmatrix} u_2 \\ u_3 \\ - \frac{1}{2}\, u_1 u_2 \end{bmatrix} , \tag{8.2} \]

which is a nonlinear vector differential equation for three unknowns.

 

Classical Picard's iteration

---

In order to apply the classical Picard's iteration, we have to convert this to a system of first-order differential equations. Let

\[ f_1 = f, \qquad f_2 = f' , \qquad\mbox{and} \qquad f_3 = f'' . \]

This leads to:

\[ \frac{\text d}{{\text d}y} \begin{bmatrix} f_1 \\ f_2 \\ f_3 \end{bmatrix} = \begin{bmatrix} f_2 \\ f_3 \\ - \frac{1}{2}\, f_1 f_3 \end{bmatrix}, \qquad f_1 (0) =0, \quad f_2 (0) =0, \quad f_2 (\infty ) =1. \]

It is not possible to specify a boundary condition at ∞ numerically, so we will have to use a large number, and verify it is "large enough". From the solution, we evaluate the derivatives at y = 0, and we have f''(0) = f₃(0).

We have to provide initial guesses for f₁, f₂ and f₃. This is the hardest part about this problem. We know that f₁ starts at zero, and is flat there (f'(0)=0), but at large y, it has a constant slope of one. We will guess a simple line of slope = 1 for f₁. That is correct at large y, and is zero at y = 0. If the slope of the function is constant at large y, then the values of higher derivatives must tend to zero. We choose an exponential decay as a guess.

Finally, we let a solver iteratively find a solution for us, and find the answer we want.

We get the initial second derivative to be f''(0) = 0.3324911095517483. For the sake of simplicity, we approximate the Blasius constant as 1/3.

Now we use Picard iteration procedure starting with the initial conditions. Since the derivative operator is unbounded, we apply its inverse to reduce the given initial value problem into fixed point form:

\[ \frac{\text d}{{\text d}y} \begin{bmatrix} f_1 \\ f_2 \\ f_3 \end{bmatrix} = \begin{bmatrix} f_2 \\ f_3 \\ - \frac{1}{2}\, f_1 f_3 \end{bmatrix}, \qquad \begin{bmatrix} f_1 (0) \\ f_2 (0) \\ f_3 (0) \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ 1/3 \end{bmatrix} \qquad \Longleftrightarrow \qquad \begin{bmatrix} f_1 (y) \\ f_2 (y) \\ f_3 (y)\end{bmatrix} = \begin{bmatrix} f_1 (0) \\ f_2 (0) \\ f_3 (0) \end{bmatrix} + \int_0^y \begin{bmatrix} f_2 (s) \\ f_3 (s) \\ - \frac{1}{2}\, f_1 (s)\,f_3 (s) \end{bmatrix} {\text d}s . \]

From this fixed point equation, we derive the iteration procedure

\[ \begin{bmatrix} f_{1,n+1} (y) \\ f_{2,n+1} (y) \\ f_{3,n+1} (y) \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ 1/3 \end{bmatrix} + \int_0^y \begin{bmatrix} f_{2,n} (s) \\ f_{3, n} (s) \\ - \frac{1}{2}\, f_{1,n} (s)\,f_{3,n} (s) \end{bmatrix} {\text d}s , \qquad n=0,1,2,\ldots . \]

The initial approximation is just the initial conditions:

\[ \phi_{1,0} =0, \quad \phi_{2,0} =0, \quad \phi_{3,0} = 1/3 \quad\mbox{however, it should be the Blasius constant}. \]

The first iteration becomes

\[ \begin{bmatrix} \phi_{1,1} (y) \\ \phi_{2,1} (y) \\ \phi_{3,1} (y) \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ 1/3 \end{bmatrix} + \int_0^y \begin{bmatrix} 0 \\ 1/3 \\ 0 \end{bmatrix} {\text d} y = \begin{bmatrix} 0 \\ y/3 \\ 1/3 \end{bmatrix} . \]

The second one is

\[ \begin{bmatrix} \phi_{1,2} (y) \\ \phi_{2,2} (y) \\ \phi_{3,2} (y) \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ 1/3 \end{bmatrix} + \int_0^y \begin{bmatrix} y/3 \\ 1/3 \\ 0 \end{bmatrix} {\text d} y = \begin{bmatrix} y^2 /6 \\ y/3 \\ 1/3 \end{bmatrix} = \frac{1}{3} \begin{bmatrix} y^2 /2 \\ y \\ 1 \end{bmatrix} . \]

The third one is

\[ \begin{bmatrix} \phi_{1,3} (y) \\ \phi_{2,3} (y) \\ \phi_{3,3} (y) \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ 1/3 \end{bmatrix} + \int_0^y \begin{bmatrix} y/3 \\ 1/3 \\ - y^2/36 \end{bmatrix} {\text d} y = \begin{bmatrix} y^2 /6 \\ y/3 \\ 1/3 - y^3/108 \end{bmatrix} = \frac{1}{3} \begin{bmatrix} y^2 /2 \\ y \\ 1 - y^3/36 \end{bmatrix} . \]

For next one, we have

\[ \begin{bmatrix} \phi_{1,4} (y) \\ \phi_{2,4} (y) \\ \phi_{3,4} (y) \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ 1/3 \end{bmatrix} + \int_0^y \begin{bmatrix} y/3 \\ 1/3 - y^3/108 \\ - \frac{1}{36} \left( y^2 - \frac{y^5}{36} \right) \end{bmatrix} {\text d} y = \begin{bmatrix} y^2 /6 \\ \frac{y}{3} - \frac{y^4}{4 \cdot 108} \\ \frac{y^6}{36^3 \cdot 6} - \frac{y^3}{6\cdot 26} \end{bmatrix} = \frac{1}{3} \begin{bmatrix} y^2 /2 \\ y - \frac{y^4}{4 \cdot 36} \\ \frac{y^6}{2\cdot 36^2} - \frac{y^3}{36} \end{bmatrix} . \]

So our fifth approximation becomes

f1[y_] = y^2/6; f2[y_] = 1/3*(y - y^4/4/36); f3[y_] = 1/3*(y^6/2/36^2 - y^3 /36); f5[y_] = Simplify[Integrate[1/9*(s - s^4 /4/36), {s, 0, y}]]

-((y^2 (-360 + y^3))/6480)

\[ f_5 (y) = \frac{y^2}{6480} \left( 360 - y^3 \right) . \]

   ◾

Example 9: Genesio equation

Example 9: The generalized Genesio differential equation is of the form

\[ \dddot{x} + a\,\ddot{x} + b\,\dot{x} + c\,x + h(x) = 0 , \]

where \( \displaystyle h(x) = \sum_{k=1}^n (-1)^k x^{2k} \) and 𝑎, b, and c are arbitrary real parameters. We consider a simple version of the Genesio equation when chaotic jerk is observed:

\[ \dddot{x} + \dot{x} - x^2 + x^4 = 0 , \]

Upon introducing new functions:

\begin{eqnarray*} x_1 &=& x , \\ x_2 &=& x' = x'_1 , \\ x_3 &=& x'' = x'_2 , \end{eqnarray*}

we reduce the Genesio equation to a system of first order equations:

\[ \frac{{\text d}}{{\text d}t} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} = \begin{bmatrix} x_2 \\ x_3 \\ - x_2 + x_1^2 - x_1^4 , \end{bmatrix} , \]

which is a nonlinear vector differential equation for three unknowns.

We make a numerical experiment and plot two solutions of the Genesio equation for small values of parameters 𝑎 = 0.1 and c = 0.1, and then compare with the case when they are zero. Click to see matlab code f = @(t,y) [y(2);y(3); y(1)^2 - y(1)^4-y(2)]; init_cond = [1, -.1, 0]; t_interval = [0,3]; [t,y] = ode45(f,t_interval , init_cond); figure plot(t, y(:,1),'linewidth',2) title('Genesio equation with a = c = 0') Copy

Then we repeat calculations with small values of parameters 𝑎 = c =0.1:

We make a numerical experiment and plot two solutions of the Genesio equation for small values of parameters 𝑎 = 0.1 and c = 0.1, and then compare with the case when they are zero. Click to see matlab code f = @(t,y) [y(2);y(3); y(1)^2 - y(1)^4-y(2)]; init_cond = [1, -.1, 0]; t_interval = [0,3]; [t,yb] = ode45(f,t_interval , init_cond); figure plot(t, yb(:,1)) title('Genesio equation with a = c = 0') %Repeat with small values of a and c and graph with 1.3b figure plot(t, yb(:,1),'b','linewidth',2) hold on g = @(t,y) [y(2);y(3); y(1)^2 - y(1)^4-y(2)-.1*y(1) - .1*y(3)]; init_cond = [1, -.1, 0]; t_interval = [0,3]; [t,y] = ode45(g,t_interval , init_cond); plot(t, y(:,1),'r',linewidth=2) title('Genesio equation with a = c = 0.1') Copy

   ◾

References

1. Brand, L., The companion matrix and its properties, The American Mathematical Monthly, 1964, Vol. 71, No. 6, pp. 629--634.
2. Eastham, M.S.P., "The spectral theory of periodic differential equations" , Scottish Acad. Press (1973)
3. Hill, G. W. "On the Part of the Motion of Lunar Perigee Which is a Function of the Mean Motions of the Sun and Moon." Acta Math. 8, 1-36, 1886.
4. Magnus,W., Winkler, S., "Hill's equation" , Dover, reprint (1979)
5. V.A. Yakubovich, V.M. Starzhinskii, "Linear differential equations with periodic coefficients" , Wiley (1975) (Translated from Russian)