Stationary solutions of random differential equations with polynomial nonlinearities

J. vom Scheidt  H.-J. Starkloff  R. Wunderlich

Abstract

The paper deals with systems of ODEs containing polynomial nonlinearities and random inhomogeneous terms. Applying perturbation method pathwise solutions are found in form of power series with respect to a parameter η controlling the nonlinearities. Under the assumption that for η = 0 the system is stable and that the inhomogeneous terms are bounded the radius of convergence of the perturbation series is estimated. Further, it is proved that the perturbation series form stationary solutions if the inhomogeneous terms are stationary.

1 Introduction

Let (Ω, A, P) be a probability space, where A denotes the σ-algebra of subsets of Ω on which is defined a probability measure P. Further, let z(t, ω) be a random function defined on R × Ω with values in C^{n}, n ≥ 1, where R and C denote the sets of real and complex numbers, respectively. The function z is required to satisfy a system of nonlinear first-order ODE

ż = Az + ηd(z) + Bh.  (1)

Here, a n × n matrix A with constant complex entries is involved in the linear part while the nonlinearities are described by a polynomial vector function d : C^{n} → C^{n} possessing the components

d_{i}(z) = ∑_{k=2}^{κ} ∑_{i_{1},...,i_{k}=1}^{n} d_{ki,i_{1},...,i_{k}}z_{i_{1}} · ... · z_{i_{k}}, i = 1, ..., n  (2)

with κ ≥ 2, some complex coefficients d_{ki,i_{1},...,i_{k}} and a non-negative parameter η. The inhomogeneous term contains the random function h with continuous paths defined on R × Ω with values in C^{r}, r ≥ 1, and a n × r matrix B with constant complex entries.

The present paper deals with the determining of a solution z of (1) to given A, B, η, d and h and the question whether there exists a stationary solution if the random function h is stationary. A random function z(t, ω) is said to be a stationary solution of (1) if z satisfies Equation (1) and if (z, h) is a stationary random function, i.e. z and h are stationarily related.

This problem arises e.g. in the investigation of the long time behaviour of the response of discrete vibration systems with a stationary random external excitation (see Soong, Grigoriu [13] and [14, 15, 16]). In the linear case, i.e. for η = 0, stationary solutions exist if the matrix A is stable, i.e. all eigenvalues of A possess strictly negative real parts (see Arnold, Wihstutz [2], Bunke [3], p. 45ff). In the present nonlinear case additional conditions to the form of nonlinearities and to the distribution of the random function h has to be taken into account.

The existence of a stationary solution (as well as of solutions with a periodic distribution) is investigated in a number of papers on qualitative theory and stability of stochastic systems, see e.g. Khas’minskij [8], Arnold, Kliemann [1], Bunke [3], Dorogovtsev [4, 5], Ito, Nisio [7]. A typical condition to the nonlinearities is global Lipschitz continuity of d(.) which is not fulfilled in the present case of a polynomial function d because of unbounded partial derivatives of d . Another approach is based on the existence of a bounded solution z.

The starting point of our approach is the consideration of pathwise solutions and the treatment of the nonlinear term ηd(z) as a perturbation of the linear term Az. Solutions of (1) are expressed as power series with respect to the parameter η, i.e.

z = ∑_{p=0}^{∞}^{p}z η^{p}.  (3)

The coefficients ^{p}z are determined in Section 2. Section 3 gives conditions for the convergence of the above perturbation series and estimates the radius of convergence η_{0}. Beside the stability of the matrix A, the essential condition used here is the boundedness of the function h. It is shown that for η < η_{0} the resulting perturbation series is a pathwise solution of system (1). As a by-product it results the pathwise boundedness of the solutions found in form of a perturbation series. In Section 4 it is proved that these solutions are stationary if h is a stationary random function.

The representation of pathwise solutions of (1) by perturbation series can be used for determining the distribution law and moment functions of z. An approximative method based on expansions with respect to η and to the correlation length of weakly correlated random functions which are involved in the random function h is described in [14], [15], [16].

2 Perturbation method

In this and in the next section Eq. (1) is considered for fixed ω ∈ Ω. For the sake of a short notation the dependence of h and z on ω is omitted during this “pathwise analysis”.

In order to find pathwise solutions of (1) for given A, η, d, B and a random function h with continuous paths the function z is represented as a power series (3) with respect to the parameter η. The coefficients ^{0}z, ^{1}z,... can be found by substituting series (3) into Equation (1) and equating the coefficients of the powers of η. First, this procedure is carried out formally. A verification of the results is given in the next section.

First, the nonlinear term d(z) is expressed as power series with respect to η. It holds for i_{1}, ..., i_{k} = 1, ..., n, k = 2, ..., κ

z_{i_{1}} · ... · z_{i_{k}} = ( ∑_{p_{1}=0}^{∞}^{p_{1}}z_{i_{1}}η^{p_{1}} ) · ... · ( ∑_{p_{k}=0}^{∞}^{p_{k}}z_{i_{k}}η^{p_{k}} )

= ∑_{p_{1},...,p_{k}=0}^{∞}^{p_{1}}z_{i_{1}} · ... ·^{p_{k}}z_{i_{k}}η^{p_{1}+...+p_{k}}
= ∑_{p=0}^{∞} η^{p} ∑_{p_{1}+...+p_{k}=p}^{p_{1}}z_{i_{1}} · ... ·^{p_{k}}z_{i_{k}}.

In the above derivation and everywhere below it is assumed that the indices p_{1}, ..., p_{k} in sums of the type ∑_{p_{1}+...+p_{k}=p} ... are nonnegative.

From representation (2) it follows for the components of d(z)

d_{i}(z) = ∑_{p=0}^{∞}^{p}d_{i}(^{0}z, ...,^{p}z) η^{p}, i = 1, ..., n,

(4)

with ^{p}d_{i}(^{0}z, ...,^{p}z) = ∑_{k=2}^{κ} ∑_{i_{1},...,i_{k}=1}^{n} d_{ki,i_{1},...,i_{k}} ∑_{p_{1}+...+p_{k}=p}^{p_{1}}z_{i_{1}} · ... ·^{p_{k}}z_{i_{k}}.

Substituting the series (3) and (4) into Equation (1) it yields

∑_{p=0}^{∞}^{p}ż η^{p} = ∑_{p=0}^{∞} A^{p}z η^{p} + η ∑_{p=0}^{∞}^{p}d(^{0}z, ...,^{p}z) η^{p} + Bh.

For the coefficients ^{p}z results an infinite sequence of linear first-order systems

pż = A^{p}z +^{p}b, p ≥ 0,
with ^{p}b = \{ \begin{matrix}Bh & for p = 0 \\ ^{p−1}d(^{0}z, ...,^{p−1}z) & for p > 0. \end{matrix}

The first order systems above possess solutions

pζ(t) = ∫_{−∞}^{t} e^{A(t−s)} ^{p}b(s)ds = ∫_{0}^{∞}e^{As} ^{p}b(t − s)ds  (5)

which can be found recursively. Now, the power series

ζ(t) := ∑_{p=0}^{∞}^{p}ζ(t) η^{p}  (6)

is considered. In Section 3 conditions for the convergence of the above series are investigated and in Section 4 it is checked if ζ(t) satisfies the system (1).

3 Convergence of perturbation series

After the investigation of the coefficients of the perturbation series conditions for the convergence of series (6) will be determined.

Stability for the linear part of the system (1) is assumed, i.e. the matrix A possesses eigenvalues with strictly negative real parts, only. The diagonalizability of A is supposed as an additional technical condition. Further, the function h is assumed to be bounded and continuous on R. Below, the notations |X| = {|x_{ij}|}_{ij} and [X]_{ij} = x_{ij} for a matrix X = {x_{ij}}_{ij} will be used.

Theorem 1 If the matrix A is stable and diagonalizable, i.e. it exists a representation A = VΛV^{−1} with Λ = diag(λ_{1}, ..., λ_{n}), Re [λ_{i}] < 0, i = 1, ..., n and if for the components of h it holds |h_{i}(t)| ≤ H, i = 1, ..., r, ∀ t ∈ R, then the series ∑_{p=0}^{∞}^{p}ζ(t) η^{p} with coefficients ^{p}ζ(t) given in (5) converges uniformly

with respect to t ∈ R for η < η_{0} := \frac{λ̂}{KD}. Here, the following notations

λ̂ := min_{1≤i≤n} {|Re [λ_{i}]|}
K := (κ − 1) (\frac{κ}{κ − 1})^{κ}
D := max_{1≤i≤n} \{V̂ ∑_{k=2}^{κ} (\frac{V̂ Ĥ}{λ̂} )^{k−1} ∑_{i_{1},...,i_{k}=1}^{n} |d_{ki,i_{1},...,i_{k}}|\}

V̂ := max_{1≤i≤n} {∑_{j=1}^{n} [|V| |V^{−1}|]_{ij}}

Ĥ := max_{1≤i≤n} {∑_{j=1}^{r} |B_{ij}|} H

have been used.

Proof: In a first step by means of mathematical induction it is proved that the components of ^{p}ζ(t) are bounded by

|^{p}ζ_{i}(t)| ≤ _{κ}c_{p} \frac{V̂ Ĥ}{λ̂} (\frac{D}{λ̂} )^{p} ∀ t ∈ R, i = 1, ..., n.  (7)

Here, (_{κ}c_{p})_{p≥0} denotes the integer sequence of generalized CATALAN numbers of order κ (see Sloane, Plouffe [12], Hilton, Pederson [6], Ledermann et.al. [9], p. 93ff. ) ^{1} defined by the recursion

κc_{0} = 1
κc_{p+1} = ∑_{p_{1}+...+p_{κ}=p}_{κ}c_{p_{1}} · ... ·_{κ}c_{p_{κ}}, p ≥ 0.

For p = 0 it holds ^{0}ζ(t) = ∫_{0}^{∞}e^{Au}Bh(t − u)du, i.e. for i = 1, ..., n it yields

|^{0}ζ_{i}(t)| = |∫_{0}^{∞}∑_{j=1}^{n} [e^{Au}]_{ij} [Bh(t − u)]_{j} du|

≤ ∑_{j=1}^{n} ∫_{0}^{∞}|[e^{Au}]_{ij}| |[Bh(t − u)]_{j}| du

≤ Ĥ ∑_{j=1}^{n} ∫_{0}^{∞}|[e^{Au}]_{ij}| du.  (8)

Using the eigenvalue decomposition A = VΛV^{−1} and the assumed property of negative real parts of the eigenvalues λ_{1}, ..., λ_{n} it follows

∫_{0}^{∞}|e^{Au}| du = ∫_{0}^{∞}|Ve^{Λu}V^{−1}| du ≤ |V| ∫_{0}^{∞}|e^{Λu}| du |V^{−1}|

= |V| diag (\frac{1}{−Re [λ_{1}]}, ..., \frac{1}{−Re [λ_{n}]}) |V^{−1}|
≤ \frac{1}{λ̂} |V| |V^{−1}|

and

∑_{j=1}^{n} ∫_{0}^{∞}|[e^{Au}]_{ij}| du ≤ \frac{1}{λ̂} ∑_{j=1}^{n} [|V| |V^{−1}|]_{ij} ≤ \frac{V̂}{λ̂} .  (9)

Applying inequalities (9) and (8) it results

|^{0}ζ_{i}(t)| ≤ \frac{ĤV̂}{λ̂} = _{κ}c_{0} \frac{V̂ Ĥ}{λ̂} (\frac{D}{λ̂} )^{0} ∀ t ∈ R, i = 1, ..., n.

1The authors want to express their thanks to Dr. Hans-Jürgen Fischer for his hints and helpful discussions on Catalan numbers.

Assuming the assertion (7) is valid for q ≤ p the assertion for q = p + 1 will be proved. It holds

p+1ζ(t) = ∫_{0}^{∞}e^{Au} ^{p}d (^{0}ζ(t − u), ...,^{p}ζ(t − u)) du

and for the components ^{p+1}ζ_{i}(t) it follows for i = 1, ..., n

|^{p+1}ζ_{i}(t)| = |∫_{0}^{∞}∑_{j=1}^{n} [e^{Au}]_{ij} ^{p}d_{j} (^{0}ζ(t − u), ...,^{p}ζ(t − u)) du|

≤ ∑_{j=1}^{n} ∫_{0}^{∞}|[e^{Au}]_{ij}| |^{p}d_{j} (^{0}ζ(t − u), ...,^{p}ζ(t − u))| du. (10)

Using representation (4) and the relation (7) for q ≤ p the terms |^{p}d_{j}(.)|, j = 1, ..., n in (10) can be estimated as follows

|^{p}d_{j} (^{0}ζ, ...,^{p}ζ)| =|∑_{k=2}^{κ} ∑_{i_{1},...,i_{k}=1}^{n} d_{kj,i_{1},...,i_{k}} ∑_{p_{1}+...+p_{k}=p}^{p_{1}}ζ_{i_{1}} · ... ·^{p_{k}}ζ_{i_{k}} |

≤ ∑_{k=2}^{κ} ∑_{i_{1},...,i_{k}=1}^{n} |d_{kj,i_{1},...,i_{k}}| ∑_{p_{1}+...+p_{k}=p}_{κ}c_{p_{1}} \frac{V̂ Ĥ}{λ̂} (\frac{D}{λ̂} )^{p_{1}} · ... ·_{κ}c_{p_{k}} \frac{V̂ Ĥ}{λ̂} (\frac{D}{λ̂} )^{p_{k}}

= (\frac{D}{λ̂} )^{p} \frac{Ĥ}{λ̂} V̂ ∑_{k=2}^{κ} (\frac{V̂ Ĥ}{λ̂} )^{k−1} ∑_{i_{1},...,i_{k}=1}^{n} |d_{kj,i_{1},...,i_{k}}| ∑_{p_{1}+...+p_{k}=p}_{κ}c_{p_{1}} · ... ·_{κ}c_{p_{k}}

≤ (\frac{D}{λ̂} )^{p} \frac{Ĥ}{λ̂} D _{κ}c_{p+1} = (\frac{D}{λ̂} )^{p+1} Ĥ _{κ}c_{p+1}.  (11)

For the above inequality the estimation

∑_{p_{1}+...+p_{k}=p}_{κ}c_{p_{1}} · ... ·_{κ}c_{p_{k}} = ∑_{p_{1}+...+p_{k}=p}_{κ}c_{p_{1}} · ... ·_{κ}c_{p_{k}} ∏_{i=1}^{κ−k} 1
≤ ∑_{p_{1}+...+p_{κ}=p}_{κ}c_{p_{1}} · ... ·_{κ}c_{p_{k}} ·_{κ}c_{p_{k+1}} · ... ·_{κ}c_{p_{κ}} = _{κ}c_{p+1}

has been used which holds for k = 2, ..., κ. Applying estimations (9), (10) and (11) it follows

|^{p+1}ζ_{i}(t)| ≤ \frac{V̂ Ĥ}{λ̂} (\frac{D}{λ̂} )^{p+1} _{κ}c_{p+1} ∀ t ∈ R

and the assertion (7) is proved.

Next, upper bounds for the generalized Catalan numbers _{κ}c_{p} involved in the estimates of |^{p}ζ_{i}(t)| are derived. The explicit representation of the p-th term of the sequence is (see Ledermann et.al. [9], p. 93 ff, Hilton, Pederson [6])

κc_{p} = \frac{(κp)!}{p! ((κ − 1) p + 1)!} p ≥ 0,

which for p > 0 can be written as _{κ}c_{p} = \frac{1}{p} \binom{κp}{p−1}. From the above representation for the ratio of two subsequent terms _{κ}c_{p} and κc_{p+1}, p ≥ 0, it follows

\frac{_{κ}c_{p+1}}{_{κ}c_{p}} = \frac{(κ(p + 1))!}{(p + 1)!((κ − 1)(p + 1) + 1)!} \frac{p!((κ − 1)p + 1)!}{(κp)!}
= \frac{1}{p + 1} \frac{(κ(p + 1))!}{(κp)!} \frac{((κ − 1)p + 1)!}{((κ − 1)(p + 1) + 1)!}
= \frac{1}{p + 1} \frac{(κp + κ)!}{(κp)!} \frac{(κp − p + 1)!}{(κp − p + κ)!}
= \frac{1}{p + 1} \frac{κp + 1}{κp − p + 2} · \frac{κp + 2}{κp − p + 3} · ... · \frac{κp + κ − 1}{κp − p + κ} · (κp + κ)
= \frac{κp + κ}{p + 1} ∏_{l=1}^{κ−1} \frac{κp + l}{κp − p + 1 + l}
= κ ∏_{l=1}^{κ−1} \frac{κp + l − (p − 1) + (p − 1)}{κp + l − (p − 1)} =κ ∏_{l=1}^{κ−1} (1+ \frac{p − 1}{p(κ − 1)+1 + l} )
≤ κ ∏_{l=1}^{κ−1} (1 + \frac{p − 1}{(p − 1)(κ − 1)})
= κ (1 + \frac{1}{κ − 1})^{κ−1}= κ (\frac{κ}{κ − 1})^{κ−1}= (κ − 1) (\frac{κ}{κ − 1})^{κ} = K.

Using _{κ}c_{0} = 1 and the derived inequality \frac{_{κ}c_{p+1}}{_{κ}c_{p}} ≤ K the estimate

κc_{p} = ∏_{l=0}^{p−1} \frac{_{κ}c_{l+1}}{_{κ}c_{l}} ≤ K^{p}  (12)

can be found. Applying inequalities (7) and (12) it results that the perturbation series

∑_{p=0}^{∞}^{p}ζ_{i}(t) η^{p} are majorized by ∑_{p=0}^{∞} \frac{V̂ Ĥ}{λ̂} (\frac{KDη}{λ̂} )^{p} for i = 1, ..., n. Since the ma-

jorizing series converges for \frac{KDη}{λ̂} < 1 a sufficient condition for the uniform convergence of perturbation series (6) is η < η_{0} = \frac{λ̂}{KD}.

Corollary 1 Under the assumption of Theorem 1 the perturbation series (6) is bounded on R and it holds

|ζ_{i}(t)| = |∑_{p=0}^{∞}^{p}ζ_{i}(t) η^{p}| ≤ \frac{V̂ Ĥ}{λ̂ − KDη}

for η < η_{0} = \frac{λ̂}{KD}, i = 1, ..., n.

In the special case of a scalar equation (1), i.e. n = 1, which is written as

ż = az + ηd(z) + bh where d(z) = ∑_{k=2}^{κ} d_{k}z^{k}

the following assertion can be derived. Corollary 2 For Re [a] < 0 and |h(t)| ≤ H, ∀ t ∈ R, the perturbation series ∑_{p=0}^{∞}^{p}ζ(t) η^{p} converges uniformly with respect to t ∈ R for

η < η_{0} = \frac{|Re [a]|}{K ∑_{k=2}^{κ} |d_{k}| (\frac{|b|H}{|Re [a]|})^{k−1}} .

If the nonlinearity is a single power function, i.e. d(z) = d_{κ}z^{κ}, κ ≥ 2 then the series converges uniformly on R for

η < η_{0} = \frac{|Re [a]|^{κ}}{K|d_{κ}| (|b|H)^{κ−1}} .

After proving the convergence of the series ζ = ∑_{p=0}^{∞}^{p}ζ(t) η^{p} for η < η_{0} = \frac{λ̂}{KD} it remains to check, if ζ satisfies Equation (1). Theorem 2 Under the assumption of Theorem 1 the perturbation series ζ = ∑_{p=0}^{∞}^{p}ζ(t) η^{p} with coefficients ^{p}ζ given in (5) satisfies Equation (1) if η < η_{0} where η_{0} is given in Theorem 1. Proof: First, it is proved that the representation ζ̇ = ∑_{p=0}^{∞}^{p}ζ̇(t) η^{p} is valid for η < η_{0}. To this end the uniform convergence of the formal differentiated series \frac{d}{dt} ( ∑_{p=0}^{∞}^{p}ζ(t) η^{p}) for η < η_{0} is checked. Using representation (5) of the coefficients ^{p}ζ formal differentiation leads to

\frac{d}{dt}(∑_{p=0}^{∞}^{p}ζ(t) η^{p}) = ∑_{p=0}^{∞}^{p}ζ̇(t) η^{p}
= A ^{0}ζ + Bh +∑_{p=1}^{∞} A ^{p}ζ(t) η^{p} +∑_{p=1}^{∞}^{p−1}d(^{0}ζ, ...,^{p−1}ζ)η^{p}.

The uniform convergence of the first series on the right hand side for η < η_{0} follows immediately from Theorem 1.

Now the second series is investigated. From (11) and (12) it follows for the components ^{p−1}d_{i}(^{0}ζ, ...,^{p−1}ζ), i = 1, ..., n,

|^{p−1}d_{i}(^{0}ζ, ...,^{p−1}ζ)| ≤ (\frac{D}{λ̂} )^{p} Ĥ _{κ}c_{p} ≤ (\frac{KD}{λ̂} )^{p} Ĥ.

It yields, that ∑_{p=1}^{∞}^{p−1}d(^{0}ζ, ...,^{p−1}ζ)η^{p} is majorized by Ĥ ∑_{p=1}^{∞} (\frac{KDη}{λ̂} )^{p} which con-

verges for |\frac{KDη}{λ̂} | < 1. From Weierstrass majorant criterion it follows the uni-

form convergence of ∑_{p=1}^{∞}^{p−1}d(^{0}ζ, ...,^{p−1}ζ)η^{p} for η < η_{0} = \frac{λ̂}{KD}.

The perturbation series ζ(t) = ∑_{p=0}^{∞}^{p}ζ(t) η^{p} satisfies equation (1) for η < η_{0}

since it holds

ζ̇ = A^{0}ζ + Bh + ∑_{p=1}^{∞} (A^{p}ζ +^{p−1}d(^{0}ζ, ...,^{p−1}ζ)) η^{p}
= Aζ + ηd(ζ) + Bh

where for d(ζ) the representation (4) was used.

4 Stationary solutions

The preceding results concerning the perturbation series representation of solutions of the nonlinear system (1) for fixed ω ∈ Ω can be used for the stochastic analysis of solutions with respect to random functions h_{i}(t, ω) with continuous and bounded paths. Analogously, it follows that a pathwise solution is given by

ζ(t, ω) = ∑_{p=0}^{∞}^{p}ζ(t, ω) η^{p}

for η < η_{0} where the random coefficients are found by (see (5))

0ζ(t, ω) = ∫_{0}^{∞}e^{As} Bh(t − s, ω)ds

pζ(t, ω) = ∫_{0}^{∞}e^{As} ^{p−1}d(^{0}ζ(t − s, ω), ...,^{p−1}ζ(t − s, ω))ds

for p > 0.

Next, the question whether ζ(t, ω) is stationary if h(t, ω) is a stationary random function is investigated. This question is related to the problem of the existence of stationary solutions of nonlinear equations with a stationary random exitation term which has been investigated in the context of stability of stochastic systems e.g. in Khas’minskij [8], Bunke [3], Arnold, Kliemann [1], Ito, Nisio [7]. Here, explicit representations of the solutions need not to be known. Some of these results in literature show the existence of stationary solutions for nonlinearities d(z) which are global Lipschitz-continuous or possess bounded partial derivatives. This assumption is not fulfilled in the present case of polynomial nonlinearities.

Another type of results proves the existence of stationary solutions under the assumption of the boundedness of solutions. In Khas’minskij [8], p. 52 ff the condition

∫_{0}^{T} P(|z(t, ω)| > R)dt → 0 as R → ∞

uniformly in T > T_{o} > 0 or T < T_{0} < 0 is used while Ito, Nisio [7] require that certain moments of the solution are bounded. These results correspond to the present case since from Corollary 1 it can be deduced that for η < η_{0} it holds

|ζ_{i}(t, ω)| ≤ \frac{V̂ Ĥ}{λ̂ − KDη} .

Provided, ζ(t, ω) is a stationary random function not only the existence but also the construction of a stationary solution (in form of a perturbation series) would be clarified. To answer the question of stationarity of ζ the assertions of the following lemmas will be used.

Lemma 1 Bunke [3], p. 38, Dorogovtsev [4] Let ϕ : R×Ω → R^{n}, n ≥ 1, be a stationary random function and J : R^{n} → R^{m}, m ≥ 1, a measurable mapping. Then the composition (ϕ, J (ϕ)) is stationary.

Lemma 2 Let ϕ : R × Ω → R^{n}, n ≥ 1, be a stationary random function with continuous paths and J : R^{n} → R^{m}, m ≥ 1, a measurable mapping. Further, let A be a m×m-matrix with eigenvalues possessing strictly negative real parts and let

ψ(t, ω) := ∫_{0}^{∞}e^{As} J (ϕ(t − s, ω)) ds.

Then (ϕ, ψ) is stationary.

Proof: The proof is analogous to the proofs of theorems 3.4 and 3.5 in Bunke [3] where the special case J : R^{n} → R^{n}, J(x) = x, is considered.

Theorem 3 The nonlinear system (1) possesses a stationary solution which can be represented as perturbation series

ζ(t, ω) = ∑_{p=0}^{∞}^{p}ζ(t, ω) η^{p}

with coefficients ^{p}ζ given in (5) if the following assumptions are fulfilled: i) the matrix A is diagonalizable and stable , i.e. it exists an eigenvalue

decomposition A = VΛV^{−1} with Λ = diag(λ_{1}, ..., λ_{n}) and Re [λ_{i}] < 0, i = 1, ..., n,

ii) the random function h(t, ω) is stationary and possesses continuous and bounded paths, it holds |h_{i}(t, ω)| ≤ H, a.s. ∀ t ∈ R,

iii) for the parameter η it holds η < η_{0}, where η_{0} is given in Theorem 1.

Proof: From Theorem 2 it follows that under assumptions (i), (ii), (iii) the function ζ(t, ω) is a solution of system (1). It remains to prove that (h, ζ) is a stationary random function.

a) By means of mathematical induction it can be proved that ∀ N ≥ 0 (h,^{0}ζ, ...,^{N}ζ) is stationary.

For p = 0 the coefficient ^{p}ζ(t, ω) is given by (see (5))

0ζ(t, ω) = ∫_{0}^{∞}e^{As} Bh(t − s, ω)ds.

Applying Lemma 2 with ϕ = h, J(ϕ) = Bϕ and ψ = ^{0}ζ it follows that (h,^{0}ζ) is stationary.

Now, assuming the stationarity of (h,^{0}ζ, ...,^{N}ζ) it will be proved that (h,^{0}ζ, ...,^{N+1}ζ) is stationary. From relation (5) it follows

N+1ζ(t, ω) = ∫_{0}^{∞}e^{As} ^{N}d(^{0}ζ(t − s, ω), ...,^{N}ζ(t − s, ω))ds.

Again, the application of Lemma 2 with

ϕ = (h,^{0}ζ, ...,^{N}ζ), J(ϕ) = ^{N}d(^{0}ζ, ...,^{N}ζ), ψ = ^{N+1}ζ,

shows the stationarity of (h,^{0}ζ, ...,^{N}ζ,^{N+1}ζ).

b) Now, the measurable mapping J : R^{(N+2)n} → R^{2n} defined by

J(\overline{x}, x_{0}, …, x_{N}) = (\overline{x}, ∑_{p=0}^{N} x_{p}η^{p}) , with \overline{x}, x_{0}, …, x_{N} ∈ R^{n}, η > 0,

is considered. From Lemma 1 it follows that for all N ≥ 0 the composition

((h,^{0}ζ, ...,^{N}ζ), J(h,^{0}ζ, ...,^{N}ζ))

is stationary. From the assertion of Theorem 1 it is known that in case

of η < η_{0} the series ∑_{p=0}^{N}^{p}ζ(t, ω)η^{p} converges a. s. and uniformly in t as

N → ∞. Consequently, the limit

lim_{N→∞} (h, ∑_{p=0}^{N}^{p}ζη^{p}) = (h, ζ)

is a stationary random function.

Remark The assertions of Lemma 1 and 2 hold analogously, if the term ”stationary” is replaced by ”periodically distributed with period T ”, see e.g. Bunke [3], Dorogovtsev [4, 5]. Moreover, an analogous assertion of Theorem 3 can be proved. That means, if the assumption of stationarity of the inhomogeneous term h(t, ω) is replaced by the property of periodicity of its distribution then the perturbation series (3) form solutions of Equation (1) which are periodically distributed, too. This seems to be useful in the analysis of equations whose inhomogeneous terms contain stationary weakly correlated random functions (see e.g. [14, 15, 16]). The stochastic simulation of those equations requires the replacement of stationary by periodically distributed weakly correlated functions, see the corresponding chapters on simulation processes in [10], p. 53 ff and [11].

5 Randomly forced nonlinear oscillator

The results of this paper are applied to the following model of the motion of a single oscillator with linear damping, nonlinear stiffness and an external random excitation

ÿ + 2γẏ + β^{2}y + ηy^{κ} = h(t, ω).

Here, the scalar random excitation function h(t, ω) is stationary and bounded by H > 0 and possesses continuous paths and it is assumed that β^{2} > γ^{2} is fulfilled. For κ = 3 the above ODE describes the so-called randomly driven Duffing oscillator.

The second-order ODE can be transformed into a system of first-order ODEs of form (1)

ż = Az + ηd(z) + Bh

with n = 2, r = 1, z = (z_{1}, z_{2})^{τ} = (ẏ, y)^{τ} and

A = ( \begin{matrix}−2γ & −β^{2} \\ 1 & 0 \end{matrix} ) , B = ( \begin{matrix}1 \\ 0 \end{matrix} ) , d(z) = ( \begin{matrix}−z_{2}^{κ} \\ 0 \end{matrix} ) .

Now, applying Theorem 1 the estimate of the radius of convergence of the perturbation series ∑_{p=0}^{∞}^{p}ζ η^{p} with coefficients^{p}ζ given in (5) can be determined. To this end, the eigenvalue decomposition A = VΛV^{−1} with

Λ = diag(λ_{1}, λ_{2}), λ_{1/2} = −γ ± i\sqrt{β^{2} − γ^{2}}
V = ( \begin{matrix}λ_{1} & λ_{2} \\ 1 & 1 \end{matrix} ) , V^{−1} = \frac{1}{λ_{1} − λ_{2}} ( \begin{matrix}1 & −λ_{2} \\ −1 & λ_{1} \end{matrix} )

is used to find λ̂ and V̂. It holds

λ̂ = min_{1≤i≤n} {|Re [λ_{i}]|} = γ

and with

|V||V^{−1}| = \frac{1}{\sqrt{β^{2} − γ^{2}}} ( \begin{matrix}β & β^{2} \\ 1 & β \end{matrix} )

it follows

V̂ = max_{1≤i≤n} \{∑_{j=1}^{n}[|V||V^{−1}|]_{ij}\} = \frac{1 + β}{\sqrt{β^{2} − γ^{2}}} max\{β, 1\}.

Further, Ĥ is found to be

Ĥ = max_{1≤i≤n} {∑_{j=1}^{r} |B_{ij}|} H = max{1, 0}H = H

Using

∑_{i_{1},...,i_{k}=1}^{n} |d_{ki,i_{1},...,i_{k}}| = \{ \begin{matrix}0 & for 2 ≤ k < κ \\ 1 & for k = κ, i = 1 \\ 0 & for k = κ, i = 2 \end{matrix}

it yields

D = max_{1≤i≤n} \{V̂ ∑_{k=2}^{κ} (\frac{V̂ Ĥ}{λ̂} )^{k−1} ∑_{i_{1},...,i_{k}=1}^{n} |d_{ki,i_{1},...,i_{k}}|\}
= max \{V̂ (\frac{V̂ Ĥ}{λ̂} )^{κ−1} , 0\} = V̂^{κ} (\frac{H}{γ} )^{κ−1} .

Then, the estimate η_{0} of the radius of convergence is found to be

η_{0} = \frac{λ̂}{KD} = \frac{γ}{(\frac{κ}{κ−1})^{κ} (κ − 1)V̂^{κ} (\frac{H}{γ} )^{κ−1}}
= (\frac{γ(κ − 1)}{κV̂ H})^{κ} \frac{H}{κ − 1} = (\frac{(κ − 1)γ\sqrt{β^{2} − γ^{2}}}{κ(1 + β) max\{β, 1\}H} )^{κ} \frac{H}{κ − 1}.

In [14], p. 201 ff the above nonlinear oscillator is considered with the parameters γ = 2, β^{2} = 12, κ = 3, and η = 50. The sufficient condition η < η_{0} for the convergence of perturbation series leads to the following condition to the bound H of the excitation

H < (\frac{1}{η(κ − 1)} (\frac{(κ − 1)γ\sqrt{β^{2} − γ^{2}}}{κ(1 + β) max\{β, 1\}})^{κ})^{\frac{1}{κ−1}} ≈ 0.012.

In [14] weakly correlated simulation processes are used possessing a variance σ^{2} = 1 and hence it is

P(|h(t, ω)| > 3σ) ≤ \frac{1}{9} and H ≈ 3σ = 3.

Although a convergence of perturbation series cannot be proved with the assumed parameters the obtained results in [14] from simulation and perturbation series can be compared well. The causes seem to be that the perturbation series converge for a larger domain than it can be established by the proof above and the utilized first terms of the perturbation series approximate well the solution.

6 Conclusions

Solutions of the nonlinear system of ODE (1) with a random inhomogeneous term Bh have been found in a perturbation series representation. Under the assumption that for η = 0 the system is stable and that h is a bounded random function with continuous paths the radius of convergence of the power series has been estimated. In case of a stationary h the perturbation series are also stationary. An open question is whether the condition of boundedness of h can be replaced by boundedness in the mean, stochastic boundedness or other conditions on the distribution of h.

Acknowledgement

The support by the Deutsche Forschungsgemeinschaft (DFG) is greatfully appreciated.

References

[1] L. Arnold and W. Kliemann. Qualitative theory of stochastic systems. In A.T. Bharucha-Reid, editor, Probabilistic analysis and related topics, Vol. 3, pages 1–79, New York, 1983. Academic Press.

[2] L. Arnold and V. Wihstutz. Stationary solutions of linear systems with additive and multiplicative noise. Stochastics, 7:133–155, 1982.

[3] H. Bunke. Gewöhnliche Differentialgleichungen mit zufälligen Parametern. Akademie-Verlag, Berlin, 1972.

[4] A.Ya. Dorogovtsev. Some remarks on differential equations perturbed by periodic random processes. Select. Transl. Math. Stat. Probab., 5:259–269, 1965.

[5] A.Ya. Dorogovtsev. Periodic solutions of differential equations perturbed by stochastic processes (translation from Ukr. Mat. Zh. 41(12):1642-1648, 1989). Ukr. Math. J., 41(12):1412–1419, 1989.

[6] P. Hilton and J. Pedersen. Catalan numbers, their generalization, and their uses. Math. Intell., 13(2):64–75, 1991.

[7] K. Ito and M Nisio. On stationary solutions of a stochastic differential equation. J. Math. Kyoto Univ, 4:1–75, 1964.

[8] R.Z. Khas’minskij. Stochastic stability of differential equations. Sijthoff & Noordhoff, Alphen aan den Rijn, 1980.

[9] W. Ledermann, E. Lloyd, S. Vajda, and C. Alexander. Handbook of applicable mathematics. Supplement. Wiley, Chichester, 1990.

[10] S. Mehlhose. Eigenwerte zufälliger Matrizen von linearen Integralfunktionalen schwach korrelierter Prozesse. Doctoral Thesis, TU Chemnitz, 1998.

[11] S. Mehlhose, J. vom Scheidt, and R. Wunderlich. Random eigenvalue problems for bending vibrations of beams. ZAMM, Z. Angew. Math. Mech., 14 pages, accepted for publication 1999.

[12] N.J.A. Sloane and S. Plouffe. The encyclopedia of integer sequences. Academic Press, San Diego, 1995.

[13] T.T. Soong and M. Grigoriu. Random vibration of mechanical and structural systems. PTR Prentice Hall, Englewood Cliffs NJ, 1993.

[14] J. vom Scheidt, B. Fellenberg, and U. Wöhrl. Analyse und Simulation stochastischer Schwingungssysteme. Leitfäden der angewandten Mathematik und Mechanik, Band 71. B. G. Teubner, Stuttgart, 1994.

[15] J. vom Scheidt and U. Wöhrl. Nonlinear vibration systems with two parallel random excitations. Journal for Analysis and its Applications, 16(1):217–228, 1997.

[16] J. vom Scheidt and R. Wunderlich. Nonlinear random vibrations of vehicles. Z. Angew. Math. Mech., 76(S5):447–448, 1996.