BACHELOR THESIS IN MATHEMATICS/APPLIED MATHEMATICS
Comparison of numerical methods for solving a system of ordinary differential equations: accuracy, stability and efficiency
by
Kolar Amir Taher
Kandidatarbete i matematik /till¨ampad matematik
DIVISION OF APPLIED MATHEMATICS M ¨ALARDALEN UNIVERISTY
Bachelor thesis in mathematics/applied mathematics Current version:
7th June 2020 Project name:
Comparison of numerical methods for solving a system of ordinary differential equa-tions: accuracy, stability and efficiency
Author :
Kolar Amir Taher Supervisor :
Siyang Wang Reviewer :
Thomas Westerb¨ack Examiner :
Doghonay Arjmand Comprising:
Abstract
In this thesis, we compute approximate solutions to initial value problems of first-order linear ODEs using five explicit Runge-Kutta methods, namely the forward Euler method, Heun’s method, RK4, RK5, and RK8. This thesis aims to compare the accuracy, stability, and efficiency properties of the five explicit Runge-Kutta methods. For accuracy, we carry out a convergence study to verify the convergence rate of the five explicit Runge-Kutta methods for solving a first-order linear ODE. For stability, we analyze the stability of the five explicit Runge-Kutta methods for solving a linear test equation. For efficiency, we carry out an efficiency study to compare the efficiency of the five explicit Runge-Kutta methods for solving a system of first-order linear ODEs, which is the main focus of this thesis. This system of first-order linear ODEs is a semi-discretization of a two-dimensional wave equation.
Acknowledgments
l would like to express my sincere gratitude to my supervisor, Siyang Wang, who has provided me with a deep understanding of MATLAB features. This understanding allowed me to be able to solve initial value problems of first-order linear ODEs in MATLAB. l would also like to thank Siyang Wang and my reviewer Thomas Westerb¨ack for their input and guidance on how to develop this thesis.
Contents
1 Introduction 6 1.1 Background . . . 6 1.2 Literature review . . . 7 1.3 Research question . . . 9 1.4 Outline . . . 9 2 Differential equations 10 2.1 Partial differential equations . . . 102.2 Ordinary differential equations . . . 13
3 Numerical methods for ODEs 15 3.1 Finite difference methods . . . 15
3.1.1 Errors . . . 16
3.1.2 Numerical properties . . . 16
3.2 Stability of explicit Runge-Kutta methods . . . 18
3.2.1 The forward Euler method . . . 19
3.2.2 Heun’s method . . . 20
3.2.3 The fourth-order Runge-Kutta method . . . 21
3.2.4 The fifth-order Runge-Kutta method . . . 22
3.2.5 The eighth-order Runge-Kutta method . . . 23
3.2.6 Comparison of stability regions . . . 27
4 Numerical experiments 29 4.1 Convergence study . . . 29 4.2 Efficiency study . . . 35 4.2.1 Stability limits . . . 36 4.2.2 Time-error efficiency . . . 37 5 Conclusion 42 Bibliography 46 A MATLAB code 46 A.1 Script 1 . . . 46 A.2 Script 2 . . . 46 A.3 Script 3 . . . 46 A.4 Script 4 . . . 46 A.5 Script 5 . . . 46 A.6 Script 6 . . . 46
A.7 Script 7 . . . 46 A.8 Script 8 . . . 46 A.9 Script 9 . . . 46 A.10 Script 10 . . . 46 A.11 Script 11 . . . 46 A.12 Script 12 . . . 46 A.13 Script 13 . . . 46 A.14 Script 14 . . . 46 A.15 Script 15 . . . 46 A.16 Script 16 . . . 46 A.17 Script 17 . . . 46 A.18 Script 18 . . . 46 A.19 Script 19 . . . 46
List of Figures
3.1 Stability regions of the five explicit Runge-Kutta methods in a com-plex plane . . . 27 4.1 Global truncation errors in the numerical solution of y0(t) = −2ty(t),
y(0) = 1, for t ∈ [0, T ] using the five explicit Runge-Kutta methods against the step size h on a log-log scale with logarithm of base 10. . 32 4.2 Global truncation errors in the numerical solution of ~y0(t) = A~y(t) +
~b(t), ~y(0) = ~y0, for t ∈ [0, 1] computed at increasing stability n number of time points using the forward Euler method. . . 37 4.3 Time (seconds) required for the forward Euler method to compute
global truncation errors in the numerical solution of the system of ODEs. . . 37 4.4 Global truncation errors in the numerical solution of ~y0(t) = A~y(t) +
~b(t), ~y(0) = ~y0, for t ∈ [0, 1] computed at increasing stability n number of time points using Heun’s method. . . 38 4.5 Time (seconds) required for the Heun’s method to compute global
truncation errors in the numerical solution of the system of ODEs. . . 38 4.6 Global truncation errors in the numerical solution of ~y0(t) = A~y(t) +
~b(t), ~y(0) = ~y0, for t ∈ [0, 1] computed at increasing stability n number of time points using RK4. . . 38 4.7 Time (seconds) required for RK4 to compute global truncation errors
in the numerical solution of the system of ODEs. . . 38 4.8 Global truncation errors in the numerical solution of ~y0(t) = A~y(t) +
~b(t), ~y(0) = ~y0, for t ∈ [0, 1] computed at increasing stability n number of time points using RK5. . . 39 4.9 Time (seconds) required for RK5 to compute global truncation errors
in the numerical solution of the system of ODEs. . . 39 4.10 Global truncation errors in the numerical solution of ~y0(t) = A~y(t) +
~b(t), ~y(0) = ~y0, for t ∈ [0, 1] computed at increasing stability n number of time points using RK8. . . 39 4.11 Time (seconds) required for RK8 to compute global truncation errors
List of Tables
4.1 Global truncation errors and approximations of the convergence rate of the forward Euler method for solving y0(t) = −2ty(t), y(0) = 1, for t ∈ [0, 1]. . . 30 4.2 Global truncation errors and approximations of the convergence rate
of Heun’s method for solving y0(t) = −2ty(t), y(0) = 1, for t ∈ [0, 1]. . 30 4.3 Global truncation errors and approximations of convergence rate of
RK4 for solving y0(t) = −2ty(t), y(0) = 1, for t ∈ [0, 1]. . . 31 4.4 Global truncation errors and approximations of the convergence rate
of RK5 for solving y0(t) = −2ty(t), y(0) = 1, for t ∈ [0, 1]. . . 31 4.5 Global truncation errors and approximations of the convergence rate
of RK8 for solving y0(t) = −2ty(t), y(0) = 1, for t ∈ [0, 1.2]. . . 31 4.6 Stability limit n of the five explicit Runge-Kutta methods for ~y0(t) =
Chapter 1
Introduction
1.1
Background
A differential equation is an equation that relates one or several functions and its derivatives [10]. Differential equations are used to model advanced systems, for in-stance, mechanical and electrical systems from fields such as biology, social sciences, engineering, and economics. Mechanical systems and electrical systems are two ex-amples of physical systems that change over time [18]. Differential equations are used to describe a continuous change of physical systems mathematically. Partial differential equations (PDEs) and ordinary differential equations (ODEs) are two classes of differential equations that are used to model and characterize the beha-vior of physical systems. Parabolic PDEs exist in physical problems such as a heat conduction problem for solid bodies. Heat equations describe how the heat tem-perature distributes in, for instance, a rode as the time changes. Wave equations exist to be used to model physical problems such as physical systems, where wave motions are considered [23].
Furthermore, the above-mentioned examples of systems are complex which ac-cording to Houcque [11] are not possible to solve using analytical methods, therefore, numerical methods are used instead. Numerical methods are important to use to solve those differential equations whose exact solutions are not possible to obtain using analytical methods. Analytical methods such as the method of separation of variables give us exact solutions of simple PDEs and ODEs, whereas numerical methods give us approximations of the exact solution. Heath [10] emphasizes that numerical methods are viewed as a significant tool to use to solve problems for ex-ample in engineering and industry. Also, Omale et al. [18] describe that numerical methods have been used to solve differential equations of weather and climate fore-casts. Numerical methods are explicit or implicit computed in one step or multiple steps. An explicit method computes the numerical solution at the next time point using the previous numerical solution at the previous time point. While an implicit method evaluates a function using the numerical solution at the next time point which is solved for.
There are various numerical methods for solving PDEs and ODEs [10]. The method of lines is an example of a numerical method used to find numerical solu-tions of hyperbolic and parabolic PDEs by transforming the PDEs to a system of first-order ODEs which approximates the original PDEs. Finite difference meth-ods, finite element methmeth-ods, and finite volume methods are examples of numerical
methods that are used to approximate the partial derivatives in the PDEs [23]. Finite difference methods such as Runge-Kutta methods are essential for finding approximate solutions of initial value problems of first-order ODEs [3]. High-order Runge-Kutta methods are used to solve ODEs because of their high accuracy and efficiency [10].
Moreover, C. Runge and M. W. Kutta developed explicit and implicit Runge-Kutta methods [1]. The developments of low-order to high-order Runge-Runge-Kutta meth-ods started in 1895. In 1895, C. Runge extended the forward Euler method to more elaborate Runge-Kutta methods with higher accuracy. Runge developed the second-order Runge-Kutta method in 1895. K. Heun introduced the third-order Runge-Kutta method in 1900 [6]. Runge-Kutta methods became important in the studies of explicit and implicit methods for solving ODEs through time discretization [7]. In 1901, W. Kutta introduced the fourth-order and the fifth-order Runge-Kutta methods. F. J. Nystr¨om developed Runge-Kutta methods to be used to solve sys-tems of second-order differential equations [27]. In 1957, R. H. Merson proposed the idea to combine Runge-Kutta methods of different orders in Butcher tableau. In 1964, Butcher introduced the sixth-order Runge-Kutta method with seven stages. The seventh-order Runge-Kutta method with nine stages was recognized from 1968. In 1970, the eighth-order Runge-Kutta method with eleven stages was introduced by Curtis [6].
Runge-Kutta methods are active in research [1] and for the past several years, most research papers have relied on the derivation of new Runge-Kutta methods with higher derivatives to obtain more accurate Runge-Kutta methods [30]. In several papers, researchers have investigated, compared lower-order and high-order Runge-Kutta methods in terms of numerical properties, namely accuracy, stability, and efficiency. In Section 1.2, we review such research papers that have been published for the past several years.
1.2
Literature review
Wusu et al. [30] published a paper on deriving a new Runge-Kutta method involving higher derivatives. In their paper, they derived a three-stage multiderivative expli-cit Runge-Kutta method with first-order and second-order derivatives using the Taylor series expansion of a function. They analyzed the multiderivative explicit Runge-Kutta method in terms of stability and consistency. They also solved two first-order ODEs to compare the accuracy of the multiderivative explicit Runge-Kutta method with Heun’s method and Goeken’s method. Heun’s method and Goeken’s method are also three-stage explicit Runge-Kutta methods but they in-volve fewer derivatives than the multiderivative explicit Runge-Kutta method. The authors computed absolute errors in the numerical solution of the two first-order ODEs using the multiderivative explicit Runge-Kutta method, Heun’s method, and Goeken’s method with different step sizes. An absolute error is the magnitude of the difference between the numerical solution and the exact solution of an ODE. Their research has shown that the absolute errors of the Goeken’s method and the mul-tiderivative explicit Runge-Kutta method do not increase as fast as Heun’s method. The multiderivative explicit Runge-Kutta method did perform better than Heun’s method and Goeken’s method. Their research has also shown that the multideriv-ative explicit Runge-Kutta method is stable, it is the most accurate and efficient
three-stage explicit Runge-Kutta method investigated in their paper.
Papadopoulos and Simos [19] have also derived an explicit Runge-Kutta method, namely the fourth-order modified Runge-Kutta Nystr¨om method. The fourth-order modified Runge-Kutta Nystr¨om method consists of four additional variable coef-ficients than the classical fourth-order Runge-Kutta method. In their research, they have used the fourth-order modified Runge-Kutta Nystr¨om method to solve a Schr¨odinger equation, which is a second-order ODE. The authors have compared the accuracy and efficiency of the fourth-order modified Runge-Kutta Nystr¨om method with three additional types of fourth-order Runge-Kutta Nystr¨om methods for solv-ing the Schr¨odinger equation. They verified the efficiency by computing the error at the final time point. Their research has shown that the fourth-order modified Runge-Kutta Nystr¨om method is more accurate and efficient than the other three types of fourth-order Runge-Kutta Nystr¨om methods.
Another research paper on the accuracy and efficiency of Runge-Kutta methods is written by Hasan [9]. In his paper, he has compared the accuracy and efficiency of the following six Runge-Kutta methods for solving an initial value problem of a first-order linear ODE: the forward Euler method, the backward Euler method, Heun’s method, the fourth-order Runge-Kutta method, ODE23, and ODE45. He computed the numerical solution and exact solution of the initial value problem using the six Runge-Kutta methods with different step sizes. He verified the ac-curacy and efficiency of the six Runge-Kutta methods by computing the average absolute errors and computational time. His research has shown that the ODE45 is the most accurate Runge-Kutta method. The ODE45 is more efficient in terms of function evaluations and computational time than the other Runge-Kutta methods investigated in this paper [9].
Moreover, Islam [12] has also written a research paper on the forward Euler method and the fourth-order Runge-Kutta method for finding numerical solutions of two first-order ODEs with four step sizes. He has compared the accuracy of the two explicit Runge-Kutta methods by computing the maximum errors. Another comparison he did was to compare the efficiency of the two explicit Runge-Kutta methods by verifying how fast the numerical solution converges to the exact solution. His research has shown that the fourth-order Runge-Kutta method is more accurate and the numerical solution obtained by this method converges faster in comparison with the forward Euler method. Also, the fourth-order Runge-Kutta method com-puted the numerical solution of the two first-order ODEs more efficiently than the forward Euler method.
Furthermore, according to Ketcheson and Ahmadia [14] it is significant to find stability polynomials of optimal explicit Runge-Kutta methods for solving first-order ODEs. In their research, they have for instance introduced an approach for how to construct stability polynomials of first-order ODEs. The authors have plotted the stability regions of optimal explicit Runge-Kutta methods such as the fourth-order Runge-Kutta method. Another research paper on stability regions is written by S´eka and Assui [24], where they have proved that the evolution of stability regions does not depend on the order of Runge-Kutta methods. They did this proof by analyzing the stability of the first-order, the second-order, the third-order, and the fourth-order Runge-Kutta method with an eighth-order Runge-Kutta method with eleven stages. Afterwards, they compared the stability regions of the first-order, the second-order, the third-order, and the fourth-order Runge-Kutta method with an
eighth-order Runge-Kutta method with eleven stages which they have introduced. Their research has shown that the size of the stability regions of the second-order, the third-order, and the fourth-order Runge-Kutta methods are larger than the stability region of the eighth-order Runge-Kutta method [24].
In the above-reviewed papers, we see that several researchers have focused to a large extent on finding the accuracy, stability, and efficiency properties of low-order explicit Runge-Kutta methods such as the forward Euler method and Heun’s method. The researchers have less focused on investigating the accuracy, stability, and efficiency properties of high-order explicit Runge-Kutta methods, especially the fifth-order and the eighth-order Runge-Kutta method. Indeed, there is less research on stability analysis of high-order Runge-Kutta methods, in particular the fifth-order Runge-Kutta method. Therefore, in this thesis, we derive the stability analysis for the following five explicit Runge-Kutta methods: the forward Euler method, Heun’s method, the fourth-order Kutta method (RK4), fifth-order Runge-Kutta method (RK5) and the eighth-order Runge-Runge-Kutta method (RK8). We also investigate the accuracy and efficiency properties of the five explicit Runge-Kutta methods.
1.3
Research question
In this thesis, we study the accuracy of the five explicit Runge-Kutta methods for solving an initial value problem of a first-order linear ODE by verifying their convergence rates. We derive the stability analysis to investigate which one of the five explicit Runge-Kutta methods has a better stability property for solving an initial value problem of a linear test equation. This thesis mainly focuses on comparing the efficiency properties of the five explicit Runge-Kutta methods to be able to answer the research question: Which one of the five explicit Runge-Kutta methods is the most efficient for solving an initial value problem for a system of first-order linear ODEs? This system is given and already obtained by the method of lines for a two-space dimensional wave equation, where the wave equation is a hyperbolic PDE. This thesis only focuses on explicit methods to solve this system.
1.4
Outline
The outline of this thesis is as follows. In Chapter 2, we present a theoretical background on PDEs and ODEs. Chapter 3 starts with a theoretical background on finite difference methods and it ends with the stability analysis for the five explicit Runge-Kutta methods. In Chapter 4, we present, analyze, and discuss the results of the convergence study and the efficiency study. Lastly, in Chapter 5 we conclude this thesis and mention future work.
Chapter 2
Differential equations
This chapter presents a theoretical background on PDEs and ODEs that are funda-mental in this thesis.
2.1
Partial differential equations
A partial differential equation (PDE) is a differential equation involving an unknown function of partial derivatives with respect to more than one independent variable [10].
PDEs are classified into hyperbolic, parabolic, and elliptic PDEs. Hyperbolic PDEs model physical processes that are time-dependent and do not develop into a steady-state. Parabolic PDEs model physical processes that are time-dependent and do develop to a steady-state. Elliptic PDEs model processes that are time-independent and have already obtained a steady-state [10].
Yang et al. [31] describe that a wave equation is an example of a hyperbolic PDE and it is used to model wave motion problems. An amplitude of a wave measures the height of the wave. A two-space dimensional wave equation for an amplitude function u(x, y, t) in position (x, y) and time t is written as
∂2u(x, y, t) ∂t2 = B ∂2u(x, y, t) ∂x2 + ∂2u(x, y, t) ∂y2 ! (2.1) for x ∈ [x0, xf], y ∈ [y0, yf], t ∈ [0, T ],
where T is a fixed final time point. Furthermore, u(x, y, t) is a solution of Equation (2.1) which we seek for. This solution defines the dependent variable u is a function of the spatial independent variables, x, y, and the independent variable t representing time. This solution describes the waves that move with the constant wave velocity B. In Equation (2.1), x ranges from the initial spatial point x0 to the final spatial point xf. Also, y ranges from the initial spatial point y0 to the final spatial point yf. And t ranges from the initial time point t = 0 to the fixed final time point t = T .
We want to find u(x, y, t) that satisfies the boundary conditions and the initial conditions. The boundary conditions are
u(x0, y, t) = bx0(y, t), u(xf, y, t) = bxf(y, t), (2.2)
and the initial conditions are u(x, y, 0) = i0(x, y), ∂u ∂t t=0(x, y, 0) = i 0 0(x, y), (2.3)
which should be given in order to solve Equation (2.1). In Equation (2.2), b is a constant which denotes boundary, u(x0, y, t) is a solution at the boundary point x0 which is equal to the boundary condition bx0(y, t) in position y and time t.
Additionally, u(xf, y, t) is a solution at the boundary point xf which is equal to the boundary condition bxf(y, t). Also, u(x, y0, t) is a solution at the boundary point y0
which is equal to the boundary condition by0(x, t) in position x and time t. Also,
u(x, yf, t) is a solution at the boundary point yf which is equal to the boundary condition byf(x, t). In Equation (2.3), u(x, y, 0) is a solution at t = 0 which is equal
to the the initial condition i0 of a spatial index i in position x and y. The first-order derivative of the initial condition i0 is i00 in position x and y [31].
Furthermore, Yang et al. [31] describe that a heat equation in two space di-mensions is an example of a parabolic PDE. A two space dimensional heat equation describes the temperature distribution u(x, y, t) which is written as
∂u(x, y, t) ∂t = B ∂2u(x, y, t) ∂x2 + ∂2u(x, y, t) ∂y2 ! (2.4) for x ∈ [x0, xf], y ∈ [y0, yf], t ∈ [0, T ].
In this context, the temperature distribution is the temperature that spreads in, for example, a wire. The temperature distribution is a function of x, y and time t. In Equation (2.4), B is a diffusion constant [31].
Equation (2.4) has the boundary conditions
u(x0, y, t) = bx0(y, t), u(xf, y, t) = bxf(y, t),
u(x, y0, t) = by0(x, t), u(x, yf, t) = byf(x, t),
and the initial condition is
u(x, y, 0) = i0(x, y).
Furthermore, u(x, y, t) is a solution of Equation (2.4) we seek for that satisfies the boundary conditions and the initial condition [31].
According to Yang et al. [31] Poisson’s equation in two space dimensions is an elliptic PDE and it is written as
∂2u(x, y) ∂x2 +
∂2u(x, y)
∂y2 = f (x, y), (2.5)
for x ∈ [x0, xf], y ∈ [y0, yf],
where f (x, y) is a function of x and y. Equation (2.5) has the boundary conditions u(x0, y) = bx0(y), u(xf, y) = bxf(y), (2.6)
u(x, y0) = by0(x), u(x, yf) = byf(x).
We seek for the solution u(x, y) of Equation (2.5) which satisfies the boundary conditions. In Equation (2.6), u(x0, y) is a solution at the boundary point x0 which
is equal to the boundary condition bx0(y) in position y. Also, u(xf, y) is a solution at
the boundary point xf which is equal to the boundary condition bxf(y). Furthermore,
u(x, y0) is a solution at the boundary point y0 in position x. Additionally, u(x, yf) is a solution at the boundary point yf [31].
Furthermore, the method of lines is a method that solves hyperbolic PDEs and parabolic PDEs that are time-dependent PDEs. The method of lines is a procedure that replaces spatial partial derivatives in a PDE with algebraic approximations numerically. Consequently, the spatial derivatives that were written in the explicit form before are not anymore in terms of x and y. Finite difference methods, finite element methods, and finite volume methods are examples of numerical methods that are used to approximate the spatial derivatives. Consequently, we obtain a system of ODEs with one independent variable t. Numerical methods are used to solve the ODEs to compute an approximate solution of the original PDE [23].
Now we consider the following two-dimensional wave equation to explain more how to proceed with the method of lines
utt = uxx+ uyy, (2.7)
The numerical procedure of the method of lines for Equation (2.7) is as follows. In Equation (2.7), we can see that the function u is dependent on t, x, and y. Furthermore, utt, uxx and uyy are also second-order partial derivatives of u with respect to t, x and y, respectively. Equation (2.7) is changing over time and in space. We want to approximate the second-order partial derivatives utt, uxx, and uyy using finite differences. Therefore, we need to put discrete points in time and space. However, in the method of lines, we do spatial discretization which means that we only discretize x and y in space and not in time. The time variable t is still continuous. We discretize the spatial derivatives uxx and uyy. This procedure is called semi-discretization and it gives a system of ODEs. We put discrete points in the space x and y on intervals. We use the discrete points to approximate the second-order partial derivative utt in time t [23].
Furthermore, after the spatial discretization we get the following system of ODEs ~
utt = D~u, (2.8)
where D is the spatial discretization matrix which according to [15] can be con-structed using, for instance, the finite difference methods, finite element methods, or finite volume methods. Equation (2.8) is in the second-order partial derivative with respect to time t which can’t be solved using the numerical methods for solv-ing ODEs with first-order derivatives. Therefore, we have to transform ~utt to a new equation with first-order derivative with respect to t by introducing a new vector ~v as the following
~
v = ~ut. (2.9)
We must take the first-order derivative with respect to t on both sides of Equation (2.9) and we get
~vt= ~utt = D~u.
We write the vectors, ~u and ~v in one vector in first-order derivative with respect to t ~u ~v t = ~v D~u = 0 I D 0 ~u ~ v , (2.10)
where I is an identity 2 × 2 matrix. The initial value problem of a system of first-order linear ODEs solved in Section 4.2 is in matrix-vector form as Equation (2.10) with
A = 0 I
D 0
,
where A is a 2 × 2 matrix. The system of ODEs in Equation (2.8) should be solved using explicit methods and not implicit methods.
2.2
Ordinary differential equations
An ordinary differential equation (ODE) is a differential equation involving an un-known function of derivatives with respect to one independent variable usually t. The order of an ODE is determined by the highest derivative in the ODE. The highest derivative in a first-order linear ODE is one, therefore this type of ODE is of first-order [10].
Consider this first-order linear ODE taken from [28]
y0(t) = 5y(t), (2.11)
where y is a dependent variable with respect to the independent variable t and y(t) is an unknown solution which we seek for. According to Wang [28] Equation (2.11) is simple, therefore we can use analytical methods to solve it analytically. The general solution which is a family of solutions of Equation (2.11) is
y(t) = Ce5t,
where C is a real constant [28]. If we instead want to find an exact solution of Equation (2.11), then we can impose the following initial condition to determine the value of C which according to [28] it is
y(0) = 2.
Wang [28] describes that the initial condition provide us with information that the initial value of y(t) at the initial time point t0 = 0 is 2. We get that y(0) = Ce0 = C = 2 and the exact solution that satisfies the initial condition of Equation (2.11) is
y(t) = 2e5t, (2.12)
where Equation (2.12) is in agreement with [28]. The first-order ODE with an initial condition is called an initial value problem [28]. Wang [28] describes that a general form of an initial value problem is
y0(t) = f (t, y(t)), (2.13)
y(a) = c, for t ∈ [a, b],
where f (t, y(t)) is a function which depend on t and y. From the initial condition in Equation (2.13), we know that the initial value is c. Also, the time interval is
t ∈ [a, b], the initial time point is t0 = a and the final time point is tn= b, where n is the number of time points.
Heath [10] describes that in most applications of differential equations there exists more than one ODE which is transformed into a system of ODEs. Jung [13] introduce a system of k first-order linear ODEs with coefficients written in this form
y10(t) = a11y1(t) + a12y2(t) + · · · + a1kyk(t) + b1(t)), y20(t) = a21y1(t) + a22y2(t) + · · · + a2kyk(t) + b2(t)), ..
.
yk0(t) = ak1y1(t) + ak2y2(t) + · · · + akkyk(t) + bk(t)),
where aij(t) and bj(t) are known functions on the time interval t ∈ [a, b], i = 1, 2, · · · , n and j = 1, 2, · · · , n. The column vector is ~y(t) of the unknown functions y1(t), · · · , yk(t) which are dependent on t, where k = 1, 2, · · · , n [13]. Jung [13] write this system of k first-order linear ODEs with an initial condition in matrix notation form as
~y0(t) = A~y(t) + ~b(t), (2.14)
~y(t) = ~y0, for t ∈ [a, b].
In Equation (2.14), A is a n × n matrix, ~y0(t) is a column vector and ~b(t) is another column vector. The column vectors depend on t and the column vectors have a length n. If ~b(t) has length zero, then it is a zero column vector, and Equation (2.14) becomes homogeneous. In Equation (2.14), we have the initial value vector ~y0. In the efficiency study, we use numerical methods for ODEs to solve an initial value problem of a homogenous system of first-order linear ODEs formulated as Equation (2.14).
Chapter 3
Numerical methods for ODEs
In this chapter, we present a theoretical background on finite difference methods and analyze the stability of the five explicit Runge-Kutta methods.
3.1
Finite difference methods
Finite difference methods solve for instance initial value problems of first-order ODEs on the time interval t ∈ [a, b]. The idea of finite difference methods for an initial value problem of a first-order ODE is to divide the time interval t ∈ [a, b] into n subintervals and approximate the first-order derivative in the initial value problem for every discrete time point tj in t ∈ [a, b]. In time discretization of initial value problems of first-order ODEs, we obtain discrete solutions that are approximate values of the exact solution of the first-order ODE at discrete time points. Time-discretization is a procedure for solving initial value problems of first-order linear ODEs [31].
Now we show how to proceed with time-discretization of Equation (2.13). We want to find an approximate solutions of Equation (2.13) by discretizing the time interval t ∈ [a, b] into n subintervals using n + 1 points [28]
a = t0 < t1 < t2 < · · · < tn = b.
Wang [28] describes that the length of each n subinterval is the step size hj = tj−tj−1, where j = 1, 2, · · · , n. In this procedure, the previous computed numerical solution yj−1 at the previous discrete time point tj−1 is used to compute the next numerical solution yj at the next discrete time point tj, where j = 0, 1, · · · , n. An approximate solution of Equation (2.13) at the discrete time point tj is yj ≈ y(tj). The numerical solution at t0is y0 = y(t0) = c. We can use the numerical solution y0 = c to compute the next numerical solution y1 at the next discrete time point t1. We can use this procedure to get the numerical solution yn at the final time point tn.
One-step methods such as explicit Runge-Kutta methods use the information from the previously computed numerical solution at a discrete time point to compute the next numerical solution at the next discrete time point. Explicit Runge-Kutta methods are finite difference methods for solving initial value problems of first-order ODEs. Explicit Runge-Kutta methods involve s stages which is equal to the number of function evaluations of the function f (t, y(t)) in Equation (2.13) needed to advance the numerical solution in one time step [10].
3.1.1
Errors
Errors occur when one performs numerical computation on a computer. Errors are classified into round-off errors and truncation errors. Round-off errors are created from the representation of numbers that is approximated. We get round-off errors from numbers such as fractional numbers, π, and the square root of numbers which can not be represented exactly in the computer. Truncation errors are created when a numerical method is approximating an exact solution. Truncation errors are the difference between the approximate solution and the exact solution. Truncation errors are classified into local truncation errors and global truncation errors [10].
Heath [10] describes that the local truncation error `j is `j = yj− uj−1(tj),
where j = 1, 2, · · · , n, the numerical solution of an initial value problem of a first-order ODE is yj and uj−1(t) is the solution of the ODE that passes through the point before, which is (tj−1, yj−1).
According to Heath [10] the global truncation error en is en = yn− y(tn),
the difference between the numerical solution yn and the exact solution y(tn) at the final time point tn.
The differences between a local truncation error and a global truncation error are that we compute the local truncation error in the numerical solution of an ODE at a discrete-time point tj. While we compute the global truncation error at the final time point tnand we obtain a total error at the final time point tn [28]. In this thesis, we only compute global truncation errors because in both the convergence study and the efficiency study we want to compute global truncation errors at the final time point tn. Heath [10] states that global truncation errors are essential to study when evaluating the performance of numerical methods.
3.1.2
Numerical properties
Accuracy, stability, and efficiency are three numerical properties that are used to determine the performance of numerical methods [10].
The following definition of accuracy is from [10].
Definition 3.1.1. The accuracy of a numerical method is said to be of order p if `j = O(hp+1j ).
The motivation for this definition, with the order of accuracy one less than the exponent of the step size in the local error, is that if the local error is O(hp+1j ), then the local error per unit step, `j
hj, is O(h
p
j), and it can be shown that under reasonable conditions the global error en is O(hp), where h is the average step size.
Leveque [16] describe that O is a big-oh notation and it is also known as an asymptotic notation which describe a function that is bounded asymptotically by upper bounds. In Definition 3.1.1, j = 1, 2, · · · , n and the global truncation error
O(hp) is independent of h. We expect the following from a numerical method that has order of accuracy p
eh = Chp+ o(hp) when h → 0,
where eh is an error on a grid of M points and C is in this case an error constant. The little-oh notation o describe a function that is not bounded asymptotically by upper bounds [16].
According to Leveque [16] if h is too small, then eh ≈ Chp. A refinement of the grid by a factor of 2, we get
eh2 ≈ C
h 2
p .
Other factors can certainly be used as well to refine the grid, but refining a grid by a factor of 2 is a standard way.
The ratio of the errors eh and eh2 is
eh eh2
≈ 2p,
where p is a positive integer and the errors decrease approximately by a factor of 2p. Hence, p ≈ log2 eh eh2 .
Furthermore, Leveque [16] describes that we can estimate p with two grid spa-cings h1 and h2 as the following
p ≈
logeeh1h2 log(h1
h2)
(3.1)
If the positive integer p > 0, then the global truncation error O(hp) → 0, when h → 0, and the numerical method converges with the convergence rate p. The rate p of convergence is known as the order of accuracy of a numerical method. The convergence rate p shows how fast the numerical method converges to the exact solution when the step size h → 0. We can measure the accuracy of a numerical method by verifying the convergence rate p of a numerical method [28].
According to S¨oderlind [26], it is not possible to obtain accuracy in numerical solutions if a numerical method is not stable. For this reason, a numerical method must be stable too. The following definition of stability is from [10]
Definition 3.1.2. A numerical method is said to be stable if small perturbations do not cause the resulting numerical solutions to diverge away without bound.
Perturbations are created in numerical computation. Small perturbations are due to round-off errors or truncation errors in the initial data. A small perturbation is when a small change in the initial value of an ODE leads to often a small change
in the numerical solutions of the ODE. We are perturbing the initial value by doing small changes to them. A stable numerical method is not sensitive to small perturb-ations, therefore the numerical solutions do not diverge from the exact solution. For a stable numerical method, the errors in the numerical solutions decrease to zero. Small perturbations to the numerical solutions of a stable ODE get diminished as the time increases because the numerical solution curves of the ODE converge to the exact solution of the ODE. Consequently, the numerical solutions will bound the exact solution.
Furthermore, an unstable numerical method is sensitive to small perturbations, where a small change in the initial value does result in a major change in the numerical solution which causes the numerical solutions to diverge from the exact solution. Small perturbations to a numerical solution of an unstable ODE will grow as the time increases because the numerical solution curves of the ODE diverge. If the numerical solutions of an unstable ODE diverge from the exact solution, then the numerical solutions are unbounded the exact solution [10].
Sabrowski [25] describes that convergence and error analysis can be used to investigate the efficiency of a numerical method. He defines efficiency as
Definition 3.1.3. Numerical efficiency means here the combination of short com-putational time and an acceptable error level.
An important aspect of numerical methods is to verify how much computational time is required for numerical methods to compute errors [22]. We can measure the efficiency of numerical methods by computing errors in the numerical solution of ODEs and computational time [25]. Several researchers have improved the efficiency of Runge-Kutta methods for solving ODEs by decreasing the number of function evaluations.
3.2
Stability of explicit Runge-Kutta methods
S¨oderlind [26] emphasizes that the stability property has a crucial role in solving first-order ODEs using numerical methods. We can verify the stability of numerical methods by analyzing how the numerical methods behave on a linear test equation. In this section, we derive the stability analysis for the five explicit Runge-Kutta methods for solving a linear test equation from [28] defined as
y0(t) = −λy(t), (3.2)
y(0) = 1, t ∈ [0, T ], where λ ∈ C.
We choose Equation (3.2) to be considered in the stability analysis for the five explicit Runge-Kutta methods because it is a simple first-order ODE which according to S¨oderlind [26] one choose in the stability analysis to be able to solve it analytically by pen and paper.
The known exact solution of Equation (3.2) is
y(t) = e−λt. (3.3)
The exact solution (3.3) is also an exponential solution of Equation (3.2) which decrease exponentially to zero when t → ∞ and Re(λ) > 0. The notation, Re(λ) > 0
means positive real values of λ on the real coordinate (Re) axis in a complex plane [28]. A complex plane also consists of an imaginary coordinate (Im) axis which visualizes complex numbers z = hλ. A numerical method is stable for different values of hλ, where h is the step size. Different values of hλ are bounded in a region called stability region which we can visualize in a complex plane. An numerical method is stable if |S(z)|< 1, where S(z) is a stability function. The stability function S(z) is a series in power of z which approximate the exponential solution of a linear test equation [20].
Furthermore, in Section 3.2.1 to Section 3.2.5, we analyze each explicit Runge-Kutta method in terms of stability for solving Equation (3.2). We obtain a stability condition for each explicit Runge-Kutta method in the stability analysis. In Section 3.2.6, we plot the stability conditions to visualize the stability regions of the five explicit Runge-Kutta methods in a complex plane. Lastly, we compare the stability regions to determine which explicit Runge-Kutta method has the smallest and the largest stability region.
3.2.1
The forward Euler method
Atkinson et al. [4] describe that the forward Euler method is a first-order Runge-Kutta method and it is defined as
yj = yj−1+ hf (tj−1, yj−1).
The forward Euler method is used to compute the numerical solution yj at the discrete time point tj. The forward Euler method is a one-stage explicit Runge-Kutta method which means that to advance the numerical solution in one time step, then one function evaluation is needed [28].
We can now derive the stability analysis for the forward Euler method for solving Equation (3.2) as follows
yj = yj−1+ h(−λyj−1), yj = (1 − hλ)yj−1,
where yj, yj−1 are numerical solutions and 1 − hλ is a constant [28]. Furthermore, we have yj−1 = (1 − hλ)yj−2, yj−2 = (1 − hλ)yj−3, .. . y1 = (1 − hλ)y0 = (1 − hλ). Consequently, we get the following numerical solution
yj = (1 − hλ)j. (3.4)
According to Wang [28] the numerical solution (3.4) of Equation (3.2) obtained by the forward Euler method is stable because the numerical solution (3.4) converges to zero. The numerical solution (3.4) decrease in time t as j → ∞. To satisfy this, therefore we need to require the following stability condition
1 − hλ < 1. (3.5)
The stability condition (3.5) constraints the value of the step size h. The forward Euler method is stable when it satisfies the stability condition (3.5). The stability condition (3.5) is an inequality less than 1, therefore we can simplify the stability condition by removing the absolute value and we get
−1 < 1 − hλ < 1 → 0 < hλ < 2 → 0 < h < 2 λ.
The forward Euler method is stable for solving Equation (3.2) if it satisfies the stability condition h < λ2 if λ ∈ R. The inequality 0 < hλ satisfies if λ > 0 according to [28].
We can now investigate how the numerical solution of Equation (3.2) changes for different values of hλ. Assume that we have h = 2, h = 0.1 and λ in Equation (3.2) is λ = 5. If h = 0.1 and λ = 5, then hλ = 0.5 and hλ < 2 which means that the forward Euler method is stable. We can say that the forward Euler method is stable when h < λ2 for λ ∈ R. If h = 2 and λ = 5, then hλ =10 and hλ > 2, which means that the forward Euler method is unstable, therefore the numerical solution diverges to infinity.
Furthermore, the forward Euler method has the following stability function S(hλ)
S(hλ) = 1 − hλ. This stability function alternates in sign.
Wang [28] has also obtained the stability condition (3.5) for the forward Euler method for solving the Equation (3.2). The stability condition for the forward Euler method presented in [2] does not equal the stability condition (3.5) because the test equation solved in [2] is not the same as Equation (3.2). The test equation solved in [2] is
y0(t) = λy(t), (3.6)
therefore we do not obtain the same stability condition for the forward Euler method as in [2].
3.2.2
Heun’s method
Witty [29] describes Heun’s method is a second-order Runge-Kutta method and it is defined as yj = yj−1+ 1 2h k1+ k2 , k1 = f (tj−1, yj−1), k2 = f (tj−1+ h, yj−1+ hk1).
Heun’s method is a two-stage explicit Runge-Kutta method which means that to advance the numerical solution in one time step, then two function evaluations are needed. Heun’s method consists of one more function evaluation than the forward Euler method.
We can now derive the stability analysis for Heun’s method for solving Equation (3.2) as the following k1 = (−λ)yj−1, k2 = − λ + 1 2hλ 2 yj−1, yj = 1 − hλ + 1 2h 2 λ2 yj−1. Consequently, we get the numerical solution as follows
yj =
1 − hλ +1 2h
2λ2j. (3.7)
The numerical solution (3.7) decrease in time t as j → ∞. To satisfy this, we require the following stability condition
1 − hλ + 1 2h 2 λ2 < 1 (3.8)
If we simplify the stability condition (3.8) by removing the absolute value we get, −1 < 1 − hλ + 1 2h 2λ2 < 1 → 0 < hλ − 1 2h 2λ2 < 2 → 0 < hλ < 2 → 0 < h < 2 λ, where h < λ2 is the stability condition for Heun’s method when λ ∈ R and h > 0. This stability condition is the same as the stability condition for the forward Euler method. The stability function of Heun’s method is
S(hλ) = 1 − hλ + 1 2(hλ)
2 .
The stability function of the forward Euler method and Heun’s method are related because the first two terms of the stability function of Heun’s method are the same as the terms of the stability function of the forward Euler method. The stability function of Heun’s method involves three terms, whereas the stability function of the forward Euler method involves two terms.
3.2.3
The fourth-order Runge-Kutta method
Roslan [21] describes that the fourth-order Runge-Kutta method (RK4) is also known as the classical Runge-Kutta method defined as
yj = yj−1+ 1 6h k1+ 2k2+ 2k3+ k4 , k1 = f (tj−1, yj−1), k2 = f tj−1+ 1 2h, yj−1+ 1 2hk1 , k3 = f tj−1+ 1 2h, yj−1+ 1 2hk2 , k4 = f (tj−1+ h, yj−1+ hk3).
RK4 is a four-stage explicit Runge-Kutta method which means that to advance the numerical solution in one time step four function evaluations are needed. For RK4
two more function evaluations are needed than for Heun’s method.
We can now derive the stability analysis for RK4 for solving Equation (3.2) and we get k1 = (−λ)yj−1, k2 = − λ + 1 2hλ 2y j−1, k3 = − λ + 1 2hλ 2− 1 4hλ 3y j−1, k4 = − λ + hλ2−1 2h 2λ3 +1 4h 3λ4y j−1, yj = 1 − hλ + 1 2h 2λ2− 1 6h 3λ3+ 1 24h 4λ4y j−1. Consequently, we get the following numerical solution
yj = 1 − hλ +1 2h 2λ2− 1 6h 3λ3+ 1 24h 4λ4j, (3.9)
which decrease in time t as j → ∞. To satisfy this, we require the following stability condition 1 − hλ + 1 2h 2 λ2− 1 6h 3 λ3+ 1 24h 4 λ4 < 1. The stability function of RK4 is
S(hλ) = 1 − hλ + 1 2(hλ) 2−1 6(hλ) 3+ 1 24(hλ) 4.
This stability function involves five terms which means two more terms than the terms of the stability function of Heun’s method. The first three terms of this stability function are also involved in the stability function of Heun’s method.
3.2.4
The fifth-order Runge-Kutta method
Gopal et al. [8] define the fifth-order Runge-Kutta method (RK5) as yj = yj−1+ 1 90h 7k1+ 32k3+ 12k4+ 32k5+ 7k6 , k1 = f (tj−1, yj−1), k2 = f tj−1+1 4h, yj−1+ 1 4hk1 , k3 = f tj−1+ 1 4h, yj−1+ 1 8hk1 + 1 8hk2 , k4 = f tj−1+ 1 2h, yj−1− 1 2hk2+ hk3 , k5 = f tj−1+ 3 4h, yj−1+ 3 16hk1+ 9 16hk4 , k6 = f tj−1+ h, yj−1− 3 7hk1 + 2 7hk2+ 12 7 hk3− 12 7 hk4+ 8 7hk5 .
RK5 is a six-stage explicit Runge-Kutta method which means six function evalu-ations are needed to advance the numerical solution in one time step. We see that RK5 consists of two more function evaluations than RK4.
We can now derive the stability analysis for RK5 for solving the Equation (3.2) and we get k1 = (−λ)yj−1 k2 = − λ + 1 4hλ 2y j−1 k3 = − λ + 1 4hλ 2− 1 32h 2λ3y j−1 k4 = − λ + 1 2hλ 2−1 8h 2 λ3+ 1 32h 3 λ4yj−1 k5 = (−λ + 12 16hλ 2− 9 32h 2 λ3+ 9 128h 3 λ4− 9 514h 4 λ5yj−1 k6 = − λ + hλ2− 1 2h 2 λ3+ 9 56h 3 λ4 − 3 112h 4 λ5+ 72 3584h 5 λ6yj−1. yj = 1 − hλ + 1 2h 2 λ2−1 6h 3 λ3+ 1 24h 4 λ4− 1 120h 5 λ5+ 1 640h 6 λ6 yj−1. As a consequence, we get the following numerical solution
yj = 1 − hλ +1 2h 2 λ2− 1 6h 3 λ3+ 1 24h 4 λ4 j , (3.10)
which decrease in time t as j → ∞. To satisfy this, we require the following stability condition 1 − hλ + 1 2h 2λ2−1 6h 3λ3+ 1 24h 4λ4− 1 120h 5λ5+ 1 640h 6λ6 < 1. The stability function of RK5 is
S(hλ) = 1 − hλ + 1 2h 2λ2−1 6h 3λ3 + 1 24h 4λ4− 1 120h 5λ5+ 1 640h 6λ6.
The stability function of RK5 involves seven terms which means three more terms than the terms in the stability function of RK4. The first three terms in this stability function are also involved in the stability function of RK4 and Heun’s method but not in the stability function of the forward Euler method.
3.2.5
The eighth-order Runge-Kutta method
The eighth-order Runge-Kutta method (RK8) is defined in [5] with Butcher tableau. The RK8 is yj = yj−1+ h 34 105k6+ 9 35k7+ 9 35k8+ 9 280k9+ 9 280k10+ 41 840k12+ 41 840k13 , k1 = f (tj−1, yj−1), k2 = f tj−1+ 2 27h, yj−1+ 2 27hk1 ,
k3 = f tj−1+ 1 9h, yj−1+ 1 36h k1+ 3k2 ! , k4 = f tj−1+ 1 6h, yj−1+ 1 24h k1+ 3k3 ! , k5 = f tj−1+ 5 12h, yj−1+ 1 48h 20k1− 75k3+ 75k4 ! , k6 = f tj−1+1 2h, yn+ 1 20h k1+ 5k4+ 4k5 ! , k7 = f tj−1+ 5 6h, yj−1+ 1 108h − 25k1+ 125k4− 260k5+ 250k6 ! , k8 = f tj−1+ 1 6h, yj−1+ h 31 300k1+ 61 225k5− 2 9k6+ 13 900k7 ! , k9 = f tj−1+ 2 3h, yj−1+ h 2k1− 53 6 k4+ 704 45 k5− 107 9 k6+ 67 90k7+ 3k8 ! , k10 = f tj−1+ 1 3h, yj−1+ h − 91 108k1+ 23 108k4− 976 135k5+ 311 54 k6− 19 60k7 + 17 6 k8− 1 12k9 ! , k11 = f tj−1+ h, yj−1+ h 2383 4100k1− 341 164k4+ 4496 1025k5− 301 82 k6+ 2133 4100k7, + 45 82k8+ 45 164k9+ 18 41k10 ! , k12 = f tj−1, yj−1+ h( 3 205k1− 6 41k6− 3 205k7 − 3 41k8− 3 41k9+ 6 41k10 ! , k13 = f tj−1+ h, yj−1+ h −1777 4100k1 − 341 164k4 + 4496 1025k5− 289 82 k6+ 2193 4100k7 + 51 82k8+ 33 164k9+ 12 41k10+ k12 ! .
RK8 is a thirteen-stage explicit Runge-Kutta method which means that thirteen function evaluations are needed to advance the numerical solution in one time step. RK8 consists of nine more function evaluations than RK5.
We can now derive the stability analysis for RK8 for solving Equation (3.2) and we get
k2 = − λ + 2 27hλ 2 yj−1, k3 = − λ +1 9hλ 2− 1 162h 2 λ3 yj−1, k4 = − λ +1 6hλ 2− 1 72h 2λ3+ 1 296h 3λ4y j−1 k5 = − λ + 5 12hλ 2− 25 288h 2λ3+ 125 10368h 3λ4− 25 20736h 4λ5y j−1 k6 = − λ +1 2hλ 2− 1 8h 2λ3+ 1 48h 3λ4− 1 384h 4λ5+ 5 20736h 5λ6y j−1 k7 = − λ +5 6hλ 2− 25 72h 2λ3+ 125 1296h 3λ4− 625 31104h 4λ5+ 875 279936h 5λ6 − 625 1119744h 6λ7y j−1 k8 = − λ +1 6hλ 2− 1 72h 2λ3+ 1 1296h 3λ4− 1 31104h 4λ5+ 43 1119744h 5λ6 + 85 10077696h 6λ7+ 325 40310784h 7λ8y j−1, k9 = − λ +2 3hλ 2− 2 9h 2λ3+ 4 81h 3λ4− 2 243h 4λ5+ 1655 559872h 5λ6+ 4279 10077696h 6λ7 + 7865 20155392h 7λ8− 325 13436928h 8λ9y j−1, k10 = − λ + 1 3hλ 2 − 1 18h 2λ3+ 1 162h 3λ4− 1 1944h 4λ5− 757 1119744h 5λ6− 439 1679616h 6λ7 − 3331 20155392h 7λ8+ 65 6718464h 8λ9 − 325 161243136h 9λ10y j−1, k11 = − λ + hλ2− 1 2h 2λ3+1 6h 3λ4 − 1 24h 4λ5+ 14767 1700352h 5λ6− 6511 5101056h 6λ7 + 53 186624h 7λ8− 7151 183638016h 8λ9 + 65 27205632h 9λ10− 325 367276032h 10λ11y j−1, k12 = − λ − 19 566784h 6λ7− 17 45909504h 7λ8− 2081 550914048h 8λ9+ 65 183638016h 9λ10 − 325 1101828096h 10λ11y j−1, k13 = − λ + hλ2− 1 2h 2λ3+1 6h 3λ4 − 1 24h 4λ5+ 14767 1700352h 5λ6− 1585 1275264h 6λ7 + 7297 22954752h 7λ8− 599 17216064h 8λ9+ 12809 2203656192h 9λ10+ 65 275457024h 10λ11 − 325 1101828096h 11λ12y j−1, yj = 1 − hλ +1 2h 2λ2− 1 6h 3λ3+ 1 24h 4λ4− 1 120h 5λ5+ 1 720h 6λ6− 1 5040h 7λ7 + 1 40320h 8λ8− 491 209018880h 9λ9+ 1333 5643509760h 10λ10+ 13 501645312h 11λ11 − 65 4514807808h 12λ12y j−1.
Consequently, we get the following numerical solution yj = 1 − hλ +1 2h 2 λ2 −1 6h 3 λ3+ 1 24h 4 λ4− 1 120h 5 λ5+ 1 720h 6 λ6− 1 5040h 7 λ7 (3.11) + 1 40320h 8λ8− 491 209018880h 9λ9+ 1333 5643509760h 10λ10+ 13 501645312h 11λ11 − 65 4514807808h 12λ12j,
which decrease in time t as j → ∞. To satisfy this, we require the following stability condition 1 − hλ + 1 2h 2 λ2− 1 6h 3 λ3+ 1 24h 4 λ4− 1 120h 5 λ5+ 1 720h 6 λ6− 1 5040h 7 λ7 + 1 40320h 8 λ8− 491 209018880h 9 λ9+ 1333 5643509760h 10 λ10+ 13 501645312h 11 λ11 − 65 4514807808h 12λ12 < 1. The stability function of RK8 is S(hλ) =1 − hλ + 1 2(hλ) 2−1 6(hλ) 3+ 1 24(hλ) 4− 1 120(hλ) 5+ 1 720(hλ) 6 − 1 5040(hλ) 7) + 1 40320(hλ) 8− 491 209018880(hλ) 9+ 1333 5643509760(hλ) 10+ 13 501645312(hλ) 11 − 65 4514807808(hλ) 12.
The stability function of RK8 involves six more terms than the terms of the stability function of RK5. The stability function of RK5 does also involve the first six terms of the stability function of RK8. The seventh term in this stability function of RK8 does not equal the last term of the stability function of RK5 because the fraction is different.
RK8 has a higher number of function evaluations than the forward Euler method, Heun’s method, RK4, RK5, and RK8. The forward Euler method has the lowest number of function evaluations. S¨oderlind [26] discusses that high-order explicit Runge-Kutta methods involving a high number of function evaluations will have higher computational effort than the low-order explicit-Runge-Kutta methods.
S´eka and Assui [24] have also done stability analysis for the forward Euler method, Heun’s method, and RK4 considering Equation (3.6). In their stability analysis, they obtained equations that almost equal to Equation (3.4), Equation (3.7), and Equation (3.9). The equations in [24] do not alternate in sign which Equation (3.4), Equation (3.7), and Equation (3.9) do. The reason for this is that the authors have solved Equation (3.6) is a nonnegative test equation which we have not done in the stability analysis for the five explicit Runge-Kutta methods. How-ever, since the terms of Equation (3.4), Equation (3.7), and Equation (3.9) equal the terms of the equations in [24], then we can state that Equation (3.4), Equation (3.7), and Equation (3.9) is correct. Furthermore, the first eight terms of Equation (3.11) equal the stability function of RK8 with eleven stages in [24]. The terms of this stability function in [24] do not alternate in sign as the terms of Equation (3.11) does. The stability analysis for RK5 and RK8 with thirteen stages is difficult to find
-6 -4 -2 0 2 4 6 Re(h ) -6 -4 -2 0 2 4 6 Im(h ) Forward Euler Heun RK4 RK5 RK8
Figure 3.1: Stability regions of the five explicit Runge-Kutta methods in a complex plane
in the literature. We did the stability analysis for RK5 (see Section 3.2.4) and RK8 (see Section 3.2.5) based on the procedure of the stability analysis for the forward Euler method.
3.2.6
Comparison of stability regions
In Section 3.2.1 to Section 3.2.5, we obtained stability conditions for the five explicit Runge-Kutta methods which are plotted in a complex plane to visualize its stability regions in MATLAB (see Script 1 in Appendix A). MATLAB is a software which we use to perform numerical computation, we visualize the data in figures and tables. We run Script 1 in MATLAB and we get Figure 3.1. We see in Figure 3.1, that the stability region of the forward Euler method is a disk of radius 1. The forward Euler method is centered at the point of hλ = 1. S¨oderlind [26] has also analyzed stability regions of the forward Euler method for solving Equation (3.6) and obtained instead that the forward Euler method is centered at hλ = −1 in a complex plane. The stability region of RK4, RK5, and RK8 have a shape as an ear lobe. The five explicit Runge-Kutta methods are stable inside their regions and unstable outside their regions.
In Figure 3.1, we see that the five explicit Runge-Kutta methods are stable on the positive real coordinate Re(hλ) of the complex plane. The stability regions are on the right-hand side of the complex plane because the five explicit Runge-Kutta methods are stable for values of Re(hλ) > 0, which are bounded in its stability
regions. We can see that RK4, RK5, and RK8 have larger stability regions than the forward Euler method and Heun’s method. The reason is that as the order of explicit Runge-Kutta methods gets higher, then the explicit Runge-Kutta methods become stable for more values of Re(hλ) and its stability regions get larger which is in agreement with [26]. This means that RK4, RK5, and RK8 are stable for more values of Re(hλ) than the forward Euler method and Heun’s method. Additionally, according to Atkinson et al. [4] explicit Runge-Kutta methods get larger stability regions when they are restricted on too small step sizes. The step size of h needs to be small for hλ to be in the stability region and to obtain stable numerical methods. The numerical method with the widest stability region has better stability. This means that RK8 is the most stable of all the five explicit Runge-Kutta methods because it has the largest stability region. Also, the forward Euler method with the smallest stability region is less stable than Heun’s method, RK4, RK5, and RK8.
Moreover, Butcher [6] has also plotted the stability regions of the forward Euler method, Heun’s method, and RK4. He obtained a similar figure as Figure 3.1 but without the stability regions of RK5 and RK8. He has also plotted the stability regions of the forward Euler method, Heun’s method, and RK4 in one complex plane as we have done in Figure 3.1. Therefore, the stability regions of the forward Euler method, Heun’s method, and RK4 is equal to the stability regions presented in [6]. In Figure 3.1, all the stability regions of the five explicit Runge-Kutta methods are plotted in one complex plane.
Furthermore, S´eka and Assui [24] present the result of their research paper that RK8 with eleven stages has a smaller stability region than Heun’s method and RK4. This result does not agree with what we see in Figure 3.1, where we see that RK8 has the largest stability region of all the five explicit Runge-Kutta methods. Even though this result does not agree with what we see in Figure 3.1, this result in [24] is correct because they have investigated the stability of RK8 with eleven stages and not thirteen stages as we have done in this thesis.
Chapter 4
Numerical experiments
This chapter presents two numerical experiments that focus on investigating the accuracy and efficiency properties of the five explicit Runge-Kutta methods. The first numerical experiment is a convergence study and the second numerical exper-iment is an efficiency study. MATLAB is also used in the convergence study and the efficiency study to solve initial value problems numerically and to visualize the results in figures. This is done by running the scripts mentioned in the convergence study and the efficiency study in MATLAB, which can be found in Appendix A.
4.1
Convergence study
In this convergence study, we measure the order of accuracy of the five explicit Runge-Kutta methods for solving an initial value problem of a first-order linear ODE by verifying the convergence rate p. In this study, the following initial value problem of a first-order linear ODE is considered
y0(t) = −2ty(t), y(0) = 1, for t ∈ [0, T ]. (4.1) A similar first-order linear ODE as Equation (4.1) can be found in [10]. Equation (4.1) is chosen to be considered in this study because it is a simple initial value prob-lem with a known exact solution which enables us to compute the global truncation error at the final time T .
The exact solution of Equation (4.1) is y(t) = e−t2.
In this study, we verify the convergence rate p through three steps. The first step is to compute a numerical solution yh
nof Equation (4.1) with the step size h. We choose the step sizes h, h2, h4 and h8, which decrease by a factor of 2 because according to Leveque [16] it is a common way to make h small to get highly accurate numerical solutions that we want to obtain. Additionally, other factors can certainly be used as well. The second step is that we want to compute the following global truncation error from [28] at the final time T
eh = |yhn− y(tn)|. (4.2)
The global truncation error eh with step size h is the difference between the the numerical solution yh
to estimate the convergence rate p by computing three approximations of the con-vergence rate p using Equation (3.1).
We solve Equation (4.1) using the forward Euler method, Heun’s method, RK4 and RK5 on the time interval t ∈ [0, T ], where T = 1 and with the step sizes h = 0.1, h
2, h 4 and
h
8. We solve Equation (4.1) using RK8 on t ∈ [0, T ], where T = 1.2 and with the step sizes h = 0.4, h2, h4 and h8. We have chosen other step sizes h for RK8 because RK8 reaches faster to the machine precision error 10−16 which is the smallest relative error we can get by a computer. Heath [10] describes a relative error as a quotient of the global truncation error and the exact solution. Furthermore, RK8 reaches faster to the machine precision error because it has a high accuracy level compared to for instance the forward Euler method. In the convergence study, we do not want to reach this limit of error because if we reach this limit, then we will not be able to decrease the errors even more. Therefore, for RK8 we need to use step sizes h that are larger than h = 0.1, h
2, h 4 and
h
8 just to make sure that the computed errors using RK8 do not decrease too fast. This is why the step sizes h = 0.4, h2, h4 and h8 are used instead. Furthermore, the time interval t ∈ [0, 1.2] was used instead of t ∈ [0, 1] because if we start to solve Equation (4.1) with h = 0.4, then h = 0.4 will not fit in t ∈ [0, 1]. However, it does fit in t ∈ [0, 1.2]. Script 2 is used to define Equation (4.1) which is solved using the forward Euler method (see Script 3), Heun’s method (see Script 4), RK4 (see Script 5), RK5 (see Script 6) and RK8 (see Script 7). The data obtained with Script 2 to Script 7 were put in a tabular format in LATEX to get Table 4.1 to Table 4.5. The results shown in Table 4.1 to Table 4.5 are also presented in Figure 4.1 which is obtained by running Script 8.
h Global truncation
er-rors Approximations of convergence rate 0.1 1.3827×10−2 — 0.05 6.5045×10−3 ≈ 1.0880 0.025 3.1569×10−3 ≈ 1.0430 0.0125 1.5554×10−3 ≈ 1.0210
Table 4.1: Global truncation errors and approximations of the convergence rate of the forward Euler method for solving y0(t) = −2ty(t), y(0) = 1, for t ∈ [0, 1].
h Global truncation
er-ror Approximation of con-vergence rate 0.1 1.1739×10−3 — 0.05 3.0109×10−4 ≈ 1.9630 0.025 7.6014×10−5 ≈ 1.9860 0.0125 1.9085×10−5 ≈ 1.9940
Table 4.2: Global truncation errors and approximations of the convergence rate of Heun’s method for solving y0(t) = −2ty(t), y(0) = 1, for t ∈ [0, 1].
h Global truncation er-rors Approximations of convergence rate 0.1 1.6252×10−6 — 0.05 1.0253×10−7 ≈ 3.9860 0.025 6.4067×10−9 ≈ 4.0000 0.0125 3.9993×10−10 ≈ 4.0020
Table 4.3: Global truncation errors and approximations of convergence rate of RK4 for solving y0(t) = −2ty(t), y(0) = 1, for t ∈ [0, 1].
h Global truncation
er-rors Approximations of convergence rate 0.1 2.8055×10−8 — 0.05 8.3665×10−10 ≈ 5.0670 0.025 2.5427×10−11 ≈ 5.0401 0.0125 7.8292×10−13 ≈ 5.0210
Table 4.4: Global truncation errors and approximations of the convergence rate of RK5 for solving y0(t) = −2ty(t), y(0) = 1, for t ∈ [0, 1].
h Global truncation
er-rors Approximations of convergence rate 0.4 1.1530×10−7 — 0.2 5.2249×10−10 ≈ 7.7860 0.1 1.9577×10−12 ≈ 8.0600 0.05 7.3552×10−15 ≈ 8.0560
Table 4.5: Global truncation errors and approximations of the convergence rate of RK8 for solving y0(t) = −2ty(t), y(0) = 1, for t ∈ [0, 1.2].
10-2 h 10-1 10-15
10-10 10-5 100
global truncation errors
Forward Euler Heun RK4 RK5 RK8
Figure 4.1: Global truncation errors in the numerical solution of y0(t) = −2ty(t), y(0) = 1, for t ∈ [0, T ] using the five explicit Runge-Kutta methods against the step size h on a log-log scale with logarithm of base 10.
Table 4.1 to Table 4.5 shows global truncation errors and three approximations of the convergence rate p of the five explicit Runge-Kutta methods. As the step sizes h in Table 4.1 to Table 4.5 is halved, the global truncation errors of the five explicit Runge-Kutta methods also decrease by a factor of approximately 2p. According to Leveque [16], if global truncation errors decrease by a factor of approximately 2p, then for a forward Euler method which is first-order accurate with p = 1 and we get that the global truncation errors of the forward Euler method decrease by approximately a factor of 2 because 21 = 2. This means that we should expect that the five explicit Runge-Kutta methods decrease by a factor of approximately 2p. We can now check if the global truncation errors in Table 4.1 to Table 4.5 do decrease by a factor of approximately 2p. We check this by investigating how much the first global truncation error in the tables decreases.
In Table 4.1, the first global truncation error decreases by a factor of approxim-ately 2 because
1.3827 × 10−2
6.5045 × 10−3 ≈ 2.1260. (4.3)
From Table 4.1, we can see that as the step sizes decrease by half, the global trun-cation errors in Table 4.1 decrease by a factor of approximately 2. We can perform this division (4.3) for all global truncation errors in Table 4.1 to check how much they decrease. In Table 4.1, we can see that the forward Euler method has the convergence rate p ≈ 1, therefore it is first-order accurate. In Figure 4.1, we can see that the global truncation errors of the forward Euler method are on a line of slope 1, which shows that the global truncation errors of the forward Euler method are equal to O(h) which is proportional to h.
In Table 4.2, we can see that the global truncation errors of Heun’s method decrease by a factor of approximately 4. The first global truncation error decrease by a factor of approximately 4 because
1.1739 × 10−3
3.0109 × 10−4 ≈ 3.8988.
In Table 4.2, we can see that Heun’s method has the convergence rate p ≈ 2, therefore it is second-order accurate. Heun’s method is on a line of slope 2, therefore the global truncation errors of Heun’s method are equal to O(h2), which is proportional to h2 (see Figure 4.1).
For RK4, the global truncation errors decrease by a factor of approximately 16, which is in agreement with [4]. According to Atkinson et al. [4] we should expect that RK4 decrease by approximately 24 = 16. The first global truncation error decrease by a factor of approximately 16 (see Table 4.3) as follows
1.6252 × 10−6
1.0253 × 10−7 ≈ 15.851.
From Table 4.3 we can see that RK4 is fourth-order accurate since it has the con-vergence rate p ≈ 4. In Figure 4.1, we can see that the global truncation errors of RK4 are on a line of slope 4 and equal to O(h4) which is proportional to h4.
RK4 is more accurate than the forward Euler method and Heun’s method because the computed global truncation errors of RK4 are smaller than the global truncation errors of the forward Euler method and Heun’s method. The global truncation errors