{'fraction_non_alphanumeric': 0.045027136185788214, 'fraction_numerical': 0.016597248516975895, 'mean_word_length': 4.356261619063715, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Model order reduction for random nonlinear dynamical systems and low-dimensional representations for their quantities of interest

We examine nonlinear dynamical systems of ordinary differential equations or differential algebraic equations. In an uncertainty quantification, physical parameters are replaced by random variables. The inner variables as well as a quantity of interest are expanded into series with orthogonal basis functions like the polynomial chaos expansions, for example. On the one hand, the stochastic Galerkin method yields a large coupled dynamical system. On the other hand, a stochastic collocation method, which uses a quadrature rule or a sampling scheme, can be written in the form of a large weakly coupled dynamical system. We apply projection-based methods of nonlinear model order reduction to the large systems. A reduced-order model implies a low-dimensional representation of the quantity of interest. We focus on model order reduction by proper orthogonal decomposition. The error of a best approximation located in a low-dimensional subspace is analysed. We illustrate results of numerical computations for test examples.


Introduction
Mathematical modelling often generates dynamical systems of ordinary differential equations (ODEs) or differential algebraic equations (DAEs) in science and engineering. A quantity of interest (QoI) is defined by the state variables or inner variables. Several sources of uncertainties may affect the included parameters like modelling errors and measurement errors, for example. In uncertainty quantification, a common approach is the substitution of uncertain parameters by random variables, see [37,40]. Numerical methods often apply orthogonal expansions with unknown coefficient functions and predetermined random-dependent basis functions. The stochastic Galerkin method results in a larger coupled dynamical system. A stochastic collocation technique using a quadrature rule or a sampling method can be written as a single large weakly coupled dynamical system, see [30].
A low-dimensional approximation of the random QoI often exists, where only a few basis functions are required for a sufficiently accurate representation. In a certain context, the representation can be interpreted as a sparse approximation. Such low-dimensional representations have been computed by different methods like least angle regression [5], sparse grid quadrature [7], compressed sensing [9] and ℓ 1 -minimisation [17,18]. Stochastic reduced bases were investigated for random linear systems of algebraic equations in [22,34].
Model order reduction (MOR) becomes favourable in the stochastic Galerkin method and the stochastic collocation approach due to the high dimensionality of the systems. In [28], it was shown that an MOR of the Galerkin system implies a low-dimensional approximation of the QoI in the case of linear dynamical systems. In [29], this strategy was carried over to an MOR of the single auxiliary system of the stochastic collocation method in the linear case. In this paper, we consider nonlinear dynamical systems with random parameters, where the QoI still depends linearly on state variables or inner variables. Thus an MOR generates a low-dimensional representation of the random QoI again.
Several efficient MOR methods exist for linear dynamical systems, see [1,3,14,35]. Error bounds are available by Hardy norms of transfer functions in the frequency domain. Yet MOR for nonlinear dynamical systems still represents a challenging task. Typically, projection-based MOR schemes are applied like the proper orthogonal decomposition (POD), see [1,19], or the trajectory piecewise linear approach, see [21,33]. We focus on the POD method, which identifies a low-dimensional representation of the QoI. Furthermore, the representation can be improved by the computation of a best approximation in the low-dimensional subspace. We prove an error bound for this best approximation. In the first place, our aim is the identification of a sufficiently accurate approximation with a dimensionality as small as possible. In the second place, the MOR methods should decrease the computational effort.
The paper is organised as follows. Section 2 incloses the problem setup. Numerical methods yield large dynamical systems formulated in Section 3. The main part is given by Section 4, where we apply MOR and show error bounds. We demonstrate the results of numerical experiments for two illustrative examples in Section 5.

Problem definition
We describe the task consisting in the identification of low-dimensional representations for QoIs from random nonlinear dynamical systems.

Nonlinear dynamical systems
Let a nonlinear descriptor system be given in the form with matrices and functions depending on physical parameters p ∈ Π ⊆ Ê q . The sizes of the matrices are A, E ∈ Ê n×n , B ∈ Ê n×n in , C ∈ Ê nout×n . The system involves a nonlinear function F : Ê n × Π → Ê n . For regular matrices E, the system consists of ODEs with the state variables x : [t 0 , t end ] × Π → Ê n . For singular matrices E, a system of DAEs is given with the inner variables x. We consider initial value problems with a predetermined function x 0 : Π → Ê n . In the case of DAEs, the initial values have to satisfy consistency conditions and typically depend on the physical parameters of the system.
An input u : [t 0 , t end ] → Ê n in is supplied to the system (1). An output y : [t 0 , t end ] × Π → Ê nout is defined by the state variables or inner variables and the matrix C. Without loss of generality, we restrict the analysis to the case of singleinput-single-output (SISO) with n in = n out = 1. On the one hand, the results also apply to dynamical systems with a general nonlinear dependence of the right-hand side on the input. Moreover, the theory is applicable to autonomous dynamical systems. On the other hand, the linear dependence of the output on the state variables or inner variables is essential in this paper.

Stochastic modelling and orthogonal expansions
We adopt a common approach in UQ, see [37,40], for example. Assuming that the parameters p ∈ Π of the system (1) are uncertain, they are replaced by independent random variables p : Ω → Π on some probability space (Ω, A, µ) with event space Ω, sigma-algebra A and probability measure µ. A joint probability density function ρ : Π → Ê is available in the case of traditional probability distributions. For a measurable function f : Π → Ê, the expected value reads as provided that the integral is finite. The Hilbert space The accompanying norm results to Concerning the dynamical system (1), we assume that x 1 (t, ·), . . . , x n (t, ·), y(t, ·) ∈ L 2 (Π, ρ) pointwise for t ∈ [t 0 , t end ].
Now let a complete orthonormal system (Φ i ) i∈N ⊂ L 2 (Π, ρ) be given. It holds that < Φ i , Φ j >= 0 for i = j and < Φ i , Φ j >= 1 for i = j. We assume that the first basis function is always the constant Φ 1 ≡ 1. It follows that the expansions where the coefficient functions v i : [t 0 , t end ] → Ê n and w i : converge in L 2 (Π, ρ) pointwise for t and j = 1, . . . , n.

Low-dimensional orthogonal representations
In practice, the series (5) have to be truncated to finite approximations. A truncation yields for an integer m ≥ 1. We investigate the output of the random dynamical system (1) as QoI. The truncation error results to for each t due to the orthonormality of the basis functions.
Often orthonormal polynomials are chosen as basis functions due to the theory of the polynomial chaos (PC) expansions, see [10,40]. Therein, the multivariate basis polynomials are just the products of the univariate orthonormal polynomials. Each traditional probability distribution implies their family of univariate orthonormal polynomials, see [38]. To obtain an initial set of basis functions, often all polynomials up to a total degree d are included, which results in basis functions, see [40,Eq. (5.24)]. Even if the total degree is moderate (say 3 ≤ d ≤ 5), the number of basis polynomials becomes large for high numbers q of random parameters.
The error of the complete approach is estimated by y(t, ·) −ȳ(t, ·) L 2 (Π,ρ) ≤ y(t, ·) − y (m) (t, ·) L 2 (Π,ρ) + y (m) (t, ·) −ŷ(t, ·) L 2 (Π,ρ) for each t. The upper bound consists of three terms: the truncation error, the error of the numerical method and the additional error of the low-dimensional approximation. We assume that the first and second term are sufficiently small due to choosing a large initial set of basis functions and a sufficiently accurate numerical method. The third term is analysed in Section 4.5.

Numerical Techniques
We apply two well-known classes of numerical methods for the computation of the unknown coefficient functions (6) in the orthogonal expansions (5) now.

Stochastic Galerkin method
The stochastic Galerkin technique (intrusive method) represents a general approach, which can be applied to all types of differential equations including random variables, see [37,40]. Linear or nonlinear dynamical systems were considered in [2,24,25,26], for example. In our problem, the Galerkin method changes the nonlinear dynamical system (1) into the larger coupled system To define the involved matrices, we employ auxiliary arrays: the symmetric matrix S(p) := (Φ i (p)Φ j (p)) ∈ Ê m×m and the column vector s(p) := (Φ i (p)) ∈ Ê m . Now the matrices of the linear parts in (11) read aŝ using the expected value (3) componentwise and Kronecker products. The nonlinear part is given byF : for i = 1, . . . , m.
The expected values (12) represent probabilistic integrals, which often cannot be calculated analytically. An exception are functions F in (1) consisting of polynomials with low degrees. We require approximations by quadrature formulas for the case of general nonlinear functions F in (1). A quadrature rule or a sampling method is determined by its nodes {p (1) , . . . , p (k) } ⊂ Π and its weights for i = 1, . . . , m. Thus the computational effort for one evaluation ofF comprises k evaluations of F , whereas the summation is negligible.
Concerning initial value problems of the system (11), the initial condition (2) of the original system (1) has to be expanded as well. Consistent initial values are required in the case of DAEs.

Stochastic collocation method
Alternatively, we consider a stochastic collocation (non-intrusive method) employing a quadrature rule or a sampling scheme, see [26,39]. This approach is determined uniquely by the associated nodes and weights. The approximations of the coefficient functions belonging to the QoI (9) becomê for i = 1, . . . , m. Thus k separate initial value problems (1),(2) have to be solved for the different nodes.
We write this approach as a single auxiliary system, which exhibits the form (11) again. Due to another meaning of the inner variables, we reformulatê The matrixĈ ∈ Ê m×kn follows from the formula (14). The meaning of the outputs is the same as in the stochastic Galerkin system (11). The nonlinear part becomesF The initial values result from (2), i.e., x(t 0 , p (ℓ) ) = x 0 (p (ℓ) ) for ℓ = 1, . . . , k. This auxiliary system was already derived and applied for linear dynamical systems in [27,30]. A similar single system was constructed in the case of a non-intrusive approach for dynamical systems with both stochastic noise and random parameters in [23].
In the stochastic collocation method, the resulting system (15) is large and weakly coupled, because it consists of the separate original systems (1), which are connected just by the supply of the same input and the definition of the outputs.

Model order reduction
Although the original dynamical system (1) may be small or medium-sized, the dynamical systems (11) and (15)

Projection-based model order reduction
The nonlinear dynamical systems (11) and (15) represent non-parametric formulations. Without loss of generality, we consider the system (11) as full-order model (FOM). MOR yields a reduced-order model (ROM) in form of a smaller nonlinear dynamical system with matricesĀ,Ē ∈ Ê r×r ,B ∈ Ê r ,C ∈ Ê m×r and a functionF : Ê r → Ê r . The state variables or inner variables arev : [t 0 , t end ] → Ê r , which require initial conditionsv(t 0 ) =v 0 . The input of (16) coincides with the input of (11). The aim is to achieve outputs satisfyingw(t) ≈ŵ(t) for all t.
In projection-based MOR, see [1,11], projection matrices T l , T r ∈ Ê mn×r of full rank are used. It follows that . (17) A transient simulation of the dynamical systems typically requires implicit integration schemes. The smaller dimensionality of the system (16) causes a lower computation work in the linear algebra algorithms. However, the definition of the reduced nonlinear function (17) implies that a function evaluation ofF still takes an evaluation of the original functionF . Significant reductions of the computational effort can be achieved only if one replaces the function evaluations ofF by cheaper approximations, which is also called hyper-reduction in this context. Several techniques have been proposed for this purpose like (discrete) empirical interpolation, missing point estimation or piecewise linear approaches. We refer to [6] and the references therein. A hyper-reduction is not within the scope of this paper. Our main aim is the identification of sufficiently accurate low-dimensional representations.
Furthermore, systems of DAEs often require more sophisticated MOR methods in comparison to systems of ODEs, see [4,15]. Indeed the situation becomes more critical for higher-index DAEs (index larger than one).

Proper orthogonal decomposition
The projection-based ROM (16), (17) is determined uniquely by the matrices T l and T r . We apply the POD to identify a projection matrix. A numerical integration scheme yields a solution of an initial value problem for the larger system (11) or (15) in [t 0 , t end ]. In grid points t 0 < t 1 < · · · < t ℓ−1 = t end , we obtain approximations, which are called snapshots in this context, and collect them in a matrix Often it holds that ℓ < mn. Now a singular value decomposition (SVD), see [13,Ch. 2.5] results in with orthogonal matrices U ∈ Ê mn×mn , Q ∈ Ê ℓ×ℓ and a diagonal matrix Σ ∈ Ê ℓ×ℓ with the singular values σ 1 ≥ σ 2 ≥ · · · ≥ σ ℓ ≥ 0. Let u 1 , . . . , u mn be the columns of U. The right-hand projection matrix is defined as for an arbitrary r ≤ ℓ. Thus just the singular vectors associated to the r dominant singular values are included in the MOR scheme. The left-hand projection matrix is chosen by a Galerkin-type approach as T l = T r .
The POD technique requires a solution of an initial value problem of the FOM to generate the snapshots. Thus saving computational effort is achieved only if the ROM can be reused for other numerical simulations. These additional numerical simulations may involve • other input signals u, • other initial valuesv(t 0 ), The third case will be examined for a test example in Section 5.

Low-dimensional subspaces and approximation
The following derivation is applicable in the case of an arbitrary projection-based MOR as discussed in Section 4.1. The output of the ROM (16) generates an approximation of the QoI in the original system (1) bȳ for t ∈ [t 0 , t end ] and p ∈ Π. Since a linear dependence of the QoI y on the state variables or inner variables x is assumed in (1), the same derivation applies as in [28]. The approximation (20) exhibits the formulation with the inner variables of (16) and the output matrixC = (c ij ). Thus new basis functions are defined. We obtain the low-dimensional representation (10) withw j =v j for j = 1, . . . , r. Hence the coefficients are already identified by the MOR scheme. Furthermore, we will analyse the best approximation in the space spanned by {Ψ 1 , . . . , Ψ r }.
The set of functions (21) is linearly independent in most of the cases. However, they are not orthogonal in general. An orthonormal basis can be constructed from the original basis using an SVD of the output matrixC. This SVD also allows for the deletion of a possible linear dependence or a numerical rank deficiency.
The details are explained in [28].

Best approximation
A projection-based MOR (16), (17) identifies a low-dimensional representation (10) including its time-dependent coefficient functions. In addition, we obtain a best approximation within the subspace spanned by the new basis (21), where just information from the right-hand projection matrix T r is used.
The best approximation with w * = (w * 1 , . . . , w * r ) ⊤ ∈ Ê r is defined by the optimisation problem min w(t) with w = (w 1 , . . . , w r ) ⊤ ∈ Ê r pointwise for t ≥ t 0 . This best approximation can be computed from the solution of the FOM by a linear least squares problem.
We obtain Using (21), basic calculations yield The degrees of freedom arew ∈ Ê r , whereasŵ ∈ Ê m is constant. A necessary condition for a minimum is a vanishing gradient. Thus we achieve which is equivalent to the normal equation, see [36,Ch. 4.8], associated to the linear least squares problem (24).
The least squares problem (24) implies the formula provided that the output matrix exhibits full rank. Due to (25), the best approximation (22) is continuous in time, because the solution of the dynamical system (11) or (15) is assumed to be smooth. The computation of the best approximation requires to solve the FOM. However, the application of the formula (25) afterwards is cheap, because the transformation matrix is time-invariant. A QRdecomposition of the output matrix can be reused at all time points.

Error analysis
The differenceŷ −ȳ between the QoI (9) from a numerical method and the approximation (10) is exactly the error of the MOR approach. Hence this error depends on the individual choice of the projection-based MOR method.
A more detailed examination is feasible for the best approximation using the POD method of Section 4.2. Again the Galerkin system (11) is considered without loss of generality, since the analysis also applies to the collocation system (15). We obtain the following result on the approximation quality with respect to the state variables or the inner variables.
Lemma 1 Letv ∈ C 1 be the exact solution of an initial value problem of the dynamical system (11). Let T r ∈ Ê mn×r be the projection matrix from the POD approach with ℓ > r snapshots. For each t ∈ [t 0 , t end ], the solution v * (t) ∈ Ê r of the linear least squares problem with the time step size ∆t = max{t j − t j−1 : j = 1, . . . , ℓ − 1} and the singular value σ r+1 from the POD.

Proof:
The intermediate value theorem of differential calculus yields componentwisê Since v * represents an optimum, we obtain for arbitraryv ∈ Ê r . The first term can be bounded by and the estimate (27) for each t ∈ [t j , t j+1 ]. For the second term, we apply the equalityv(t j ) = V e j with the matrix (18) and a canonical unit vector e j ∈ Ê ℓ . Let V r ∈ Ê mn×ℓ be the optimal rank-r approximation by the SVD, see [13,Ch. 2.5].
It holds that where the diagonal matrix Σ r ∈ Ê ℓ×ℓ contains only the r dominant singular values. Furthermore, we obtain V r = T r Σ r Q ⊤ with a modified diagonal matrix Σ r . Thus we choose the vectorv = Σ r Q ⊤ e j . It follows that due to the error estimate for the spectral norm in [13,Thm. 2.5.3].
The solutions of the least squares problems (26) are smooth again, because C 1solutions of the FOM are assumed. However, we do not require this smoothness of the optimum in the following. Now we show an error bound on the QoI.
Theorem 2 Letv ∈ C 1 andŷ be the solution and the QoI of an initial value problem of the dynamical system (11). The best approximation y * with respect to the subspace, which is identified by the POD with ℓ > r snapshots, satisfies the error bound for each t ∈ [t 0 , t end ] with the time step size ∆t = max{t j −t j−1 : j = 1, . . . , ℓ−1} and the singular value σ r+1 from the POD.
The spectral norm of the output matrix is often not critical in the estimate (28). For example, if the single output of (1) is just a state variable or inner variable, then it follows that Ĉ 2 = 1. The error bound (28) becomes small in the case of a fast decay of the singular values and a sufficiently small time step size. On the one hand, the singular values are computed a priori, which yields the first term of the estimate. On the other hand, the second term is not computable, because the involved derivatives are unknown.

Illustrative examples
Now we apply an MOR approach to the stochastic Galerkin system as well as a stochastic collocation formulation for two test examples. All numerical calculations are performed within the software package MATLAB (version R2016b).

Scrapie model
We use a model of the scrapie disease from [8, p. 37]. Reaction kinetics yields an autonomous system of ODEṡ for three molecular concentrations x 1 , x 2 , x 3 . The system (29)  In the stochastic modelling, we replace the parameters by independent uniformly distributed random variables with 10% variation around the nominal values. The PC expansion (5) includes multivariate polynomials, which are the products of univariate Legendre polynomials. Our truncated expansion (7) involves all polynomials up to total degree three, which implies m = 56 basis functions. The stochastic Galerkin method yields a dynamical system (11) with mn = 168 state variables. Since just quadratic nonlinearities are given in the right-hand side of the system (29), the right-hand sides of (11) can be evaluated exactly except for roundoff errors. In a stochastic collocation approach, we consider a tensor-product formula of the Gauss-Legendre quadrature, see [36,Ch. 3.6], with k = 3 5 = 243 nodes. The accompanying dynamical system (15) exhibits kn = 729 state variables.
The trapezoidal rule yields the numerical solutions of all initial value problems in this test example. Variable time step sizes are determined by a local error control with relative tolerance ε rel and absolute tolerance ε abs . In the stochastic Galerkin method, we apply the choice ε rel = 10 −4 , ε abs = 10 −6 for the computation of the snapshots, which results in 149 steps and thus 150 snapshots including the initial values. These snapshots are also used to obtain an approximation of the expected values as well as the standard deviations of the three random concentrations shown in Figure 1. Now the POD method requires the SVD (19). illustrates the resulting singular values, which exhibit typically a rapid decay. The projection matrices T l , T r and the ROM (16) follow from the POD technique for arbitrary reduced dimensions. Since the dimensionality of the FOM (11) is relatively small, we cannot expect saving computational effort by an MOR. Nevertheless, we discuss the MOR to show the feasibility of the approach for the identification of a low-dimensional representation.
Now both the FOM (11) and the ROMs (16) are integrated with accuracy requirements ε rel = 10 −3 , ε abs = 10 −6 . Initial values are transformed viav(0) = T ⊤ rv (0). We obtain the approximation (20) from solving the ROM as well as the best approximation (22) for each concentration as QoI. The L 2 (Π, ρ)-error between FOM and ROM can be determined pointwise in time. We evaluate the approximations at 200 equidistant time points in the total time interval, where interpolation is used in time. Figure 3 depicts the maximum L 2 (Π, ρ)-errors occurring at the time points for different reduced dimensions. We recognise that the MOR error decreases rapidly for increasing dimensions at the beginning and then stagnates. The error becomes dominated by the quality of the snapshots, which depends on the density of the associated time points, and the error of the time integration. In contrast, the error of the best approximation decreases further for increasing dimensions and achieves a much lower magnitude.
We repeat the strategy for the stochastic collocation method now. Each system (29) is integrated separately for the nodes of the quadrature with different step sizes for tolerances ε rel = 10 −4 , ε abs = 10 −6 . Interpolation yields 200 snapshots of the weakly coupled system (15) at equidistant time points. The result of the POD is shown in Figure 2 (right), which demonstrates a fast decay of the singular values again. We impose the accuracy requirements ε rel = 10 −3 , ε abs = 10 −6 on the integration of both the FOM (15) and the ROMs, where the same time  step size has to be used for all systems (29) in the FOM now. Initial values are obtained byx(0) = T ⊤ rx (0). Figure 4 illustrates the maximum L 2 (Π, ρ)-errors for the approximations (20) and (22) of the QoIs in the case of different reduced dimensions. A comparison to Figure 3 shows that the quality of the approximations is slightly worse in this stochastic collocation approach. We remark that the accuracies of the stochastic Galerkin method or the stochastic collocation method with respect to the exact random QoI in (1) are not compared, because this item is out of scope. We determine the error of the MOR separately in the Galerkin technique and the collocation scheme.

Transistor amplifier circuit
We consider the electric circuit of a transistor amplifier depicted in Figure 5. A mathematical modelling yields a dynamical system (1) of DAEs (n = 5) with differential index one, which is given in [16, p. 377]. The system is nonlinear, where the mapping F models the bipolar transistor of the circuit and thus involves exponential functions. Three capacitances and six resistances are included in the matrices E and A, respectively. The matrix B ∈ Ê 5×2 inserts the inputs u = (U in , U op ) ⊤ . The operating voltage U op is constant. We put the value U op into the matrix B and the input becomes u = (U in , U op ) ⊤ with U op ≡ 1, because U op is a parameter in the matrix B now. A single time-varying input is supplied by the voltage source U in , which we select as the sinusoidal input The unknowns of the system consist of the five node voltages. The output is defined as the voltage U out at the fifth node. Thus the matrix C in (1) is just a unit vector. Figure 6 shows the numerical solution of an initial value problem (1), (2) for constant physical parameters from [16].
In this test example, all numerical solutions of initial value problems are computed by the backward differentiation formulas (BDF), see [36,Ch. 7.2]. A local error control with tolerances ε rel and ε abs yields adaptive time step sizes as well as adaptive orders.
For an uncertainty quantification, we choose all capacitances, all resistances and the operating voltage as random parameters (q = 10) with independent uniform distributions varying 1% around the constant parameters from above. Concerning the orthogonal expansion (7), we apply all basis polynomials up to total degree three again and obtain m = 286 terms. The stochastic Galerkin method generates a DAE system (11) with mn = 1430 inner variables. We require a quadrature formula (13) to evaluate the nonlinear right-hand side approximatively. We use a sparse grid quadrature, see [12], which is adapted from the Gauss-Legendre rule with level 3 and k = 2441 nodes. Initial values for the system (11) are determined as in the previous example. Now the numerical solution with tolerances ε rel = 10 −5 , ε abs = 10 −6 yields 203 snapshots including the initial values in the time interval [0, T ]. On the one hand, the snapshots imply an approximation of the coefficient functions in the truncated expansion (7) of the output voltage. Figure 7 illustrates the maximum coefficients occurring in the discrete time points. We observe different orders of magnitudes, which indicates a potential for a sufficiently accurate low-dimensional approximation of the QoI. On the other hand, a POD of the snapshots reveals the singular values within (19) shown by Figure 8 (left). We obtain the associated projection matrices T l , T r for arbitrary reduced dimensions.
We discuss the ROM of dimension r = 25 in the longer time interval [0, 3T ] now. Unfortunately, smaller dimensions resulted in a failure of the time integration. Both the FOM (11) and its ROM (16) are integrated using the tolerances ε rel = 10 −3 , ε abs = 10 −6 . The expected value as well as the standard deviation of the output voltage are illustrated by Figure 9. The difference between FOM and ROM is tiny for the expected value and thus not visible. We show the FOM only. However, the standard deviation from the ROM exhibits an incorrect oscillation at some subintervals in time, whereas most parts are matched well. Furthermore, we observe a correct peak of the standard deviation at the beginning t 0 = 0,  because consistent initial values (2) depend on the physical parameters in the case of DAEs. Now we investigate the difference of the approximation (20) from the ROM and the best approximation (22) to the solution (9) of the FOM. Figure 10 (left) depicts the L 2 (Π, ρ)-errors in time. Of course the best approximation always yields a smaller error. Both errors are relatively small in the time interval of the snapshots. The MOR error amplifies successively at later times. The error of the best approximation also increases after the first cycle, while its magnitude becomes stable at a relatively low level. Now we repeat the procedure for a stochastic collocation technique, where we apply the sparse grid quadrature from above. The auxiliary system (15) becomes a DAE with differential index one and kn = 12205 inner variables. Concerning the time integrations, the same tolerances are used as in the stochastic Galerkin approach. The DAEs (1) are solved separately for the nodes of the quadrature, which yields snapshots in the first cycle [0, T ] by interpolation at 200 equidistant time points. Figure 8 (right) shows the singular values from (19) associated to the POD method.
The ROM is investigated for dimension r = 10 in the time interval [0, 3T ] now. Unfortunately, larger dimensions cause a failure of the time integration. The statistics of the FOM (15) and its ROM are portrayed in Figure 11. Again the error between FOM and ROM is small in the case of the expected value. Yet an incorrect oscillation emerges in the standard deviation from the ROM. Figure 10 (right) illustrates the L 2 (Π, ρ)-errors in time for both the representation (20) by MOR and the best approximation (22). We recognise that the behaviour of the MOR error is roughly the same as in the stochastic Galerkin approach. However, the error of the best approximation is worse in the time interval of the snapshots, whereas this error does not increase significantly at later times. In this test example, the collocation method is more efficient than the Galerkin method, because the reduced dimension is smaller.

Conclusions
Orthogonal expansions were applied to the solution of random nonlinear dynamical systems. An MOR for the stochastic Galerkin system or an auxiliary system of a stochastic collocation method implied low-dimensional approximations of the expansions. Firstly, a transient MOR method yields a low-dimensional representation directly. Secondly, the projection of a solution of the full-order model onto the subspace, which is identified by the MOR, generates a best approximation. Numerical simulations demonstrated that both the Galerkin method and collocation techniques are feasible to determine low-dimensional approximations. However, a transient simulation of a reduced-order model becomes critical in the case of complex nonlinear dynamical systems. In contrast, the best approximation is robust and identifies an accurate low-dimensional approximation, while still the full-order model has to be solved.