## 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.

### Keywords

- NWP
- time-spectral
- chaotic
- error analysis
- initial condition

## 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*,

*, and*spatial

*domains can be included in the algebraic equations.*parameter

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

*. 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.*real

### 2.1. Generalized weighted residual method

Consider the following set of differential equations,

where

Here,

where

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

where the Chebyshev weight is * integral minimum*in calculus of variation, which states that if

Continuing, we expand Eq. (5),

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

where

and

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

where

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

For a more exhaustive description of the model, see [8]. Suffice it to say the variables

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)

where

An interesting and useful property of Chebyshev polynomials is

Thus, the exact solution at the end state * 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 * conv*” of the solution is larger than the pre-set minimum error, then the current subinterval is halved and the

*” 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.*conv

### 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

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,

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

It is of interest to provide a quantitative measure of the effect of

where

The standard deviations are supplied with indices to indicate the chosen parameter and variable for analysis. For example,

Parameter dependency example, using a single GWRM computation | |||
---|---|---|---|

GWRM sol. | |||

It is immediately seen that the

### 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 | |||
---|---|---|---|

GWRM sol. | |||

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

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

It would also be of interest to study the effect of allowing the initial conditions for all variables

In the table below, initial condition interval average values, as well as standard deviations, are provided for all the variables

Initial condition dependency example, using a single GWRM computation | |||
---|---|---|---|

GWRM sol. | |||

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.

It is found that the ICPD technique with

Numerically, Figure 4a is based on

## 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

*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.*prediction horizon

## 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

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

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

Algorithm 1. GWRM time Intervals. |
---|

1: |

2: |

3: |

4: |

5: |

6: |

7: |

8: |

9: |

10: |

11: |

12: |

13: |

14: |

15: |

16: |

17: |

18: |

19: |

20: |