Spectral Representation of Time and Physical Parameters in Numerical Weather Prediction

Abstract

Numerical weather prediction (NWP) is a difficult task in chaotic dynamical regimes because of the strong sensitivity to initial conditions and physical parameters. As a result, high numerical accuracy is usually necessary. In this chapter, an accurate and efficient alternative to the traditional time stepping solution methods is presented; the time-spectral method. The generalized weighted residual method (GWRM) solves systems of nonlinear ODEs and PDEs using a spectral representation of time. Not being subject to CFL-like criteria, the GWRM typically employs time intervals two orders of magnitude larger than those of time-stepping methods. As an example, efficient solution of the chaotic Lorenz 1984 equations is demonstrated. The results indicate that the method has strong potential for NWP. Furthermore, employing spectral representations of physical parameters and initial values, families of solutions are obtained in a single computation. Thus, the GWRM is conveniently used for studies of system parameter dependency and initial condition error growth in NWP.

1. Introduction

Numerical methods are routinely used when modeling complex systems such as meteorological and atmospheric systems. These systems can be described by nonlinear ordinary and partial differential equations. Nonlinear systems in particular, as opposed to linear systems, constitute specific numerical challenges in terms of processing and memory requirements. In order to efficiently use the computer’s resources, it is important to choose a suitable numerical method [1].

Typically, when solving initial value ODEs, it is common practice to use finite difference methods (FDM) that discretize the temporal domain into finite steps. Information from the initial state and previously computed steps are used to estimate the solution at a “future” discrete time step. This is an intuitive and powerful method because we can easily visualize the cause and effect of each step. However, the problem with these techniques when solving complex systems is that small temporal steps are required in order for a solution to be found. In order to acquire a stable solution, the Courant-Friedrichs-Lewy (CFL) condition needs to be satisfied when employing more accurate, explicit FDM, thus limiting time step length.

Another approach is to use a time-spectral method, which takes a bird’s eye view of the problem [2, 3, 4, 5, 6, 7]. Instead of focusing on small local steps in time, an approximate solution that could fit the entire temporal domain is postulated. This approximate solution is set in the form of a finite series of basis functions. Thus, the unknowns of the equations are thereby changed from the physical variables to the coefficients of the solution ansatz. The procedure allows for highly accurate and global temporal solutions.

In the following, we will present application of the GWRM [7] to the Lorenz 1984 [8] differential equations. This is a set of three coupled, nonlinear, and chaotic ODEs, providing a basic model of the important process of Hadley tropical atmospheric circulation. In Section 2, Method, we present a short derivation of the GWRM. Section 3 will introduce the Lorenz equations and how they can be solved by the time-spectral method. In this section, we also show how the parameter dependence of the solution can be found from a single GWRM computation. We also show that, similarly, the dependence on initial conditions can be found from one computation only. A discussion is given in Section 4, followed by conclusion in Section 5.

2. Method

The time-spectral method employed here is the generalized weighted residual method (GWRM). The physical equations of this Galerkin method are projected onto a space of weighted orthogonal basis functions resulting in a set of algebraic equations. The residual of the differential equation is necessarily zero for the exact solution, thus the numerical method seeks to solve for the coefficients that minimize the residual in the interval where orthogonality holds.

In this procedure, all dimensions can be included, thus the term “generalized.” This means that the dependence on all temporal, spatial, and parameter domains can be included in the algebraic equations.

A key characteristic of Chebyshev polynomials, and a main reason for choosing them as basis functions, is their minimax property [9]. This property states that a Chebyshev series of order $n$ is the best approximation of the same Chebyshev series of order $n+1$ . The implication is that the highest mode can be neglected without deteriorating the approximate solution. This enables computation of nonlinear products entirely within spectral space. Actually, all GWRM numerical computations are carried out in spectral space. In contrast to, for example, Fourier series, Chebyshev series are also conveniently real. Any type of boundary condition can be employed, such as Dirichlet, Neumann, and periodic boundary conditions. Finally, Chebyshev polynomial series are known to converge rapidly when approximating smooth functions.

2.1. Generalized weighted residual method

Consider the following set of differential equations,

$$\frac{\partial \mathbf{u}}{\partial t}-F\left[\mathbf{u}\right]=s,$$ E1

where $\mathbf{u}\left(t,\mathbf{x};\mathbf{p}\right)$ is the unknown variable with a temporal, spatial, and parameter dependence. $F$ is a linear or nonlinear operator, and $s$ is a source term. The GWRM solution ansatz takes the form, in the case of a singular physical dimension and one parameter,

$$u\left(t,x;p\right)\approx \mathrm{U}\left(\tau ,\xi ;P\right)=\underset{k=0}{\overset{K}{{\sum }^{\prime }}}\underset{l=0}{\overset{L}{{\sum }^{\prime }}}\underset{m=0}{\overset{M}{{\sum }^{\prime }}}{a}_{\mathit{klm}}{T}_k\left(\tau \right)T_l\left(\xi \right)T_m\left(P\right).$$ E2

Here, $T_n\left(x\right)=\text{cos}\left(n\text{arccos}\left(x\right)\right)$ is the Chebyshev polynomial of the first kind, $a_{\mathit{klm}}$ are coefficients, and the prime ( ${}^{\text{'}}$ ) denotes that the first term of each summation should be multiplied by $1/2$ . Since Chebyshev polynomials are orthogonal in the interval $\left[-1,1\right]$ , variable transformations are applied for general intervals,

$$\tau =\frac{t-A_t}{B_t},\xi =\frac{x-A_x}{B_x},P=\frac{p-A_p}{B_p},$$ E3

where $A_z=\left(z_1+z_0\right)/2$ and $B_z=\left(z_1-z_0\right)/2$ , $z$ being the variable $t$ , $x$ , and $p$ with their respective upper and lower bounds. The Picard integral is then applied to (1), after which the residual of the differential equation is formed,

$$R=u\left(t,x;p\right)-\left[u\left(t_0,x;p\right)+{\int }_{t_0}^t\left(F\left[u\right]+s\right)d{t}^{\prime }\right].$$ E4

The residual is multiplied by weighted Chebyshev polynomials and then integrated over the computational domain, whereafter it is set to zero. The equation takes the form

$${\int }_{p_0}^{p_1}{\int }_{x_0}^{x_1}{\int }_{t_0}^{t_1}{\mathit{RT}}_q\left(\tau \right)T_r\left(\xi \right)T_s\left(P\right){\mathrm{w}}_t{\mathrm{w}}_x{\mathrm{w}}_p\mathit{\text{dtdxdp}}=0,$$ E5

where the Chebyshev weight is ${\mathrm{w}}_{\zeta }={\left(1-{\zeta }^2\right)}^{-1/2}$ , $\zeta$ being the transform variables $\tau$ , $\xi$ , and $P$ . Eq. (5) has the form of an integral minimum in calculus of variation, which states that if $R$ is smooth and differentiable, and the basis functions and weights are nonzero, then $R$ will be minimized in the domain [10].

Continuing, we expand Eq. (5),

$${\int }_{p_0}^{p_1}{\int }_{x_0}^{x_1}{\int }_{t_0}^{t_1}\left\{u\left(t,x;p\right)-\left[u\right(t_0,x;p)+{\int }_{t_0}^t\left(F\left[u\right]+s\right)d{t}^{\prime }]\right\}{T}_q\left(\tau \right)T_r\left(\xi \right)T_s\left(P\right){\mathrm{w}}_t{\mathrm{w}}_x{\mathrm{w}}_p\mathit{\text{dtdxdp}}=0,$$ E6

yielding the following relations for each term of (6):

$${\int }_{p_0}^{p_1}{\int }_{x_0}^{x_1}{\int }_{t_0}^{t_1}u\left(t,x;p\right)T_q\left(\tau \right)T_r\left(\xi \right)T_s\left(P\right){\mathrm{w}}_t{\mathrm{w}}_x{\mathrm{w}}_p\mathit{\text{dtdxdp}}=B_t{B}_x{B}_p{\left(\frac{\pi }{2}\right)}^3{a}_{\mathit{qrs}}$$ E7

$${\int }_{p_0}^{p_1}{\int }_{x_0}^{x_1}{\int }_{t_0}^{t_1}u\left(t_0,x;p\right)T_q\left(\tau \right)T_r\left(\xi \right)T_s\left(P\right){\mathrm{w}}_t{\mathrm{w}}_x{\mathrm{w}}_p\mathit{\text{dtdxdp}}=B_t{B}_x{B}_p{\left(\frac{\pi }{2}\right)}^2{\mathit{\pi \delta }}_{q0}{b}_{\mathit{rs}}$$ E8

$${\int }_{p_0}^{p_1}{\int }_{x_0}^{x_1}{\int }_{t_0}^{t_1}\left({\int }_{t_0}^{t}F\left[u\right]d{t^{\prime }}^{\text{'}}\right)T_q\left(\tau \right)T_r\left(\xi \right)T_s\left(P\right){\mathrm{w}}_t{\mathrm{w}}_x{\mathrm{w}}_p\mathit{\text{dtdxdp}}=B_t{B}_x{B}_p{\left(\frac{\pi }{2}\right)}^3{A}_{\mathit{qrs}}$$ E9

$${\int }_{p_0}^{p_1}{\int }_{x_0}^{x_1}{\int }_{t_0}^{t_1}\left({\int }_{t_0}^t\mathit{sd}{t}^{\prime }\right)T_q\left(\tau \right)T_r\left(\xi \right)T_s\left(P\right){\mathrm{w}}_t{\mathrm{w}}_x{\mathrm{w}}_p\mathit{\text{dtdxdp}}=B_t{B}_x{B}_p{\left(\frac{\pi }{2}\right)}^3{S}_{\mathit{qrs}}.$$ E10

where

$${\int }_{t_0}^{t}F\left[u\right]d{t}^{\prime }\approx \underset{k=0}{\overset{K}{{\sum }^{\prime }}}\underset{l=0}{\overset{L}{{\sum }^{\prime }}}\underset{m=0}{\overset{M}{{\sum }^{\prime }}}{A}_{\mathit{klm}}{T}_k\left(\tau \right)T_l\left(\xi \right)T_m\left(P\right),$$ E11

and

$${\int }_{t_0}^t\mathit{sd}{t}^{\prime }\approx \underset{k=0}{\overset{K}{{\sum }^{\prime }}}\underset{l=0}{\overset{L}{{\sum }^{\prime }}}\underset{m=0}{\overset{M}{{\sum }^{\prime }}}{S}_{\mathit{klm}}{T}_k\left(\tau \right)T_l\left(\xi \right)T_m\left(P\right).$$ E12

Combining Eqs. (7a)–(7d), we obtain the final set of algebraic equations

$$a_{\mathit{qrs}}=2{\delta }_{q0}{b}_{\mathit{rs}}+A_{\mathit{qrs}}+S_{\mathit{qrs}}$$ E13

where $a_{\mathit{qrs}}$ are the solution coefficients, ${\delta }_{q0}$ is the Kronecker delta function, $b_{\mathit{rs}}$ represents the initial conditions, $A_{\mathit{qrs}}$ is the spectral representation of the time-integrated linear or nonlinear operator $F$ , and $S_{\mathit{qrs}}$ represents the time-integrated source term. Eq. (10) is here defined for the truncated domain $0\le q\le K$ , whereas strictly the time integration renders $K+1$ terms. For the spatial domain, we have $0\le r\le L_{\mathit{BC}}$ , where $L_{\mathit{BC}}=L-N_{\mathit{BC}}$ , with $N_{\mathit{BC}}$ representing the total number of boundary conditions for this spatial dimension. Boundary condition equations are thus added to (10) for problems including spatial derivatives. Here, we will be solving a set of time-dependent ODEs, thus only initial conditions are required to find a unique solution. Finally, it holds that $0\le s\le M$ .

Eq. (10) can be solved iteratively with a suitable root solver. Here, we use the semi-implicit root solver (SIR) [11] because it shows superior global stability compared to the Newton method (NM) with line-search while still maintaining a fast convergence.

For more details regarding GWRM basics, please see references [7, 12].

3. The GWRM, exemplified by the Lorenz 1984 equations

The Lorenz 1984 model is a set of three nonlinear ordinary differential equations that features chaotic behavior similar to that of meteorological systems [8]. It is consequently one of the simplest models of Hadley circulation due to the substantial amount of approximations postulated in order to arrive at this low order model. Thus, whereas the Lorenz 1984 model is not an accurate NWP model, it allows for a rigorous numerical analysis of simple, yet nonlinear chaotic behavior. This is the reason for employing them in the present work. The equations are

$$\begin{array}{l}\frac{d\mathrm{X}}{\mathit{dt}}=-{\mathrm{Y}}^2-{\mathrm{Z}}^2-\alpha \mathrm{X}+\mathit{\alpha Q},\\ \frac{d\mathrm{Y}}{\mathit{dt}}=\mathrm{XY}-\beta \mathrm{XZ}-\mathrm{Y}+W,\\ \frac{d\mathrm{Z}}{\mathit{dt}}=\beta \mathrm{XY}+\mathrm{XY}-\mathrm{Z}.\end{array}$$ E14

For a more exhaustive description of the model, see [8]. Suffice it to say the variables $\mathrm{X}$ , $\mathrm{Y}$ , and $\mathrm{Z}$ represent certain meteorological systems such as wind currents and large-scale eddies. The coefficients $\alpha$ , $\beta$ , $Q$ , and $W$ are chosen within certain limits to act as damping, coupling, and amplification of the physical processes of the system.

Lorenz, in his 1984 article, posed the question, “What can such a simple model possibly tell us about the real atmosphere?”. His answer was that postulating certain hypothesis about the behavior of the real atmosphere, we can strengthen the reasoning behind the hypothesis qualitatively by the help of these models. However, that is not what we are concerned with here, but rather the model is used as a test to show how different numerical methods handle chaotic behavior. Furthermore, the model is useful to develop new techniques that can analyze such quantities as error growth, parameter uncertainties, and variable perturbations in a more efficient way.

3.1. Time intervals

The GWRM is a global method where a finite number of Chebyshev modes are used to represent the physics in a given temporal domain. Instead of using a large number of Chebyshev modes to obtain a single solution for the entire temporal domain, it is more efficient to split the domain into smaller subintervals, so that lesser Chebyshev modes can be used in each interval.

For controlling and maintaining the same numerical accuracy within each time interval, an automatic time-adaptive algorithm is employed. Time-adaptive techniques are easily implemented with Chebyshev polynomials. Since the absolute value of the Chebyshev coefficients decrease with increasing modes for well-resolved functions, the ratio of the absolute values of the highest modes to those of the lowest modes (see Eq. (12)) gives a direct estimate of how “good” the solution is. Thus, it holds that (for the example that (10) represents an ODE)

$$R_K=\frac{\mid a_K\mid +\mid a_{K-1}\mid }{\mid a_0\mid +\mid a_1\mid }\to 0,K\to \infty$$ E15

where $a_k$ , $k=0..K$ , are the Chebyshev coefficients with $K$ as the highest mode. By requiring that $R_K<\mathrm{ϵ}$ , where $\epsilon$ is a stipulated accuracy, the time subintervals to achieve this accuracy are automatically computed. We have found it efficient to use a balanced algorithm that tries to expand the time interval at about every 10th time interval by a factor 1.5, and that halves the interval whenever $\mathrm{ϵ}$ accuracy could not be achieved.

An interesting and useful property of Chebyshev polynomials is $T_k\left(1\right)=1$ , so that

$$u\left(t_1\right)\approx \underset{k=0}{\overset{K}{{\sum }^{\prime }}}{a}_k{T}_k\left(1\right)=\underset{k=0}{\overset{K}{{\sum }^{\prime }}}{a}_k.$$ E16

Thus, the exact solution at the end state $u\left(t_1\right)$ can be approximated by a sum of the Chebyshev coefficients. Eq. (13) conveniently allows for the previous end conditions, represented as Chebyshev coefficients, to be applied as initial conditions for the next interval. All computations of a series of connected intervals are thus calculated in spectral space. A GWRM time-interval pseudocode is provided in Appendix A. Only one domain, the temporal, is presented for the sake of clarity.

The pseudocode shows the following: the algebraic equations are collected in a vector $\varphi$ , along with the initial conditions in the first element of the vector. The time domain is divided into subintervals in which the $\varphi$ equations are solved by SIR. The initial condition is equal to the end-state of the previous solution, see (13). If the convergence parameter “conv” of the solution is larger than the pre-set minimum error, then the current subinterval is halved and the $\varphi$ equations are solved by SIR. On the other hand, if the “conv” convergence criterion is satisfied, the algorithm immediately proceeds to solve the next time interval, the length of which is increased by a factor 1.5 typically every 10th interval. This procedure is done until the entire temporal domain is solved for. When all coefficients are solved for, the semianalytical solution ansatz is easily computed.

3.2. Parameter dependency

Instruments are continually measuring the weather’s temperature, wind speeds, and topological height and width to name a few parameters. However, the measurement devices cannot, at present, give a global coverage. Thus, meteorological parameters are approximate and interpolated.

The meteorological parameters are represented in the Lorenz model by $\alpha$ , $\beta$ , $Q$ , and $W$ . The uncertainty of the measured parameters will inevitably lower the confidence of the solution, specifically important in weather prediction.

When studying parameter dependencies, the GWRM enables an interesting solution. Instead of assigning the parameter a single value, we can introduce a spectrum of values. Thus, we introduce a parameter dimension. This technique was first presented in [7], where it was applied to the viscosity parameter in the 1D nonlinear Burger equation. Later, it was applied to a two-dimensional magnetohydrodynamic problem in the paper by Riva et al. [13].

For the Lorenz ODEs, where we have temporal and parameter dependencies, the solution ansatz then takes the form,

$$\mathrm{u}\left(t;\alpha \right)\approx \mathrm{U}\left(\tau ;P\right)=\underset{k=0}{\overset{K}{{\sum }^{\prime }}}\underset{m=0}{\overset{M}{{\sum }^{\prime }}}{a}_{\mathit{km}}{T}_k\left(\tau \right)T_m\left(P\right).$$ E17

In Eq. (14), a single parameter has been included; however, the procedure can be expanded to any number of parameters. The impact of the parameter $\alpha$ from the Lorenz equation (11) on the $X$ variable is demonstrated in Figure 1. The result from a single GWRM computation is displayed, using the ansatz (14). It can be seen that varying $\alpha$ with the other parameters of Eq. (11) held fixed, the solution is strongly dependent on $\alpha$ for longer times. For lower values of $\alpha$ , the solution becomes more stable and fluctuates with a higher frequency than for higher values of $\alpha$ .

Figure 1.

GWRM solution $X\left(t;\alpha \right)$ showing the parameter $\alpha$ dependency (interval [0.2, 0.3]) on variable $X$ in the time interval $t\in \left[0,5\right]$ . The Lorenz parameters used are $\beta =4.0$ , $Q=8.0$ , and $W=1.0$ . GWRM parameters used are $M=8$ and $K=8$ .

It is of interest to provide a quantitative measure of the effect of $\alpha$ on the solution. One approach would be to compare the deviations from a “base” run for different values of $\alpha$ . For simplicity, however, we have chosen to monitor the standard deviations as a measure of the deflection from the average value in the entire parameter interval. Standard deviations can be computed from the semianalytical solutions with the formula

$$\sigma =\sqrt{\frac{{\sum }_{i=1}^N{\left({\mathrm{U}}_i-\overline{\mathrm{U}}\right)}^2}{N-1}},$$ E18

where ${\mathrm{U}}_i$ is the GWRM solution for a specific value of the parameter. The evaluations have been carried out at $N$ parameter points $i$ at a time of interest. The parameter interval averaged solution is denoted by $\overline{\mathrm{U}}$ . A table is presented below where comparisons are made with a base run with specified parameter values ${\alpha }^{\ast }$ , ${\beta }^{\ast }$ , $Q^{\ast }$ , and $W^{\ast }$ . The solution at $t=1.0$ is used. Average values $\overline{\mathrm{U}}$ and standard deviations $\sigma$ , representing the entire parameter interval are also given.

The standard deviations are supplied with indices to indicate the chosen parameter and variable for analysis. For example, ${\sigma }_{\alpha ,x}$ states that $\alpha$ is the parameter domain and the standard deviation of $X$ is being analyzed. One parameter at a time has been varied by $\pm 20\%$ , while the others are kept constant.

Parameter dependency example, using a single GWRM computation $t=1.0$ , $K=8$ , $M=8$
GWRM sol.	$X\left(1;P^{\ast }\right)=1.311$	$Y\left(1;P^{\ast }\right)=-0.046$	$Z\left(1;P^{\ast }\right)=1.541$
	$\alpha =\left[\mathrm{0.2,0.3}\right]$ , ${\alpha }^{\ast }=0.25$
$\overline{U}\pm \sigma$	$\overline{X}\pm {\sigma }_{\alpha ,x}$	$\overline{Y}\pm {\sigma }_{\alpha ,y}$	$\overline{Z}\pm {\sigma }_{\alpha ,z}$
	$1.307\pm 0.079$	$-0.061\pm 0.353$	$1.498\pm 0.073$
$\beta =\left[\mathrm{3.2,4.8}\right]$ , ${\beta }^{\ast }=4.0$
$\overline{U}\pm \sigma$	$\overline{X}\pm {\sigma }_{\beta ,x}$	$\overline{Y}\pm {\sigma }_{\beta ,y}$	$\overline{Z}\pm {\sigma }_{\beta ,z}$
	$1.322\pm 0.068$	$-0.036\pm 0.584$	$1.401\pm 0.128$
$Q=\left[\mathrm{6.4,9.6}\right]$ , $Q^{\ast }=8.0$
$\overline{U}\pm \sigma$	$\overline{X}\pm {\sigma }_{Q,x}$	$\overline{Y}\pm {\sigma }_{Q,y}$	$\overline{Z}\pm {\sigma }_{Q,z}$
	$1.308\pm 0.096$	$-0.071\pm 0.414$	$1.482\pm 0.087$
$W=\left[\mathrm{0.8,1.2}\right]$ , $W^{\ast }=1.0$
$\overline{U}\pm \sigma$	$\overline{X}\pm {\sigma }_{W,x}$	$\overline{Y}\pm {\sigma }_{W,y}$	$\overline{Z}\pm {\sigma }_{W,z}$
	$1.310\pm 0.010$	$-0.044\pm 0.110$	$1.536\pm 0.015$

It is immediately seen that the $Y$ variable is strongly affected by changes in all parameter values, whereas the $X$ and $Z$ solutions are more robust.

3.3. Initial condition uncertainty

Meteorological models generally include time-dependent PDEs and ODEs that require initial conditions as starting points. In mathematical terms, initial conditions are required to close the system of equations. If we are interested in computing a family of scenarios employing different initial conditions, the standard approach to determine the final states is to restart the computations from scratch, using a new initial condition each time.

Time-spectral methods offer a more convenient approach. Similarly, as when we computed parameter dependency in a single GWRM computation in the previous section, we will now study the effect of a spectrum of initial values on the solution of Eq. (11), employing a single GWRM computation.

Initial condition dependency example, using a single GWRM computation $t=1.0$ , $K=8$ , $M=8$
GWRM sol.	$X\left(1;P^{\ast }\right)=1.311$	$Y\left(1;P^{\ast }\right)=-0.046$	$Z\left(1;P^{\ast }\right)=1.541$
$X\left(0;P\right)\in \left[\mathrm{0.768,1.152}\right]$ , $X^{\ast }=0.96$
$\overline{U}\pm \sigma$	$\overline{X}\pm {\sigma }_{x,x}$	$\overline{Y}\pm {\sigma }_{x,y}$	$\overline{Z}\pm {\sigma }_{x,z}$
	$1.306\pm 0.079$	$-0.047\pm 0.205$	$1.521\pm 0.087$
$Y\left(0;P\right)\in \left[-1.32,-0.88\right]$ , $Y^{\ast }=-1.1$
$\overline{U}\pm \sigma$	$\overline{X}\pm {\sigma }_{y,x}$	$\overline{Y}\pm {\sigma }_{y,y}$	$\overline{Z}\pm {\sigma }_{y,z}$
	$1.319\pm 0.156$	$-0.021\pm 0.440$	$1.460\pm 0.092$
$Z\left(0;P\right)\in \left[\mathrm{0.4,0.6}\right]$ , $Z^{\ast }=0.5$
$\overline{U}\pm \sigma$	$\overline{X}\pm {\sigma }_{z,x}$	$\overline{Y}\pm {\sigma }_{z,y}$	$\overline{Z}\pm {\sigma }_{z,z}$
	$1.309\pm 0.008$	$-0.041\pm 0.132$	$1.535\pm 0.008$

An interesting, albeit challenging, feature of meteorological models is the inherent uncertainty of the initial conditions. As a result, how certain can we be of our numerical results? To be able to predict a chaotic system with absolute precision, one would need infinitely accurate measuring devices, computer, and numerical method. Since none of these are granted, error growth analysis is important to gauge how far into the future we can be confident of our predictions.

A classical error growth analysis employed for the Lorenz equations (11) would be the following. A base scenario ${\mathrm{U}}_{k,1}$ (where ${}^{"}{k}^{"}$ denotes the variable) with initial conditions $\mathbf{v}=<X\left(0\right),Y\left(0\right),Z\left(0\right)>=<0.96,-\mathrm{1.10,0.50}>$ with a time window of $t\in \left[0,50\right]$ is solved for. Subsequently, a number of perturbed scenarios $N$ are solved, denoted with superscript $\left({}^{\prime }\right)$ , where the initial conditions are perturbed an amount $\delta$ , taking values $<0.001$ , for example. The error growth is then computed with the formula,

$$E_k\left(t\right)=\frac{1}{N}\sum _{n=2}^N{\left({\mathrm{U}}_{k,n}^{\prime }-{\mathrm{U}}_{k,1}\right)}^2,k=1,2,\text{and}3.$$ E19

Thus, in this traditional analysis, a large number of computations are needed, where the ODEs are solved for different perturbed initial conditions.

We suggest another approach. By use of an ansatz of the type (14), a spectrum of initial conditions is allowed in a single computation. In the table above, results are presented for three cases where Eq. (11) has been solved to time $t=1.0$ . For the three cases, the base initial conditions $X^{\ast }$ , $Y^{\ast }$ , and $Z^{\ast }$ as well as intervals for $X\left(0,P\right)$ , $Y\left(0,P\right)$ , and $Z\left(0,P\right)$ are shown. Also, results for the base initial conditions at $t=1$ are provided. For the solutions, where the initial conditions are allowed to vary in an interval, both averages over the interval are shown as well as corresponding standard deviations. The analysis shows that the $Y$ variable value at $t=1$ is strongly dependent on the initial condition. In Figure 2, the time and initial condition dependence of the $X$ variable is shown in a 3D diagram. Here, the run time has been extended to $t=5$ .

Figure 2.

GWRM solution $X\left(t;P\right)$ of the Lorenz equations in the time interval $t\in \left[0,5\right]$ , where $P$ is the initial condition parameter varying in the interval $\left[\mathrm{0.768,1.152}\right]$ . The Lorenz parameters used are $\alpha =0.25$ , $\beta =4.0$ , $Q=8.0$ , and $W=1.0$ . GWRM parameters used; $M=8$ and $K=8$ .

It would also be of interest to study the effect of allowing the initial conditions for all variables $X^{\ast }$ , $Y^{\ast }$ , and $Z^{\ast }$ to vary simultaneously. In order to accomplish this in a single GWRM run, the natural thing to do would be to extend the ansatz (14) to include three parameters, corresponding to the different initial conditions. Instead, we have chosen a simpler approach. We let $X^{\prime }\left(0;\delta \right)=X^{\ast }+\delta$ , $Y^{\prime }\left(0;\delta \right)=Y^{\ast }+\delta$ , and $Z^{\prime }\left(0;\delta \right)=Z^{\ast }+\delta$ , where the perturbation $\delta$ is a single parameter that is applied to all three variables. An example computation is shown in Figure 3. The $X$ variable here features a quite a different time evolution than in Figure 2, where only the $X$ variable initial condition was perturbed.

Figure 3.

This figure shows the result of simultaneously varying the initial conditions of all the variables $X$ , $Y$ , and $Z$ . The plot is of the GWRM solution $X\left(t;\delta \right)$ in the time interval $t\in \left[0,5\right]$ , where $\delta$ is the perturbation parameter with interval $\left[-\mathrm{0.1,0.1}\right]$ applied to the initial condition $X^{\prime }\left(0;\delta \right)$ . The Lorenz parameters used are $\alpha =0.25$ , $\beta =4.0$ , $Q=8.0$ , and $W=1.0$ . GWRM parameters used are $M=8$ and $K=8$ .

In the table below, initial condition interval average values, as well as standard deviations, are provided for all the variables $X$ , $Y$ , and $Z$ for a run extending to $t=20$ .

Initial condition dependency example, using a single GWRM computation $t=20$ , $K=8$ , $M=8$ , $\delta <0.001$
GWRM sol.	$X\left(20;0\right)=1.779$	$Y\left(20;0\right)=0.483$	$Z\left(20;0\right)=-0.167$
$\overline{U}\pm \sigma$	$\overline{X}\pm {\sigma }_x$	$\overline{Y}\pm {\sigma }_y$	$\overline{Z}\pm {\sigma }_z$
	$1.787\pm 0.063$	$0.460\pm 0.062$	$-0.099\pm 0.226$

Returning to error growth analysis as suggested by Eq. (16), reliable results can be computed if a sufficiently large domain has been spanned by the perturbed initial conditions. We have compared the classical approach of perturbed iterations (PI) with the new initial condition parameter dependency (ICPD) technique in Figure 4a and b. Of particular interest is the potential gain in CPU time, using the ICPD.

Figure 4.

Comparison of error growth for (a) X and (b) Y (16) for the ICPD and PI schemes, applied to the Lorenz equation (11). (Solid) ICPD computation requiring $t_{\mathit{CPU}}=2.7\text{min}$ and (dashed) PI computation requiring $t_{\mathit{CPU}}=12.7\text{min}$ . The Lorenz parameters used are $\alpha =0.25$ , $\beta =4.0$ , $Q=8.0$ , and $W=1.0$ . GWRM parameters used are $K=8$ and $M=4$ .

It is found that the ICPD technique with $K=8$ and $M=4$ achieves the same result of error growth as the PI scheme. The ICPD scheme computed the error growth in $t_{\mathit{CPU}}=2.7\text{min}$ as compared to the PI scheme which required $t_{\mathit{CPU}}=12.7\text{min}$ . Thus, a near fivefold increase in efficiency was obtained for this relatively simple case.

Numerically, Figure 4a is based on $N=250$ runs with the PI scheme, whereas Figure 4b was computed in a single run, where $N=1000$ different initial condition values was used.

4. Discussion

Is it more efficient to solve a system of differential equations using a spectral representation of the parameter domain instead of solving the system multiple times with different parameter values? Some points related to this discussion are that for the ICPD approach, using the GWRM,

• the parameter dependencies can easily be analyzed since GWRM solutions are analytical

• the time-spectral method is a high-order method, leading to high accuracy solutions

• all parameters in an interval are included, whereas certain regions of critical parameter dependence could be missed out with traditional PI methods

• the parameter domain can be split into intervals, that can be solved for separately, potentially spanning a larger parameter space more efficiently

In this chapter, a univariate analysis has been performed; however, a multivariate analysis can also be implemented in the same manor. It should be noted, however, that the more parameters introduced in the analysis, the more Chebyshev modes need to be solved for, which results in larger matrices and memory demands.

It is also of interest to address the accuracy of the present computations. The ICPD with $K=8$ and $M=4$ did not achieve high accuracy of the exact solution at the end time $t=50$ as can be seen in Figure 5. To be more precise, the solution at higher times ( $\sim t>25$ ) does not represent the real dynamics because the introduced error has grown to an extent where the real solution is lost (see Figure 5). For long run times, the effect of uncertain initial conditions is muddled by the effect of numerical inaccuracy. On top of this, the introduction of the parameter dimension itself in the ansatz (14) is a source for inaccuracy, since the GWRM algorithm has to handle a larger set of Chebyshev polynomials in this case. The argument could be made that if a lower order (e.g., fourth order) solution ansatz cannot accurately represent the parameter “physics,” then the solution exhibits no predictive power if the prediction horizon has been exceeded. How could then the error growth be more accurately represented? The parameter Chebyshev coefficient could perhaps be increased, or the parameter domain could be split into intervals. These are interesting questions for future studies of ICPD methods related to NWP modeling.

Figure 5.

Density plot of $X\left(t;\delta \right)$ in the time interval $t\in \left[0,50\right]$ , where $\delta$ is perturbation parameter with interval $\left[-\mathrm{0.001,0.001}\right]$ applied to the initial condition $X^{\prime }\left(0;\delta \right)$ . GWRM parameters used; $K=8$ and $M=4$ .

5. Conclusion

Spectral methods have a long history when applied to the spatial domain in PDE modeling. The time-spectral method GWRM, demonstrated here, provides similar accuracy for the temporal domain. High-order methods, like the GWRM, may achieve resolutions much higher than that of finite difference methods for similar amount of work. Furthermore, when applied to meteorological systems, the time-spectral method can accurately and efficiently compute the physics in all spatial, temporal, and parameter domains.

A recent study, comparing the time-spectral method GWRM with commonly used time-stepping methods such as Runge-Kutta and other, high-order implicit methods, was carried out in [1]. Since then, the GWRM numerical algorithms have been further streamlined; see for example [14]. Moreover, it has been found that the Jacobian of the algebraic system of equations (10) need only be computed once by including the time interval length analytically. Thus, during the computations, the new time interval length can be substituted into the Jacobian before use. This decreases the computation time from $t_{\mathit{CPU}}=17.9\text{min}$ to $t_{\mathit{CPU}}=12.7\text{min}$ , roughly $30\%$ , for $250$ perturbed scenarios with $K=8$ and $t\in \left[0,50\right]$ .

We have, in this chapter, also presented a new approach to error analysis using the initial condition parameter dependency (ICPD) technique, that is, by spectral representation of the parameter domain. In a single computation, the same error growth was reproduced as for the classical case (where a large number of runs were carried out for different initial conditions), with a computational time of $t_{\mathit{CPU}}=2.7\text{min}$ . This amounts to a near fivefold gain in CPU time efficiency.

Finally, it was shown that the ICPD technique is also successfully applicable for determining the effect of perturbed physical variables on the solution. In one GWRM computation only, it was found that the solutions to the Lorenz 1984 equations sensitively depend on the system parameter $\alpha$ .

Algorithm 1. GWRM time Intervals.
1: procedure
2:        $\varphi \leftarrow A_k\#\text{Algebraic equations}$
3:        $\varphi \left[0\right]\leftarrow A_0+2.0\cdot \mathrm{IC}\#\text{Apply initial condition}\mathrm{IC}$
4:        $t_{\mathit{acc}}\leftarrow 0$
5:        $j\leftarrow 1$
6:       if $t_{\mathit{acc}}>\text{Time}$ them
7:             $\mathit{\text{conv}}\leftarrow 1$ .
8:            while $\text{conv}>\mathrm{rel}.\text{error}$ do
9:                  if $j>1$ them
10:                         $\varphi \left[0\right]\leftarrow A_0+2.0\left(x_{j-1}\left[0\right]/2+x_{j-1}\left[1\right]+\dots +x_{j-1}\left[K\right]\right.)00$ .
11:                   $x_j\leftarrow {SIR}\left(\varphi \right.)00$ .
12:                   $\mathit{\text{conv}}\leftarrow \left(|x\left[K\right]|+|x\left[K-1\right]|\right)/\left(|x\left[0\right]|+|x\left[1\right]|\right)$ .
13:                  if $\text{conv}<\mathrm{rel}.\text{error}$ them
14:                         $a_k\leftarrow x_k\#\text{Save solution}$
15:                  else
16:                         $\mathrm{\Delta }t\leftarrow \mathrm{\Delta }t/2$
17:             $t_{\mathit{acc}}\leftarrow t_{\mathit{acc}}+\mathrm{\Delta }t$
18:            if $\text{modp}\left(j,10\right)=0$ them
19:                    $\mathrm{\Delta }t\leftarrow 1.5\mathrm{\Delta }t$
20:             $j\leftarrow j+1$