*3.7. Function Complexity*

As is to be expected, the computational complexity increases with increasing state dimension. It is therefore of interest to develop an expression that bounds the required number of function evaluations as a function of number of states and expansion order. Due to the many different methods of calculating inner products, all with different computational requirements, the number of functional inner products is what will be enumerated.

Let *x* ∈ S*<sup>P</sup>* be a *P*-dimensional state vector consisting of angular variables, and let *q* ∈ N*<sup>P</sup>* be the expansion order of each element in *x*, where N is the set of natural numbers, including zero. The number of inner products required to calculate the chaos coefficients in Equation (20b) for element *xi* is <sup>2</sup>(*qi* + <sup>1</sup>), where {*i* ∈ N : *k* ≤ *P*} and *qi* is the *i*th element of *q*.

Assume that the mean, variance, and covariance are desired for/between each element. The mean does not require any extra inner products, since the mean is simply the zeroth coefficient. The variance from Equation (23) requires an additional (*qi* + 1)<sup>2</sup> inner products for a raw moment, or *q*2*i* inner products for a central moment. Similarly, the covariance from Equation (24) between the *i*th and *j*th elements requires (*qi* + <sup>1</sup>)(*qj* + 1) additional evaluations for a raw moment and *qiqj* for a central moment. Combining these into one expression, the generalized number of inner product evaluations for raw moments with *P* ≥ 2 is

$$2(q\_1+1)+(q\_1+1)^2+\sum\_{i=2}^{P} \left(2(q\_i+1)+(q\_i+1)^2+\sum\_{j=1}^{i-1}(q\_i+1)(q\_j+1)\right),$$

and for central moments is

$$2(q\_1+1) + q\_1^2 + \sum\_{i=2}^{P} \left( 2(q\_i+1) + q\_i^2 + \sum\_{j=1}^{i-1} q\_i q\_j \right) \dots$$

.

It should be noted that this is the absolute maximum number of evaluations that is required for an entirely angular state. In many cases inner products can be precomputed, the use of orthonormal polynomials reduces the coefficient calculation inner products by two, and expansions using real valued polynomials do not require these inner product calculations at all.

#### **4. Numerical Verification and Discussion**

To test the estimation methods outlined in Section 3.5, a system with two angular degrees of freedom is considered. The correlated, nonlinear dynamics governing this body, the initial mean directions *φ*/*θ*, initial USTDs, and constant angular velocities *φ*˙/ ˙ *θ* are given in Table 2.

**Table 2.** Initial conditions and governing equations of the dynamical system in Section 4.


For every simulation, the run time is 4 s with Δ *T* being 0.05 s; this equates to 81 time steps in each simulation. In Figure 4, the joint pdf propagates from the initial conditions (bottom left) to the final state (top right). The initial joint pdf clearly reflects an uncorrelated bivariate wrapped normal distribution. After being transformed by the dynamics, the final joint pdf exhibits not just translation and scaling, but also rotation: indicating a non-zero correlation between the two angles, which is desired.

**Figure 4.** Evolution of the joint pdf from uncorrelated to correlated bivarite wrapped normal pdf.

For a practical application, the mean and standard deviation/variance of each dimension, as well as the covariance between dimensions is desired. When dealing with angles, the mean direction and the standard deviation can be obtained from the first moment, omitting the second moment. Therefore, only the first moment and the covariance will be discussed. Recall the equations for the first moment and covariance are in terms of chaos coefficients and are given generally in Equations (22) and (24). Because two angles are being estimated, the supports of the integrals in Equations (22) and (24) are set as [− *π*, *π*): but it should be noted that the support is not rigidly defined this way, the only requirement is that the support length be 2 *π*.

Rather than exploit the computational efficiency of methods such as quadrature integral approximations on the unit circle [50–52], the integrals are computed as Riemann sums. Therefore, it is necessary to determine an appropriate number of elements that provides an adequate numerical approximation, while remaining computationally feasible.

Figure 5 show the settling effect that decreasing the elemental angle has on the estimation of the covariance. Note that this figure is used to show the sensitivity of the simulation to the integration variable rather than the actual estimation of the covariance, which will be discussed later in this section. Both plots show the relative error of each estimate with respect to a Monte Carlo simulation of the joint system described in Table 2. Clearly, as the number of elements increases, the estimates begin to converge until the difference between *dθ* = 0.01 rad (629 elements) and *dθ* = 0.005 rad (1257 elements) is perceivable only near the beginning of the simulation. Because of this, it can reasonably be assumed that any error in the polynomial chaos estimate with *dθ* = 0.005 is not attributed to numerical estimation of the integrals in Equations (22) and (24). Additionally, these results should also indicate the sensitivity of the estimate to the integration variable. Even though the dynamics used in this work's examples result in a joint pdf that somewhat resembles wrapped normal pdfs, the number of elements used in the integration must still be quite large. The final numerical element that must be covered is the Monte Carlo. For these examples, 5 × 10<sup>7</sup> samples are used in each dimension.

(**b**) **Figure 5.** Effects of integration variable on estimated covariance relative error. Each line represents a different elemental angle ranging from 0.04 to 0.005 rad. (**a**) Angle Component. (**b**) Length Component.

In each of the examples in Section 4.1, the polynomial chaos estimate first evaluates the Rogers-Szeg˝o polynomials at each of the 1257 uniformly distributed angles (*ξ*), solves Equation (20b) for the chaos coefficients, and uses Equations (22) and (24) to estimate the mean and covariance. After this, the 1257 realizations of the state (*ε*(*ξ*)) are propagated forward in time according to the system dynamics. At each time step the system is expanded using polynomial chaos to estimate the statistics.

## *4.1. Simulation Results*

The estimations of the first moment and covariance of the system described by the simulation parameters in Table 2 are shown in Figures 6 and 7. In both cases, the angle and length of the estimate are presented, rather than the exponential form used in the polynomial chaos expansion. This representation much more effectively expresses the mean and concentration components of the estimate, which are directly of interest when examining statistical moments.

 **Figure 6.** Estimation of the first moment of both angles using a fifth order polynomial chaos expansion compared against Monte Carlo simulation. (**a**) Angle Component. (**b**) Length Component.

Beginning with the mean estimate from a fifth order polynomial chaos expansion in Figure 6, the mean direction of the angle *θ* is nearly identical to the Monte Carlo estimate, while the estimate of *φ* begins to drift slightly as the simulation progresses. Of the two angles, this makes the most since; recalling Table 2, only the dynamics governing *φ* are correlated with *θ*, the dynamics governing *θ* are only dependent on *θ*. In comparison, the estimates of the lengths are both much closer to the Monte Carlo result. Looking closely at the end of the simulation, it can be seen that, again, *θ* is practically identical, and there is some small drift in *φ* downwards, indicating that the estimate reflects a smaller concentration. Effectively, the estimation of the mean reflects some inaccuracy; however, this inaccuracy is partly reflected in the larger dispersion of the pdf.

Similarly to the mean, a small drift can be seen in the estimate of the covariance in Figure 7. In both cases the initial estimate is nearly identical to the Monte Carlo result; however, throughout the simulation a small amount of drift becomes noticeable. While this drift is undesirable, the general tracking of the polynomial chaos estimate to the Monte Carlo clearly shows that the correlation between two angles can be approximated using a polynomial chaos expansion.

 **Figure 7.** Estimation of the covariance between the two angles using a fifth order polynomial chaos expansion compared against Monte Carlo simulation. (**a**) Angle Component. (**b**) Length Component.

#### 4.1.1. Unwrapped Standard Deviation and Joint PDF Assumptions

From the discussion of the generating function for the Rogers-Szeg˝o polynomials Equation (25), it is clear that these polynomials are dependent on the USTD. Unfortunately this means that the polynomials are unique to any given problem, and while they can still be computed ahead of time and looked up, it is not as convenient as problems that use polynomials that are fixed (e.g., Hermite polynomials).

Additionally, the inner product in Equation (23), which describes the calculation of the covariance, requires the knowledge of the joint pdf between the two random variables. In practice, there is no reasonable way of obtaining this pdf; and if there is, then the two variables are already so well know, that costly estimation methods are irrelevant.

It is therefore of interest to investigate what effects making assumptions about the USTD and the joint pdf have on the estimates. The basis polynomials are evaluated when solving for the chaos coefficients Equation (20b) and when estimating the statistical moments Equations (22)–(25) at every time step. If no assumption is made about the USTD, then the generating function in Equation (25) or the three step recursion in Equation (5) must be evaluated at every time step as well. In either case, the computational burden can be greatly reduced if the basis polynomials remain fixed, requiring only an initial evaluation. Additionally, if the same USTD is used for both variables, than the simplification from two to one integrals in Equation (25) can be made.

While only used in the estimation of the covariance, a simplification of the joint pdf will also significantly reduce computation and increase the feasibility of the problem. The most drastic of simplifications is to use the initial, uncorrelated joint pdf. Note that the pdf used in the inner product is mean centered at zero (even for Askey chaos schemes); therefore, the validity of the estimation will not be effected by any movement of the mean.

Assuming the USTD to be fixed at 0.1 for both random variables and the joint pdf to be stationary throughout the simulation led to estimates that are within machine precision of the unsimplified results in Figures 6 and 7. This is to be expected when analyzing Askey-chaos schemes (like Hermite-chaos) that are problem invariant. In instances where the USTD of the wrapped normal distribution is

low enough that probabilities at ± *π* are approximately zero, the wrapped normal distribution is effectively a segmen<sup>t</sup> of the unwrapped normal distribution, because the probabilities beyond ± *π* are approximately equal to zero. However, in problems where the USTD increases, the wrapped normal distribution quickly approaches a wrapped uniform distribution, this makes the time-invariant USTD a poor assumption. While a stationary USTD assumption may not hold as well for large variations in USTD, highly correlated, or poorly modeled, dynamical systems, it shows that some assumptions and simplifications can be made to ensure circular polynomial chaos is a practical method of estimation.

#### 4.1.2. Chaos Coefficient Response

The individual chaos coefficients are not always inspected for problems using Askey-chaos simply due to the commonality of Askey-chaos. The adaptation of polynomial chaos to use the Szeg˝o polynomials, and thus expanding from real to complex valued estimates presents a case that warrants inspection of the chaos coefficients. Figure 8 show the time evolution of the first 13 chaos coefficients (including the zeroth coefficient) that describe the random variable *φ*.

**Figure 8.** Time evolution of the first 13 chaos coefficients describing the random variable *φ*. (**a**) Real Coefficient Evolution. (**b**) Imaginary Coefficient Evolution. (**c**) Imaginary and Real Coefficient Evolution. (**d**) Complex Coefficient Evolution.

What becomes immediately apparent is that the coefficients are roughly anti-symmetrically paired until the length of the coefficient begins to approach zero. In this specific case, the eighth coefficient in Figure 8 initiates this trend. This is the first coefficient that does not have an initial estimate coinciding with lower ordered coefficients. All coefficients following this one show very little response to the system. This is to be expected. Recall the calculation of the chaos coefficient includes the product of polynomial and pdf as well as a division by the self inner product of each polynomial order (i.e., Ψ2*k* (*ζ*)*p*(*ζ*)). The polynomial and pdf product have opposing behaviors when approaching ±*π* from 0. Whereas the polynomial oscillation amplitude increases, the tails of the pdf approach probability values of zero. This ensures the growth in the higher order polynomials is mitigated.

For brevity, only the coefficients from the variable *φ* are shown. These have a much more interesting response than *θ* due to the nature of the dynamics. The most notable part of the coefficients from *θ* is that none of the coefficients ever move beyond the complex unit circle, which from Figure 8c, is clearly not the case for *φ*. In fact, the coefficients describing *θ* stay close to the complex unit circle and just move clockwise about it. Similarly, the eighth and higher coefficient lengths begin collapse to zero rad. For this problem (and presumably most problems) almost all of the information is coming from the first two coefficients. Comparing the estimates using two, three, and ten coefficients yields the same results to within machine precision. This is not surprising when considering the inner products (Table 3) that are required to estimate the covariance; each of the inner products are effectively zero when compared with *φ*0, *φ*0)*pwn* and *φ*1, *φ*1)*pwn* . While having to only compute two significant chaos coefficients makes computation easier, it also limits the amount of information that is used in the estimate; however, for simple problems such as this one, two significant coefficients are satisfactory.


**Table 3.** Rogers-Szeg˝o Inner Products.

#### **5. Conclusions and Future Work**

One method of quantifying the uncertainty of a random variable is a polynomial chaos expansion. For variables that exist only on the real line, this type of expansion has been well studied. This work developed the alterations that must be made for a polynomial chaos expansion to be valid for random variables that exist on the unit circle, specifically the complex unit circle (where the *y* coordinate becomes imaginary).

Previous work has shown that polynomial chaos can be used with the Rogers-Szeg˝o polynomials to estimate the raw moments of a random variable with a wrapped normal distribution. A generalized set of expressions for the mean and covariance of multi-dimensional systems for both real and complex systems has been presented that do not make the assumption that each variable has been expanded with the same set of basis polynomials. An example of two angular random variables—one with correlated dynamics, and one without—has been presented. The mean of each random variable as well as the covariance between them is estimated and compared with Monte Carlo estimates. In the case of the uncorrelated random variable, the mean estimates are highly accurate. For the correlated random variable, the estimate is found to slowly diverge from the Monte Carlo result. A similar small divergence is observed in the covariance estimate; however, the general trend is similar enough to indicate the error is not in the formulation of the complex polynomial chaos. Additionally, an approximation to the basis polynomials and time-varying joint probability density function (pdf) is made, without loss of accuracy in the estimate. From the estimates of the mean and covariance, it is clear that the Rogers-Szeg˝o polynomials can be used as an effective basis for angular random

variable estimation. However, for more complex problems, different polynomials should be considered, specifically polynomials with an appropriate number of non-negligible self inner products.

**Author Contributions:** Conceptualization, C.S.; Methodology, C.S.; Formal Analysis, C.S.; Writing—Original Draft Preparation C.S.; Writing—Review & Editing, K.J.D.; Supervision, K.J.D. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Graduate Assistance in Areas of National Need fellowship.

**Conflicts of Interest:** The authors declare no conflict of interest.
