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 timespectral method. The generalized weighted residual method (GWRM) solves systems of nonlinear ODEs and PDEs using a spectral representation of time. Not being subject to CFLlike criteria, the GWRM typically employs time intervals two orders of magnitude larger than those of timestepping 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
 timespectral
 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 CourantFriedrichsLewy (CFL) condition needs to be satisfied when employing more accurate, explicit FDM, thus limiting time step length.
Another approach is to use a timespectral 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 timespectral method. In this section, we also show how the parameter dependence of the solution can be found from a
2. Method
The timespectral 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
A key characteristic of Chebyshev polynomials, and a main reason for choosing them as basis functions, is their
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
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 semiimplicit root solver (SIR) [11] because it shows superior global stability compared to the Newton method (NM) with linesearch 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
where
An interesting and useful property of Chebyshev polynomials is
Thus, the exact solution at the end state
The pseudocode shows the following: the algebraic equations are collected in a vector
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
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 timedependent 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.
Timespectral 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
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
the parameter dependencies can easily be analyzed since GWRM solutions are analytical
the timespectral method is a highorder 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
5. Conclusion
Spectral methods have a long history when applied to the spatial domain in PDE modeling. The timespectral method GWRM, demonstrated here, provides similar accuracy for the temporal domain. Highorder 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 timespectral method can accurately and efficiently compute the physics in all spatial, temporal, and parameter domains.
A recent study, comparing the timespectral method GWRM with commonly used timestepping methods such as RungeKutta and other, highorder 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:

References
 1.
Scheffel J, Lindvall K, Yik HF. A timespectral approach to numerical weather prediction. Computer Physics Communications. 2018; 226 :127135  2.
Morchoisne Y. Resolution of the NavierStokes equations by a spacetime pseudospectral method. La Recherche Aerospatiale. 1979; 5 :293306  3.
Peyret R, Taylor TD. Computational Methods for Fluid Flow. New York: Springer; 1983  4.
TalEzer H. Spectral methods in time for hyperbolic equations. SIAM Journal on Numerical Analysis. 1986; 23 :126  5.
TalEzer H. Spectral methods in time for parabolic problems. SIAM Journal on Numerical Analysis. 1989; 26 :111  6.
BarYoseph P, Moses E, Zrahia U, Yarin A. Spacetime spectral element methods for onedimensional nonlinear advectiondiffusion problems. Journal of Computational Physics. 1995; 119 :62  7.
Scheffel J. Timespectral solution of initialvalue problems. In: Jang CL, editor. Partial Differential Equations: Theory, Analysis and Applications. New York: Nova Science Publishers, Inc.; 2011. pp. 149  8.
Lorenz EN. Irregularity: A fundamental property of the atmosphere. Tellus A: Published online 2016. 1984; 36 :98110  9.
Mason JC, Handscomb DC. Chebyshev Polynomials. London: Chapman and Hall; 2003  10.
Forsyth AR. Calculus of Variations. New York: Dover Publications Inc.; 1960  11.
Scheffel J, Lindvall K. SIR—An efficient solver for systems of equations. SoftwareX. 2018; 7 :5962  12.
Scheffel J. A spectral method in time for initialvalue problems. American Journal of Computational Mathematics AJCM. 2012; 2 :173193  13.
Riva F, Milanese L, Ricci P. Uncertainty propagation by using spectral methods: A practical application to a twodimensional turbulence fluid model. Physics of Plasmas. 2017; 24 :112  14.
Scheffel J, Lindvall K. J. Scheffel and K. Lindvall. Optimizing TimeSpectral Solution of InitialValue Problems. American Journal of Computational Mathematics AJCM. 2018; 8 :726