Chapter Name

9.1-1. Survey. Differential equations are used to express relations between changes in physical quantities and are thus of great importance in many applications.Sections 9.1-2 to 9.3-10 present a straightforward classical introduction to ordinary differential equations, including some complex-variable theory. Sections 9.4-1 to 9.4-8 introduce the linear differential equations with constant coefficients used in the analysis of vibrations, electric circuits, and control systems, with emphasis on solutions by Laplace-transform methods. Sections 9.5-1 to 9.5-6 deal with nonlinear second-order equations. Sections 9.6-1 and 9.6-2 introduce Pfaffian differential equations, although these are not ordinary differential equations.

Some naturally related material is treated in other chapters of this handbook, particularly in Chap. 8 and Secs. 13.6-1 to 13.6-7. Boundary value problems, eigenvalue problems, and orthogonal-function expansions of solutions are discussed in Chap. 15, and a number of differential equations defining special functions are treated in Chap. 21.

The notation used in the various subdivisions of this chapter has been chosen so as to simplify reference to standard textbooks in different special fields. Thus the usually real variables in Secs. 9.2-1 to 9.2-5 are denoted by x, y = y(x); the frequently complex variables encountered in the general theory of linear ordinary differential equations (Secs. 9.3-1 to 9.3-10) are denoted by z, w = w(z). The variables in Sees. 9.4-1 to 9.5-6 usually represent physical time and various mechanical or electrical variables and are thus introduced as t, yk = yk(t).

9.1-2. Ordinary Differential Equations. An ordinary differential equation of order r is an equation

to be satisfied by the function y = y(x) together with its derivatives y({x), y′′{x), . . . , y^(r)(x) with respect to a single independent variable x. To solve (integrate) a given differential equation (1) means to find functions (solutions, integrals) y(x) which satisfy Eq. (1) for all values of x in a specified bounded or unbounded interval (a, b). Note that solutions can be checked by resubstitution.

The complete primitive (complete integral, general solution) of an ordinary differential equation of order r has the form

where C₁, C₂, . . . , C_r are r arbitrary constants (constants of integration, see also Sec. 4.6-4). Each particular choice of these r constants yields a particular integral (2) of the given differential equation. Typical problems require one to find the particular integral (2) subject to r initial conditions

which determine the r constants C₁, C₂, . . . , C_r. Alternatively, one may be given r boundary conditions on y (x) and its derivatives for x = a and x = b (see also Sec. 9.3-4).*

Many ordinary differential equations admit additional solutions known as singular integrals which are not included in the complete primitive (2) (see also Sec. 9.2-26).

A differential equation is homogeneous if and only if ay(x) is a solution for all a whenever y(x) is a solution (see also Secs. 9.1-5 and 9.3-4).

* Strictly speaking, initial and boundary conditions refer to unilateral derivatives (Sec. 4.5-1).

Given an r-parameter family of suitably differentiable functions (2), one can eliminate C₁, C₂, . . . , C_r from the r + 1 equations y⁽ⁱ⁾ = y⁽ⁱ⁾(x, C₁, C₂, . . . , C_r) (j = 0, 1, 2, . . . , r) to obtain an r^th-order differential equation describing the family.

NOTE: An ordinary differential equation is a special instance of a functional equation imposing conditions on the functional dependence y = y(x) for a set of values of x. OTHER EXAMPLES OF FUNCTIONAL EQUATIONS: y(x₁x₂) = y(x₁) + y(x₂) [logarithmic property, satisfied by y(x) = A log, x], partial differential equations (Sec. 10.1-1), integral equations (Sec. 15.3-2), and difference equations (Sec. 20.4-3).

9.1-3. Systems of Differential Equations(see also Secs. 13.6-1 to 13.6-7). A system of ordinary differential equations

involves a set of unknown functions y₁= y₁(x), y₂ = y₂(x), . . . and their derivatives with respect to a single independent variable x. The order r_i of each differential equation (4) is that of the highest derivative occurring. In general, one will require n differential equations (4) to determine n unknown functions y_k(x); and the general solution y₁ = y₁(x), y₂ = y₂(x) will involve a number of arbitrary constants equal to r = r₁+ r₂+ . . . +r_n.

The solution of a system (4) can be reduced to that of a single ordinary differential equation of order r through elimination of n — 1 variables yk and their derivatives. More importantly, one can reduce every system (4) to an equivalent system of r first-order equations by introducing higherorder derivatives as new variables.

9.1-4. Existence and Desirable Properties of Solutions. A properly posed differential-equation problem requires an existence proof indicating the construction of a solution subject to the given type of initial or boundary conditions. The existence of physical phenomena described by a given differential equation may suggest but does not prove the existence of a solution; an existence proof checks the self-consistency of the mathematical model (see also Secs. 4.2-1b and 12.1-1; see Secs. 9.2-1 and 9.3-5 for examples of existence theorems).

It is desirable to design mathematical models involving differential equations so that the solutions are continuous functions of numerical coefficients, initial conditions, etc., so as to avoid excessive errors in solutions due to small errors in numerical data (see also Sec. 9.2-la).

9.1-5. General Hints. (a) Substitution of a Taylor series (Sec. 4.10-4) or other series expansion for y(x) in a given differential equation may yield equations for the unknown coefficients (see also Secs. 9.2-5b and 9.3-5). Many differential equations can be simplified through transformation of variables (Secs. 9.1-5b, 9.1-3, and 9.3-8c). Every differential equation or system of differential equations can be reduced to a system of first-order equations; so that the methods of Sec. 9.2-5 apply.

(b)The following special types of differential equations reduce easily to equations of lower order (see also Secs. 9.2-3 and 9.5-6):

If a given differential equation F(x, y, y′, y′′, . . . , y^r) = 0 is homogeneous in the arguments y, y′, y′′, . . . , y^(r) (Sec. 4.5-5; this does not necessarily imply that the differential equation is homogeneous in the sense of Sec. 9.1-2), introduce ȳ = y’y.

9.2. FIRST-ORDER EQUATIONS

9.2-1. Existence and Uniqueness of Solutions. (a) A given first-order differential equation expressible in the form

has a solution y = y(x) through every “point” (x = x₀, y = y₀) with a neighborhood throughout which f(x, y) is continuous. More specifically, let D be a region of “points” (x, y), (x, n) where f(x, y) is single-valued, bounded, and continuous and

for some real M independent of y and n. Then the given differential equation (1) has a unique solution y = y(x) through every point (x = x₀, y = y_o) of D, and y(x) is a continuous function of the given value yo = y(x₀). Each solution extends to the boundary of D.

The Lipschitz condition (2) is satisfied, in particular, whenever f(x, y) has a bounded and continuous derivative ∂f/∂y in D.

(b) (See also Sec. 9.1-3). An analogous existence theorem applies to systems of first-order differential equations

if the Lipschitz condition (2) is replaced by

9.2-2. Geometrical Interpretation. Singular Integrals (see also Secs. 17.1-1 to 17.1-7). (a) If x, y are regarded as rectangular cartesian coordinates, a first-order differential equation

describes a “field” of line elements (x, y, p) or elements of straight lines through (x, y) with slope p = dy/dx = f(x, y). Each line element is tangent to a curve of the one-parameter family of solutions

where λ is a constant of integration.

A plot of the field of tangent directions permits at least rough graphical determination of solutions; the general character of the family of solutions may be further discussed in the manner of Sec. 9.5-2. It may be helpful to know that the curves F(x, y, p₁) = 0 or f(x, y) = P₁ are isoclines where the solution curves have a specified fixed slope p₁.

The curves are loci of points of inflection (see also Sec. 9.5-2).

(b)Singular Integrals(see also Sec. 9.1-2). Let F(x, y, p) be twice continuously differentiate with respect to x and y, and let ∂F/∂y ≠ 0. Elimination of p from

yields a curve or set of curves called the p discriminant of the given differential equation (locus of singular line elements). A curve defined by Eq. (7) is a singular integral of the given differential equation if on this curve, unless both ∂F/∂x and ∂F/∂y vanish at a point of the curve. Geometrically, such singular integrals are frequently envelopes of the family of solution curves (6), and may thus be obtained from the complete primitive (6) in the manner of Sec. 17.1-7.

9.2-3. Transformation of Variables(see Sec. 9.2-4 for examples). (a) A suitable continuously differentiable transformation

will transform the given differential equation (1) or (5) into a new differential equation relating and . The new equation may be simpler, or a solution may be known. Once is found, y = y(x) is given implicitly, or by inverse transformation.

(b) Contact Transformations (see also Secs. 10.2-5, 10.2-7, and 11.6-8). A set of twice continuously differentiable transformation equations

With the special property

defines a contact transformation associating line elements (Sec. 9.2-2a) (x, y, p) and so that line elements forming regular arcs are mapped onto regular arcs, and contact of regular arcs is preserved. It is then legitimate to write , and to use suitable contact transformations (9) to simplify differential equation (1) or (5). Once a solution of the transformed equation is known, y = y(x) is given implicitly or by inverse transformation.

In particular, g(x, y, p) =1 yields the easily reversible contact transformation

which transforms a given differential equation (5) into

Equation (11) may be a simpler differential equation or, indeed, an ordinary equation relating x̄ and ȳ.

9.2-4. Solution of Special Types of First-order Equations. (a) The following special types of first-order equations are relatively easy to solve.

1. The variables are separable: y' = f₁(x)/f₂(y). Obtain the solution from ∫f₂(y) dy = ∫f₁(x) dx + C.

2. “Homogeneous” first-order equations:* y′ = f(y/x). Introduce = y/x to reduce to type 1.

3. Exact differential equations can be written in the form

where the expression on the left is an exact differential , Sec. 5.7-1 Obtain the solution from

* Note that the expression “homogeneous” differential equation is here not used in the sense defined in Sec. 9.1-2.

If the expression on the left of Eq. (12) is not an exact differential , one may be able to find an integrating factor µ = µ(X, y) such that multiplication of Eq. (12) by µ(x, y) yields an exact differential equation. The integrating factor n(x, y) satisfies the partial differential equation

4. The linear first-order equation y′ + a(x)y = f(x) (see also Sees. 9.3-1 and 9.3-3) admits the integrating factor

The complete primitive is then

Many first-order equations can be reduced to one of the above types by transformation of variables (Sec. 9.2-3). In particular

y' = f(ax + βy) reduces to type 1 if one introduces = ax + βy.

reduces to type 2 by a coordinate translation if α₁β₂ - α₂β₁ ≠ 0; otherwise introduce to separate the variables.

y′ = f₁(x)y + f₂(x)yⁿ (BERNOULLI’s DIFFERENTIAL EQUATION) reduces to a linear equation if one introduces .

(b) Given a first-order equation of the form

it may be advantageous to differentiate both sides with respect to x. The resulting differential equation

might be easy to solve for y′ = y′(x) or y′ = y′(y, respectively; substitution of this result into the given Eq. (16) yields the desired relation of x and y. If the solution of Eq. (17) takes the form u(x, y′) = 0 or u(y, y′) = 0, the desired relation of x and y is given in terms of a parameter p = y′.

EXAMPLES: Clairaut’s differential equation y = y′x + f(y′) yields the complete primitive y = Cx + f(C) and the singular integral (in parametric representation) x = —ƒ′(p), y = —pƒ′(p) + f(P). Lagrange’s differential equation y = xf₁(p) + f₂(p) is solved in the same manner.

(c)Riccati Equations. The differential equation

is sometimes simplified by the transformation y = 1/ȳ alternatively,

leads to a homogeneous second-order equation for :

If a particular integral y₁(x) of Eq. (18) is known, the transformation

yields a linear differential equation. If one knows two particular integrals y₁, y₂ or three particular integrals y_h yz, y_z, one has, respectively,

For any four particular integrals y₁, y₂, y₃, y₄, the double ratio (y₁ — y₂) (y₃ — y₄)/(y₁ — y₃) (y₂ — y₄) is constant.

The special Riccati equation

can be reduced to type 1 if m = 4k/(l — 2k) (k = 0, ±1, ±2, . . .). For k > 0, the transformation reduces Eq. (21) to a similar equation

The procedure is repeated until (after k steps) the right side of the differential equation is constant.

Similarly, for k < 0, the transformation x = x̄^—1/(m+1), yields a differential equation of the form (22a) with

9.2-5. General Methods of Solution. (a) Picard’s Method of Successive Approximations. To solve the differential equation y′ = f(x, y) for a given initial value y(x_o) = y₀, start with a trial solution y^[0](x) and compute successive approximations

to the desired solution y(x). The process converges subject to the conditions of Sec. 9.2-1. Picard’s method is useful mainly if the integrals in Eq. (23) can be evaluated in closed form, although numerical integration can, in principle, be used.

A completely analogous procedure applies to systems (3) of first-order differential equations.

(b)Taylor-series Expansion (see also Sec. 4.10-4). If the given function f(x,y) is suitably differentiable, obtain the coefficients y^(m)(x_o)/mǃ of the Taylor series

by successive differentiations of the given differential equation:

with x = x_o, y = y(x₀) = y₀.

An analogous procedure applies to systems of first-order equations.

9.3.LINEAR DIFFERENTIAL EQUATIONS

9.3-1.Linear Differential Equations. Superposition Theorems (see also Secs. 10.4-2, 13.6-2, 13.6-3, 14.3-1, and 15.4-2). A linear ordinary differential equation of order r relating the real or complex variables z and w = w(z) has the form

where the a_k(z) and f(z) are real or complex functions of z. The general solution (Sec. 9.1-2) of a linear differential equation (1) can be expressed as the sum of any particular integral and the general solution of the homogeneous linear differential equation (Sec. 9.1-2)

For any given nonhomogeneous or “complete” linear differential equation (1), the homogeneous equation (2) is known as the complementary equation or reduced equation, and its general solution as the complementary function.

Let w₁(z) and w₂(z) be particular integrals of the linear differential equation (1) for the respective “forcing functions” f(z) ≡ f₁(z) and f(z) ≡ f₂(z). Then αw₁(z) + β w₂(z) is a particular integral for the forcing function f(z) ≡ αf₁(z) + βf₂(z) (Superposition Principle). In particular, every linear combination of solutions of a homogeneous linear differential equation (2) is also a solution.

The superposition theorems often represent some physical superposition principle. Mathematically, they permit one to construct solutions of Eq. (1) or (2) subject to given initial or boundary conditions by linear superposition.

Analogous theorems apply to systems of linear differential equations (see also Sec. 9.4-2).

9.3-2.Linear Independence and Fundamental Systems of Solutions (see also Secs. 1.9-3, 14.2-3, and 15.2-la). (a) Let w₁(z), w₂(z), . . . , w_r(z) be r — 1 times continuously differentiate solutions of a homogeneous linear differential equation (2) with continuous coefficients in a domain D of values of z. The r solutions w_k(z) are linearly independent in D if and only if in D implies λ₁ = λ₂ = . . . = λ_r = 0 (Sec. 1.9-3). This is true if and only if the Wronskian determinant (Wronskian)

differs from zero throughout D. W = 0 for any z in D implies w = 0 for all z in D.*

(b) A homogeneous linear differential equation (2) of order r has at most r linearly independent solutions, r linearly independent solutions w₁(z), w₂(z), . . . , w_k(z) constitute a fundamental system of solutions

whose linear combinations include all particular integrals of Eq. (2).

(c) Use of Known Solutions to Reduce the Order. If m < r linearly independent solutions w₁(z), w ₂(z), . . . , w_m(z) of the homogeneous equation (2) are known, then the transformation w̄ = W[w₁, w₂, . . . , w_m, w]φ( x) reduces Eq. (2) to a homogeneous linear differential equation of order r — m for any conveniently chosen φ(x).

9.3-3. Solution by Variation of Constants. Green’s Functions. (a) Given r linearly independent solutions w₁(z), w₂(z), . . . , w_r(z) of the homogeneous linear differential equation (2), the general solution of the complete nonhomogeneous equation (1) is

* Note that the theorem in this simple form does not apply to every set of r — 1 times continuously differentiable functions W_k(z); they must be solutions of a suitable differential equation (2).

After solving the r simultaneous equations (5) for the r unknown derivatives C′_k(z), one obtains each C_k(z) = ∫ C′_k(z) dz + K_k by a simple integration. In principle, this procedure reduces the solution of any linear ordinary differential equation to the solution of a homogeneous linear differential equation.

(b) Assuming real variables z ≡ x and w ≡ w (x) for simplicity, particular integrals of the complete differential equation (1) can often be written as

where G(x, ξ) is known as the Green’s function (sometimes called the weighting function, Sec. 9.4-3) yielding the specific particular integral in question. The complete integral of Eq. (1) is then

where the w_k(z) are r linearly independent solutions of Eq. (2), and the A_k are r constants of integration to be determined by suitable initial or boundary conditions.

Any given set of r linearly independent solutions w_k(x) of the complementary equation (2) permits one to construct a particular integral (6) with

where the C′ _k(x) are obtained from Eq. (5), and U(x) is the unit-step function defined in Sec. 21.9-1.

For linear differential equations of order r = 2, the complete integral is given by Eq. (4) or (7) with

(c) While the general solution (7) obtained with the aid of the particular Green’s function (8) is only another way of writing Eq. (4), it is often possible to construct a Green’s function G (x, ξ) such that the particular integral (6) satisfies the specific initial or boundary conditions of a given problem. Assuming boundary conditions linear and homogeneous in w(x) and its derivatives, the required Green’s function G(x, ξ ) must satisfy the given boundary conditions and

for x in (a, b), with ∂^r—2G/∂x^r—2 continuous in (a, b), and

The existence and properties of such Green’s functions are discussed from a more general point of view in Sec. 15.5-1; see also Sec. 9.4-3. Table 9.3-1 lists the Green’s functions for a number of boundary-value problems.

9.3-4. Reduction of Two-point Boundary-value Problems to Initial-value Problems. The general theory of boundary-value problems and eigenvalue problems involving ordinary differential equations is treated in Secs. 15.4-1 to 15.5-2 (see also Secs. 9.3-3, 20.9-2, and 20.9-3). The following method is often useful in connection with numerical solution methods.

Given an r^th-order linear differential equation Lw = f(z) with r suitable boundary conditions to be satisfied by w(z) and its derivatives for z = a, z = b, write the solution as

Where the w_k(z) are defined by the r + 1 initial-value problems

Apply the r given boundary conditions to the general solution (12a) to obtain r simultaneous equations for the r unknown coefficients α_k.

NOTE: Given a nonlinear boundary-value problem like

one can often calculate w(b) for two or three trial values of the unknown initial value w′(a); the correct value of w′(a) is then approximated by interpolation.

9.3-5. Complex-variable Theory of Linear Differential Equations. Taylor-series Solution and Effects of Singularities. (a) A given

Table 9.3-1. Green’s Functions for Linear Boundary-value Problems

Each boundary-value problem listed has the solution.

Use G(x, ξ) to obtain solutions for other initial or boundary conditions from Eq. (9.3-7) (see also Secs. 9.3-3, 9.4-3, 10.4-2, 15.4-8, and 15.5-1). The table yields solutions for other intervals (a, b) with the aid of suitable coordinate transformations.

^* This is a modified Green’s function in the sense of Sec. 15.5-1b and does not satisfy Eq.(9.3-10)

linear differential equation (1) has an analytic solution w = w(z) at every regular point z where the functions a_k(z) and f(z) are analytic (see also Sec. 7.3-3). If these functions are single-valued, and if D is a singly connected region of regular points, a given set of values w (z₀), w′(z_o), . . . , w^r-1(z₀) for some point z_o in D defines a unique solution w(z) in D.

To obtain this solution in Taylor-series form (see also Sec. 4.10-4), substitute

and

into the given differential equation (1); comparison of coefficients will yield recurrence relations for the unknown coefficients (k = r, r + 1, . . .). The series (14) converges absolutely and uniformly within every circle |z — z₀| < R in D.

(b) Analytic continuation of any solution w(z) of a linear differential equation around singularities of one or more coefficients a_k(z) will, in general, yield different branches of a multiple-valued solution (see also Secs. 7.4-2, 7.6-2, and 7.8-1).

In particular, one complete circuit around a singularity will transform a fundamental system of solutions w₁(z), W₂(z), . . . , w_r(z) of the homogeneous linear differential equation (2) into a new fundamental system . The two fundamental systems are necessarily related by a nonsingular linear transformation

The eigenvalues λ_k of the matrix [a_ik] (Sec. 13.4-2) are independent of the particular fundamental system W₁(z), W₂(z), . . . , w_r(z) in question.

9.3-6. Solution of Homogeneous Equations by Series Expansion about a Regular Singular Point. (a) A singularity z = Z₁ of one or more a_k(z) is an isolated singularity of the homogeneous linear differential equation (2) if and only if z₁ has a neighborhood containing no other singular point. An isolated singularity z = z₁ is a regular singular point of the homogeneous linear differential equation if and only if none of its solutions w(z) has an essential singularity at z₁; otherwise z = Z₁ is an essential singularity of the given differential equation. z = z₁ is a regular singular point if and only if a_k(z)/a₀(z) has at worst a pole of order k at (k = 1, 2, . . . , r) (Fuchs’s Theorem). In this case,

Eq. (2) can be rewritten as

where all p_k(z) are analytic in a neighborhood D₁ of z₁

It follows that a given homogeneous linear differential equation (2) admits a solution of the form

whenever z = Z₁ is a regular point or a regular singular point. The exponent µ must satisfy the r^th-degree algebraic equation

The first coefficient a₀ may be chosen at will, and the other coefficients a_k are found successively from a set of recurrence relations obtained on substitution of the series (17) into Eq. (2) or (16). The series converges absolutely and uniformly within every circle |z — z₁| < R in D₁.

Different roots µ = µ₁, µ₂, . . . , µ_r of the indicial equation (18) yield linearly independent solutions (17) of the given differential equation, unless two roots µk coincide or differ by an integer. In such cases, one may use the known solutions to reduce the order of the given differential equation in the manner of Sec. 9.3-2c or 9.3-7a, or use Frobenius’s method (Ref. 9.14); see also Sec. 9.3-8a.

The exponents µk are related to the eigenvalues λ obtained from Eq. (15) for one circuit about the regular singular point z₁ with λ_k = e^2πiµ_k.

(b)Regular Singular Points at Infinity (see also Secs. 7.2-1 and 7.6-3). z = ∞ is a regular singular point of Eq. (2) if and only if the transformation

yields a differential equation having a regular singular point at z̄ = 0. In this case, one may obtain solutions of the transformed equation in the manner of Sec. 9.3-6a.

(c) Generalization. If z = z₁ is not a regular singular point (e.g., if the functions P_k have poles at z = z₁ one may still write solutions similar to Eq. (17) by replacing each power series by a Laurent series (Sec. 7.5-3) admitting negative as well as positive powers of (z — z₁).

9.3-7. Integral-transform Methods. The solution of a linear differential equation (1) with polynomial coefficients and given initial conditions w⁽ⁱ⁾(0) = w₀(i) is often simplified through the use of the unilateral Laplace transformation in the manner of Sec. 9.4-5. Apply the formula

to obtain a new and possibly simpler differential equation for the Laplace transform £[w(z); s] of the solution. Boundary-value problems can be transformed to initialvalue problems by the method of Sec. 9.3-4.

More general integral transformations (Table 8.6-1) may be similarly employed in various special cases (see also Sec. 10.5-1 and Ref. 9.14).

9.3-8. Linear Second-order Equations(see also Secs. 15.4-3c and 15.5-4). (a) The theory of Secs. 9.3-1 to 9.3-4 applies to linear second-order equations, so that one is mainly interested in the solution of homogeneous linear second-order equations

Equation (21) is equivalent to

If any solution w₁(z) of Eq. (21) or (22) is known, the complete primitive is

(b) Series Expansion of Solutions (see also Secs. 9.3-5, 9.3-6, 9.3-9, and 9.3-10). In the important special case of a second-order equation expressible in the form

where p₁(z) and p₂(z) are analytic at z = z₁, the indicial equation (18) reduces to

The indicial equation has two roots. The root µ = µ₁ having the larger real part yields a solution of the form

The second root µ ₂ yields a similar linearly independent solution, which is replaced by

If µ₁ and µ₂ are identical or differ by an integer. Substitute each solution (26) or (27) into the given differential equation (24) to obtain recurrence relations for the coefficients.

(c)Transformation of Variables. The following transformations may simplify a given differential equation (21) or reduce it to a differential equation with a known solution.

1. w = w̄ exp ∫φ( z) dz yields

One attempts to choose φ(z) so as to simplify the new differential equation; in particular, φ(z) = — α₁(z)/2 eliminates the coefficient of dw̄/dz. A suitable substitution z = z(z̄) may.also yield a simpler differential equation.

2. The substitution

transforms Eq. (21) into a first-order differential equation of the Riccati type (Sec. 9.2-4c).

(d)Existence and Zeros of Solutions for Real Arguments. Let x be a real variable. The homogeneous linear differential equation

has a solution w = w(x) in every interval [a, b] where a₁(x) and a₂(x) are real and continuous; the solution is uniquely determined by the values w(x_o), w′(x₀) for some x₀ in [a, b]. Unless w ≡ 0 in [a, b], w (x) has at most a finite number of zeros in any finite interval [a, b]; the zeros of any two linearly independent solutions alternate in [a, b].

9.3-9. Gauss’s Hypergeometric Differential Equation and Riemann’s Differential Equation(see also Secs. 9.3-5 and 9.3-6). (a) The homogeneous linear differential equation

has regular singular points at z = ∞ (exponents µ = a, µ = b),z = 1 (exponents µ = 0, µ = c — a — b), and z = 0 (exponents µ = 0, µ = 1

— c), and no other singularities. The solutions of Eq. (31) include many elementary functions, as well as many of the special transcendental functions of Chap. 21 as special cases.

Series expansion about z = z₁ = 0 yields solutions (hypergeometric functions) (17) for µ = 0 and µ = 1 — c, with

For µ = 0 one obtains the special hypergeometric function

The series converges uniformly and absolutely for |z| < 1; the convergence extends to the unit circle if Re (a + b — c) < 1, except for the point z = 1 if Re (a + b — c) ≥ 0. The series reduces to a geometric series (Sec. 4.10-2) for a = 1, b = c, and to a Jacobi polynomial (Sec. 21.7-8) if a and/or b equals zero or any negative integer. The function (32) is undefined if one of the denominators c, c + 1, . . . equals zero and does not cancel out.

A second (linearly independent) solution of Eq. (31) may be obtained in the manner of Sec. 9.3-8b; in particular, the hypergeometric function of the second kind

is a solution whenever c is not an integer.

Note the following relations (|z| < 1):

The following formulas serve for analytic continuation of F(a, b; c; z) outside the unit circle:

See Table 9.3-1 and Refs. 21.9 and 21.11 for additional formulas.

(b)Riemann’s Differential Equation. Paperitz Notation. Equation (31) is a special case of the linear homogeneous differential equation

whose only singularities are distinct regular singular points at z = z₁ (exponents ∞, ∞′), z = z₂ (exponents β, β′), and z = z ₃ (exponents γ, γ′); note

A solution of Eq. (40), written in the so-called Paperitz notation, is

which reduces to Eq. (32) for z₁ = 0, z₂ = ∞, z₃ = 1 and α = 0, β = α, γ=0; α′ = 1 — c, β′ = b, γ′ = c — a — b. See Table 9.3-2 and Ref. 21.11 for additional formulas.

9.3-10. Confluent Hypergeometric Functions. One can move the singularity z = 1 of the hypergeometric differential equation (31) to z = b by substituting z/b for z; the singularity at z = b will then approach the original singularity at z = ∞ as b → ∞ (confluence of singularities). One thus obtains the new differential equation

whose only singularities are a regular singular point at z = 0 and an essential singularity at z = ∞. Many special transcendental functions are solutions of Eq. (40) for special values of a and c (Chap. 21).

Series expansion about z = z_x. = 0 yields solutions (confluent hypergeometric functions) (17) for ε = 0 and ε = 1 — c, with

Table 9.3-1. Additional Formulas Relating to Hypergeometric Functions

(Continued)

Table 9.3-2. Additional Formulas Relating to Confluent Hypergeometric Functions

For µ = 0 one obtains Kummer’s confluent hypergeometric function

(see also Sec. 21.7-5).

A second solution may be obtained in the manner of Sec. 9.3-86; in particular, the confluent hypergeometric function of the second kind

is a solution whenever c is not an integer.

Note the following relations:

See Refs. 21.9 and 21.11 for additional formulas.

9.3-11. Pochhammer’s Notation. The infinite series (32) and (43) are special cases of

where (x)k ≡ x (x + 1) . . . (x + k — 1). In this notation the hypergeometric function (32) becomes ₂F₁ (a, b, c; z), and the confluent hypergeometric function (43) is written as ₁F₁ (a; c; z).

9.4. LINEAR DIFFERENTIAL EQUATIONS WITH CONSTANT COEFFICIENTS

9.4-1. Homogeneous Linear Equations with Constant Coefficients (see also Sec. 9.3-1). (a) The first-order differential equation

has the solution

For a_o/a₁ > 0

a_o/a₁ is often referred to as the time constant.

(b)The second-order equation

has the solution

If a_o, a₁, and a₂ are real, s₁ and s₂ become complex for a₁² — 4a₀a₂ < 0; in this case, Eq. (4a) can be written as

where the quantities

are respectively known as the damping constant and the natural (characteristic) circular frequency. The constants C₁, C₂, A, B, R, and a are chosen so as to match given initial or boundary conditions (see also Sec. 9.4-5a).

If a₀a₂ > 0, the quantity is called the damping ratio; for ζ > 1, ζ = 1, 0 < ζ < 1 one obtains, respectively, an overdamped solution (4a), a critically damped solution (4b), or an underdamped (oscillatory) solution (4c). In the latter case, the logarithmic decrement 2πrσ₁/ω_N is the natural logarithm of the ratio of successive maxima of y(t).

Equation (3) is often written in the nondimensional form

with

is called the undamped natural circular frequency; for weak damping (ζ² << 1), ω ≈ ω_N (see also Fig. 9.4-1).

FIG. 9.4-1. Solution of the second-order differential equation

for y(0) = 0, dy/dt]_o = 1. Response is overdamped for fζ > 1, critically damped for ζ = 1, and underdamped for 0 < ζ < 1.

find the roots of the rth-degree algebraic equation

obtained, for example, on substitution of a trial solution y = e^st. If the r roots s₁, s₂, . . . of the characteristic equation (6) are distinct, the given differential equation (5) has the general solution

If a root s_k is of multiplicity m_k, replace the corresponding term in Eq. (7a) by

The various terms of the solution (7) are known as normal modes of the given differential equation. The r constants C_k and C_kj must be chosen so as to match given initial or boundary conditions (see also Sec. 9.4-5a).

If the given differential equation (5) is real, complex roots of the characteristic equation appear as pairs of complex conjugates σ ± iw. The corresponding pairs of solution terms will also be complex conjugates and may be combined to form real terms:

where A and B, or R and α, are new real constants of integration.

(d) Given a system of n homogeneous linear differential equations with constant coefficients

where the are polynomials in d/dt, each of the n solution functions y_k = y_k (t) (k = 1, 2, . . . , n) has the form (7); the s_k are now the roots of the algebraic equation

The constants of integration must again be matched to the given initial or boundary conditions (see also Secs. 9.4-5b and 13.6-2).

9.4-2. Nonhomogeneous Equations. Normal Response, Steadystate Solution, and Transients (see also Sec. 9.3-1). (a) The superposition theorems and solution methods of Secs. 9.3-1 to 9.3-4 apply to all linear ordinary differential equations. Thus the general solution of the nonhomogeneous differential equation

can be expressed as the sum of the general solution (7) of the reduced equation (5) and any particular integral of Eq. (10).

If, as in many applications, f(t) = 0 for t ≤ 0,* the particular integral y = y_N(t) of Eq. (10) with y_N = y′_N = y′′_N = . . . = y_N^(r—1) = 0 for t ≤ 0 will be called the normal response to the given forcing function f(t). To solve Eq. (10) for t> 0 with given initial values for y, y′, y′′ , . . . y^(r—1), one adds the solution of the corresponding initial-value problem for Eq. (5) to the normal response y_N(t).

In many applications (stable electric circuits, vibrations), all roots of the characteristic equation (9) have negative real parts, and the complementary function (7) dies out more or less rapidly (stable “transient solution”). In such cases, one is often mainly interested in a suitable nontransient particular integral y = yss(t), the “steady-state solution” due to the given forcing function f(t). In other cases, yss(t) is not uniquely defined by the given differential equation but depends on the initial conditions. The normal response yN(t) may or may not include a transient term.

(b) In the same manner, each solution function y_k = y_k(t) of a system of linear differential equations with constant coefficients,

can be expressed as the sum of the corresponding solution function of the complementary homogeneous system (8) and a particular solution function of the given system (11). The normal response of the system (11) to a set of forcing functions f_j(t) equal to zero for t ≤ 0 is the particular solution such that all y_k vanish for t ≤ 0 together with all derivatives which can be arbitrarily chosen (see also Sec. 13.6-2).

(c) If a forcing function contains a periodic term whose frequency equals that of an undamped sinusoidal term in the complementary function (7), then the differential equation or system may not have a finite solution (resonance, see also Secs. 9.4-5 and 15.4-12).

9.4-3. Superposition Integrals and Weighting Functions (see also Secs. 9.3-3, 9.4-7c, and 15.5-1). (a) Physically Realizable Initial-

* This means one considers only forcing functions of the type f(t) ≡ f(t)U+(t), where U+(t) is the asymmetrical unit-step function defined in Sec. 21.9-1 (see also the footnote to Sec. 9.4-3).

value Problems. Application of the Green's-function method of Sec. 9.3-3b to the differential equation (10) yields the normal-response solution (Sec. 9:4-2) as a weighted mean over past values of f(t) in the form

if one assumes (1) f(t) = 0 for t ≤ 0 (initial-value problem), and (2) h₊(t — τ) = 0 for t ≤ τ, so that “future” values of f(t) cannot affect “earlier” values of y(t), and “instantaneous” effects are also ruled out (physically realizable systems). More general problems are considered in Sec. 9.4-3d.

If the derivative on the right exists,

The weighting function h₊(t — τ) is the special Green’s function defined by

h₊(t — τ) is the normal response to an asymmetrical unit impulse δ₊(t — τ) (Sec. 21.9-6); Eq. (14a) may be rewritten as a “symbolic differential equation”

Note that dτ is the normal response to the asymmetrical unit-step function U₊(t — τ) (Sec. 21.9-1). The “symbolic differential equation” (14c) is often easily solved for h₊(t — τ) by the Laplace-transform method of Sec. 9.4-5 (see also Secs. 8.5-1 and 9.4-7); alternatively, h₊(t) can be found as that solution of the homogeneous differential equation

which satisfies the initial conditions

EXAMPLE: For Ly ≡ a(dy/dt) + y, one has.

(b) Under similar conditions, the normal-response solution of a system of linear differential equations (11) can be expressed in the form

where (h₊)_kj(t — τ) is the k^th solution function obtained on substitution of f_j (t) = δ₊(t — τ) (i = 1, 2, … , n) in Eq. (11). The weighting-function matrix [(h₊)_kj(t — τ)] is often called a state-transition matrix (Sec. 13.6-2).

(c) A superposition integral (12) also yields the normal response if the “input function” f(t) and the “output function” y(t) are related by a differential equation of the form

Such relations may result, in particular, if a system (11) is reduced to a single differential equation by elimination of all but one of the unknown functions.

Any (h₊)_kj(t), and also h₊ (t) if relations of the type (16) are considered, may contain delta-function-type singularities (modified Green’s function, Sec. 15.5-16). Thus, if h₊(t) = c₁δ₊(t — t₁) + c₂ δ ₊+(t — t₂) + . . . + h₀(t), then Eq. (12) yields the normal response

(d) More General Problems. “Symmetrical” vs. “Asymmetrical” Weighting Functions. To deal with forcing functions different from zero for t ≤ 0, one may introduce the “symmetrical” weighting function h(t — τ) defined by

or Lh(t — τ) = δ( t — τ) (see also Sec. 21.9-2) with suitable initial or boundary conditions. The resulting solution

will, in particular, satisfy Eq. (10) for t ≥ 0 with

If f( t ) ≡ f ( t ) U (t ) (Sec. 21.9-1), and one adds the condition

h₊(t) and h(t), and the solutions (12) and (18), are easily confused. The “asymmetrical” weighting function h₊(t) is particularly convenient for use with the unilateral Laplace transforms employed by most engineers, while h(t) fits the context of Fourier analysis or bilateral Laplace transforms (see also Sec. 18.10-5). In the usual physical applications, Eq. (19a) holds since “future” values of f(t) cannot affect the solution; h₊(t) and h (t) are then identical wherever they are continuous. Frequently, forcing functions cannot even affect the solution instantaneously, so that h(t — τ) satisfies the stronger condition

and h(t) and h₊ (t) are identical.

EXAMPLE: In a purely resistive electric circuit, the current y(t) and the voltage f(t) are related by y(t) = f (t)/R; so that h₊(t) = δ ₊(t)/R and h(t) = δ( t)/R. But for

9.4-4. Stability. A linear differential equation (10) or a system (11) will be called completely stable if and only if all roots of the corresponding characteristic equation (6) or (9) have negative real parts, so that effects of small changes in the initial conditions tend to zero with time (refer to Sec. 13.6-5 for a more general discussion of stability). The nature of the roots may be investigated with the aid of Secs. 1.6-6 and 7.6-9 (stability criteria for electric circuits and control systems). A differential equation (10) is completely stable if and only if [or equivalently exists; a similar condition for every weighting function of a system (11) is necessary and sufficient for complete stability of the system.

9.4-5. The Laplace-transform Method of Solution (see also Secs. 8.1-1, 8.4-1 to 8.4-5, 9.3-7, and 13.6-2). (a) To solve a linear differential equation (10) with given initial values y(0 + 0), y′(0 + 0), y′′(0 + 0), . . . , y^(r—1)(0 + 0), apply the Laplace transformation (8.2-1) to both sides, and let . The resulting linear algebraic equation (subsidiary equation)

is easily solved to yield the Laplace transform of the desired solution y(t) in the form

Here the first term is the Laplace transform Y_N(s) of the normal response yN(t) (Sec. 9.4-2a), and the second term represents the effects of nonzero initial values y(t) and its derivatives. The solutions y(t) and yN(t)are found as inverse Laplace transforms by reference to tables (Appendix D), or by one of the methods of Secs. 8.4-2 to 8.4-9. In particular, each of the r terms in the partial-fraction expansion of G(s)/(a_os^r + a₁s^r—1 + . . . + a_r) (Sec. 8.4-5) yields a corresponding term of the force-free solution (7).

This solution method applies without essential changes to differential equations of the type (16).

(b) In the same manner, one applies the Laplace transformation to a system of linear differential equations (11) to obtain

where the functions G_j(s) depend on the given initial conditions. The linear algebraic equations (22) are solved by Cramer’s rule (1.9-4) to yield the unknown solution transforms

where A_jk(s) is the cofactor of φ_jk(s) in the system determinant D(s) ≡ det [φ_jk(s)] (see also Sec. 1.9-2). The first sum in Eq. (23) is the Laplace transform of the normal-response solution, while the second sum represents the effect of the initial conditions.

The desired solutions y_k(t) are obtained from Eq. (23) by inverse Laplace transformation.

In problems involving unstable differential equations (Sec. 9.4-4) and/or impulse-type forcing functions, the solutions may contain delta-function-type singularities (see also Secs. 8.5-1 and 21.9-6).

9.4-6. Periodic Forcing Functions and Solutions. The Phasor Method. (a) Sinusoidal Forcing Functions and Solutions. Sinusoidal Steady-state Solutions. Every system of linear differential equations (11) with sinusoidal forcing functions of equal frequency,

admits a unique particular solution of the form

In particular, if all roots of the characteristic equation (9) have negative real parts (stable systems, Sec. 9.4-4), the sinusoidal solution (24b) is the unique steady-state solution obtained after all transients have died out (Sec. 9.4-2).

(b) The Phasor Method. Given a system of linear differential equations (11) relating sinusoidal forcing functions and solutions (24), one introduces a reciprocal one-to-one representation of these sinusoids by corresponding complex numbers (vectors, phasors)

The absolute value of each phasor equals the root-mean-square value of the corresponding sinusoid, while the phasor argument defines the phase of the sinusoid. The phasors (25) are related by the (complex) linear algebraic equations (phasor equations)

which correspond to Eq. (11) and may be solved for the unknown phasors

(see also Sec. 9.4-5b). In the case of resonance (Sec. 9.4-2c), the expression (27) may not exist (may become “infinitely large”).

(c) Rotating Phasors. A set of sinusoidal functions (24) satisfies any given system of linear differential equations (11) if and only if the same is true for the corresponding set of complex exponential functions (rotating phasors)

which are often more convenient to handle than the real sinusoids (24).

(d) More General Periodic Forcing Functions (see also Sees. 4.11-4, 4.11-5, and 9.4-5c). Given a stable system (11) with more general periodic forcing functions expressible in the form

one can apply the phasor method of Sec. 9.4-6a separately for each sinusoidal term and superimpose the resulting sinusoidal solutions to obtain the steady-state periodic solution. This procedure may be more convenient than the Laplace-transform method if only a few harmonics of the periodic solution are needed.

9.4-7. Transfer Functions and Frequency-response Functions. (a) Transfer Functions. The function

in Eq. (21) is known as a transfer function. The transfer function “represents” a linear operator (Sec. 15.2-7) which operates on the forcing function (input) to yield the normal response (output; see also Fig. 9.4-1).

FIG. 9.4-1. Transfer-function representation of linear differential equations with constant coefficients. If y_N(t) in turn serves as the forcing function for a second differential equation to produce the normal response z_N(t), the two transfer functions multiply, i.e., Z_N(s)/F(s) = H₁(s)H₂(s).

More generally, each function A_jk(s)/D(s) in Eq. (23) is the transfer function relating the normal-response “output” y_k(t) of the system (11) to the “input” f_j(t) when all other forcing functions vanish identically. The transfer functions A_jk(s)/D(s) together constitute the transfer matrix.

The transfer f unction corresponding to Eq. (16) is

(b)Frequency-response Functions (see also Sec. 9.4-6a). The frequency-response functions H(iω) and A_jk(iω)/D(iω) similarly relate the phasors representing sinusoidal forcing functions and steady-state solutions of given circular frequency ω. Specifically, the absolute value and the argument of a frequency-response function respectively relate the amplitudes and the phases of input and output sinusoid; thus for f(t) = B sin (ωt + β), y(t) = A sin (ωt + a)

If frequency-response functions are “cascaded” in the manner of Fig. 9.4-1, the amplitude responses |H(iω)| multiply, and the phase responses arg H(iω) add.

(c) Relations between Transfer Functions or Frequency-response Functions and Weighting Functions (see also Sees. 4.11-4e, 9.4-3, and the convolution theorem of Table 8.3-1). The transfer function H(s) is the unilateral Laplace transform of the asymmetrical weighting function h₊(t), and the bilateral Laplace transform (Sec. 8.6-2) of the symmetrical weighting function h(f):

Hence the frequency-response function H(iω) is related to the symmetrical weighting function h(t) by the Fourier transformation*

Equations (33) and (34) indicate the possibility of obtaining weighting functions as inverse Laplace or Fourier transforms of rational functions.

9.4-8. Normal Coordinates and Normal-mode Oscillations. (a) Free Oscillations. Small oscillations of undamped mechanical or electrical systems are often described by a set of n linear second-order differential equations of the form

where the matrices [a_jk] and [b_jk] are both symmetric, positive-definite (Sec. 13.5-2), and such that the resulting characteristic equation (9) has 2n distinct, nonzero, purely imaginary roots ±iω₁, ±iω ₂, . . . , ±iω_n. Pairs of these roots correspond to sinusoidal free oscillations at the n normal-mode frequencies ω₁/2π, ω ₂/2π, . . . , ω_n/2π .

One may introduce normal coordinates ̄y₁ ̄y₂, . . . , ̄y_n for the given system (35) by a linear transformation

with coefficients t_kh chosen in the manner of Secs. 13.5-5 and 14.8-7 so as to diagonalize the matrices [a_jk] and [b_jk] simultaneously; the transformed system takes the simple form

The resulting free sinusoidal normal-mode oscillations

do not affect one another (are “uncoupled”). The normal coordinates (38) may have an intuitive physical interpretation.

The problem is a generalized eigenvalue problem involving sets of n functions [ȳ₁(t), ȳ₂ (t), . . . , ȳ _n(t)] as eigenvectors (see also Sees. 13.6-2a, 14.8-7, and 15.4-5).

* See the footnote to Sec. 4.11-2.

EXAMPLE: For a pair of similar coupled oscillators described by

the normal coordinates are simply Given y₁ = 1, y₂ = dy₁/dt = dy₂/dt = 0 for t = 0, the normal-mode equations

For α² << ω₀² (weak coupling), this solution describes the so-called resonance phenomenon.

(b) Forced Oscillations. The corresponding forced-oscillation problem

can, in principle, be solved in the manner of Sec. 14.8-10 through normal-mode expansion of the forcing functions f_i(t). The Laplace-transform method of Sec. 9.4-5 is usually more convenient.

9.5. NONLINEAR SECOND-ORDER EQUATIONS

9.5-1. Introduction. Sections 9.5-2 to 9.5-5 introduce the general terminology and the most easily summarized approximation method of the theory of nonlinear oscillations. References 9.15 to 9.17, 9.22, and 9.23 are recommended for further study; for better or for worse, many solution methods are closely tied to specific applications.

The perturbation method of Sec. 10.2-7c is often used to simplify nonlinear problems, especially in celestial mechanics. See Secs. 20.7-4 and 20.7-5 for numerical methods of solution.

9.5-2. The Phase-plane Representation. Graphical Method of Solution (see also Sec. 9.2-2). The second-order differential equation

is equivalent to the system of first-order equations

The general solution y = y (t), of Eq. (1) or (2) can be represented geometrically by a family of directed phase-trajectory curves in the yy plane or phase plane. The phase-plane representation is most useful if the given function f(t, y, dy/dt) does not involve the independent variable t explicitly (e.g., “free“ oscillations). In this case, the system (2) is of the general form

and the phase trajectories satisfy the first-order differential equation

which specifies the slope of the solution curve through at that point. The resulting field of tangent directions (“phase-plane portrait” of the given differential equation) permits one to sketch ( y) and hence y(t) for given initial values of y and one may begin by drawing loci of constant slope dy/dy = m(isoclines, Fig. 9.5-1).

9.5-3. Critical Points and Limit Cycles (see also Sec. 9.5-4). (a) Ordinary and Critical Phase-plane Points. Given a differential equation (1) reducible to a system (3), a phase-plane point (y, ẏ) is an ordinary point if and only if P(y, ẏ) and Q(y, ẏ) are analytic and not both equal to zero; there exists a unique phase trajectory through each ordinary point. Phase-plane points (y₀, y_Q) such that

are critical points (or singular points) where the trajectory is not uniquely determined. Critical points are classified according to the nature of the phase trajectories in their neighborhood; Fig. 9.5-2 illustrates the most important types. Physically, critical points are equilibrium points admitting stable or unstable equilibrium solutions y = y_o (Sec. 9.5-4).

(b) Periodic Solutions and Limit Cycles. Periodic solutions y = y(t) correspond to closed phase-trajectory curves, and vice versa. A closed phase trajectory C is called a limit cycle if each trajectory point has a neighborhood of ordinary points in which all phase trajectories spiral into C (stable limit cycle, see also Sec. 9.5-4) or out of C (unstable limit cycle), or into C on one side of C and out of C on the other side (half-stable limit cycle). For an example, see Secs. 9.5-4c and 9.5-5 (see also Fig. 9.5-3).

FIG. 9.5-1. Isoclines, tangent directions, and some solutions of the differential equation

corresponding to Van der Pol’s differential equation

with dy/dt = , ε = 1. Only the right half-plane is shown.

FIG. 9.5-2. Phase trajectories in the neighborhood of six types of critical points (Sees. 9.5-3 and 9.5-4).

FIG. 9.5-3. (a) A stable limit cycle enclosing an unstable critical point at the origin. “Soft” self-excitation of oscillations for arbitrarily small initial values of y and y. (b) A stable limit cycle enclosing a stable critical point at the origin and an unstable limit cycle (shown in broken lines). “Hard” self-excitation of oscillations for initial values outside the unstable limit cycle.

(c) Poincare’s Index and Bendixson’s Theorems. The possible existence of limit cycles (stable oscillations) is of interest in many applications. In addition to the analytical criteria of Sec. 9.5-4, the following theory is sometimes helpful.

For any given phase-plane portrait, the index of a closed curve C containing only ordinary points is the number of revolutions of the solution tangent as the point (y, y) makes one complete cycle around C. The index of an isolated critical point P is the index of any closed curve enclosing P and no other critical point. Then

1. The index of any closed curve C equals the sum of the indices of all the (isolated) critical points enclosed by C; if C encloses only ordinary points, its index is zero.

2. The index of a nodal point, focal point, or vortex point is 1; the index of a saddle point is —1 (see also Fig. 9.5-2).

3. The index of every closed phase trajectory is 1; hence a limit cycle must enclose at least one critical point other than a saddle point (Fig. 9.5-3).

No closed phase trajectories exist within any phase-plane domain where ∂P/dy + ∂Q/dy is of one sign (Bendixson’s First Theorem). A trajectory which remains in a bounded region of the phase plane and does not approach any critical point for 0 ≤ t ≤ ∞ is either closed or approaches a closed trajectory asymptotically {Bendixson’s Second Theorem).

9.5-4. Poincaré-Lyapunov Theory of Stability (see also Sees. 13.6-5 to 13.6-7). (a) The solution y_i = y_1i(t) (i = 1, 2, . . . , n) of a system of ordinary differential equations reduced to the form

is stable in the sense of Lyapunov if and only if each y_i(t) → y_1i(t) as y_i(t_o) → y_1i(t₀) (i = 1, 2, . . . , n), the convergence being uniform for t > 0 (Sec. 4.4-4); sufficiently small initial-condition changes then cannot cause large solution changes. A stable solution is asymptotically stable if and only if there exists a bound r(t_Q) such that |y_i (t_o) < r(t_o)| < r(t₀) i implies for all i. An asymptotically stable solution is asymptotically stable in the large (completely stable) if and only if r(t₀) = ∞. This is, in particular, true for every solution of a completely stable system of linear differential equations with constant coefficients (all roots of characteristic equation have negative real parts, Secs. 9.4-4 to 9.4-7).

(b) Stability of Equilibrium. For a system of the form

with suitably differentiable f_i, an equilibrium solution

is asymptotically stable whenever the linearized system

is completely stable, where the partial derivatives are computed for y₁= y₁₁, y₂ = y₁₂, . . . , y_n = y_1n. This is true whenever all roots s of the characteristic equation

have negative real parts. The equilibrium is unstable if Eq. (9) has a root with positive real part; if the real parts of all roots are negative or zero, one requires a more detailed stability investigation (Sec. 13.6-7).

In particular, for a second-order differential equation (1) reducible to a system (3) the characteristic equation is

where all derivatives are computed at the equilibrium point (y₀, ₀). The equilibrium point is

A stable or unstable nodal point if both roots Si and s₂ of Eq. (10) are real and negative or positive, respectively

A saddle point if si and s₂ are real and of opposite sign

A stable or unstable focal point if si and s₂ are complex conjugates with negative or positive real parts, respectively

A vortex point if s₁ and s₂ are purely imaginary

(see also Sec. 9.5-3a and Fig. 9.5-2).

Examples for the six types of equilibrium points listed are most easily obtained from the linear differential equation (9.4-3) for different values of the coefficients (Sec. 9.4-16).

(c) Stability of a Periodic Solution. The stability of a periodic solution y = y_p(t) ≠ 0, y = y_p (t) of Eq. (3) depends on that of the linearized system

which is satisfied by small variations (Sec. 11.4-1) δy, δẏ of the periodic solution. This is a system of linear differential equations with periodic coefficients having the same period T as the given solution; Eq. (11) admits two linearly independent solutions of the form

where the h_ik(t) are periodic functions. The periodic solution is stable if

is less than zero, and unstable if λ > 0; the case λ = 0 requires additional investigation.

See Sec. 13.6-6 for additional stability criteria.

EXAMPLE: The differential equation

yields a stable limit cycle in the neighborhood of the approximate periodic solution y = y_p(t) ≈ 2 cos t, = _p(t) ≈ —2 sin t, with λ = —μ[1 + 0(μ)] (see also Sec. 9.5-5).

9.5-5. The Approximation Method of Krylov and Bogolyubov. (a) The First Approximation. Equivalent Linearization. To solve a differential equation of the form

where ω is a given constant, and the last term is a small nonlinear perturbation, write

Assuming that errors of the order of μ² are negligible, the “amplitude” r(t) and the “total phase” φ(t) are then obtained from the first-order differential equations

For a given initial value r (0) = r₀, the solution (9.4-4) of the equivalent linear differential equation

approximates the solution of the given differential equation (15) with an error of the order of μ². For a periodic solution such as a limit cycle (Sec. 9.5-3b), the approximate amplitude r_L is obtained from a₁(r_L) = 0, and the circular frequency is approximated by. The limit cycle is stable if da₁/dr]_r=rL > 0 and unstable if da₁/dr]_{r=r_L} < 0. For self-excitation from rest, one must have a₁(0) < 0.

The first approximation is of considerable interest in connection with periodic nonlinear oscillations. In such cases, the equivalent linear differential equation (18) yields the same energy storage and dissipation per cycle as the given nonlinear equation (15). The equivalent linear equation can therefore be used in many investigations of nonlinear resonance phenomena (Ref. 9.17).

(b) The Improved First Approximation. An improved first-order approximation is given by

where r (t) and φ( t) are given by Eq. (17), and

EXAMPLE: In the case of Van der Pol’s differential equation (14), Eq. (17) yields

There is a stable limit cycle for r = r_L= 2. The coefficients (20) all vanish except for β₃; the improved first approximation is

The Krylov-Bogolyubov approximation (19) is an improvement on an earlier method due to Van der Pol, who derived an approximation of the form

in an analogous manner.

The Krylov-Bogolyubov approximation method can be extended to apply in the case of a periodic forcing function on the right side of the differential equation (15) (nonlinear forced oscillations, subharmonic resonance, entrainment of frequency). For this and other methods, see Refs. 9.15 to 9.17.

9.5-6. Energy-integral Solution. Differential equations of the form

(9.5-23)

which are of considerable interest in dynamics, can be reduced to first-order equations through multiplication by dy/dt and integration:

9.6. PFAFFIAN DIFFERENTIAL EQUATIONS

9.6-1. Pfaffian Differential Equations (see also Sees. 3.1-16 and 17.2-2). A Pfaffian differential equation (first-order linear total differential equation)

with continuously differentiable coefficients P, Q, R may be interpreted geometrically as a condition P . dr= 0 on the tangent vector dr ≡ (dx, dy, dz) of a solution curve (integral curve) described by two equations f(x, y, z) = 0, g(x, y, z, C) = 0, where C is a constant of integration. To find the integral curves lying on an arbitrary regular surface

solve the ordinary differential equation obtained by elimination of z and dz from Eq. (1) and df(x, y, z) = 0.

9.6-2. The Integrable Case (see also Sec. 9.2-4). The Pfaffian differential equation (1) is integrable if and only if there exists an integrating factor μ = μ(x,y,z) such that μ(P dx + Q dy + R dz) is an exact differential dφ( x,y, z)>; this is true if and only if

In this case, every curve on an integral surface

orthogonal to the family of curves described by

is a solution. It follows that the solutions found in the manner of Sec. 9.6-1 from a conveniently chosen family of surfaces (usually planes)

lie on an integral surface (4) obtainable by elimination of X from the resulting solution f(x, y, z; λ) = 0, g(x, y, z, C; λ) = 0 (Mayer’s Method of Solution).

To find integral surfaces (4) by another method, hold z constant and obtain the solution of the ordinary differential equation P dx + Q dy = 0 in the form u(x, y, z) — K = 0. Then

describes an integral surface, where φ( z ) is the solution of the ordinary differential equation obtained by elimination of x and y from

with the aid of Eq. (7). Note that the function (8) is the required integrating factor ε( x , y , z ).

An important application is found in thermodynamics, where the reciprocal of the absolute temperature T is an integrating factor for the adiabatic condition δq = 0 of the form (1); δq / T is the (exact) differential of the entropy. See Ref. 10.23 for a discussion of total differential equations involving more than three variables.

The Laplace transformation Chap. 8

Calculus of variations Chap. 11

Matrix notation for systems of differential equations Chap. 13

Boundary-value problems and eigenvalue problems Chap. 15

Numerical solutions Chap. 20

Special transcendental functions Chap. 21

9.7-2. References and Bibliography.

9.1. Agnew, R. P.: Differential Equations, 2d ed., McGraw-Hill, New York, 1960.

9.2. Andronow, A. A., and C. E. Chaikin: Theory of Oscillations, Princeton University Press, Princeton, N.J., 1949.

9.3. Bellman, R.: Stability Theory of Differential Equations, McGraw-Hill, New York, 1953.

9.4. Bieberbach, L.: Differentialgleichungen, 2d ed., Springer, Berlin, 1965.

9.5. Birkhoff, G., and G. Rota: Ordinary Differential Equations, Blaisdell, New York, 1962.

9.6. Coddington, E. A.: An Introduction to Ordinary Differential Equations, Prentice-Hall, Englewood Cliffs, N.J., 1961.

9.7. and N. Levinson: Theory of Ordinary Differential Equations, McGraw- Hill, New York, 1955.

9.8. Ford, L. R.: Differential Equations, 2d ed., McGraw-Hill, New York, 1955.

9.9. Golomb, M., and M. E. Shanks: Elements of Ordinary Differential Equations, 2d ed., McGraw-Hill, New York, 1965.

9.10. Hale, J. K.: Oscillations in Nonlinear Systems, McGraw-Hill, New York, 1963.

9.11. Hartman, P.: Ordinary Differential Equations, Wiley, New York, 1964.

9.12. Hochstadt, H.: Differential Equations, Holt, New York, 1964.

9.13. Hurewicz, W.: Lectures on Ordinary Differential Equations, M.I.T., Cambridge, Mass., 1958.

9.14. Kamke, E.: Differentialgleichungen, Losungsmethoden und Losungen, vol. I, Chelsea, New York, 1948.

9.15. Krylov, N., and N. Bogolyubov: Nonlinear Oscillations, translated by S. Lefschetz, Princeton University Press, Princeton, N.J., 1943.

9.16. Lefschetz, S.: Differential Equations: Geometric Theory, 2d ed., Interscience, New York, 1963.

9.17. Minorski, N.: Nonlinear Oscillations, Van Nostrand, Princeton, N.J., 1962.

9.18. Petrovsky, I. G.: Lectures on Partial Differential Equations, Interscience, New York, 1955.

9.19. Pontryagin, L. S.: Ordinary Differential Equations, Addison-Wesley, Reading, Mass., 1962.

9.20. Saaty, T. L., and J. Bram: Nonlinear Mathematics, McGraw-Hill, New York, 1964.

9.21. Sansone, G., and R. Conti: Nonlinear Differential Equations, Macmillan, New York, 1964.

9.22. Stoker, J. J.: Nonlinear Vibrations, Interscience, New York, 1950.

9.23. Struble, R.: Nonlinear Differential Equations, McGraw-Hill, New York, 1961.

9.24. Tenenbaum, M., and H. Pollard: Ordinary Differential Equations, Harper and Row, New York: 1963.

9.25. Tricomi, F. G.: Differential Equations, Hafner, New York, 1961.