# Novel Exponentially Fitted Two-Derivative Runge-Kutta Methods with Equation-Dependent Coefficients for First-Order Differential Equations

Authors: Yanping Yang; Yonglei Fang; Xiong You; Bin Wang
Journal: Discrete Dynamics in Nature and Society
 (2016)

---

## Abstract

The construction of exponentially fitted two-derivative Runge-Kutta (EFTDRK) methods for the numerical solution of first-order differential equations is investigated. The revised EFTDRK methods proposed, with equation-dependent coefficients, take into consideration the errors produced in the internal stages to the update. The local truncation errors and stability of the new methods are analyzed. The numerical results are reported to show the accuracy of the new methods.


---

## Body

## 1. Introduction

In this paper, we are interested in the numerical integration of initial-value problems (IVPs) of first-order differential equations in the form(1)y′x=fx,y,yx0=y0,whose solutions exhibit a pronounced exponential behaviour [1, 2]. For the numerical solution of problem (1), classical general-purpose methods such as Runge-Kutta (RK) methods and linear multistep methods (LMMs) can not produce satisfactory results due to the special structure of the problems. The possible ways to construct the numerical methods adapted to the character of the solutions can be obtained by using exponential fitting (EF) technique [3–9]. There have been a large number of results concerning this topic in [10–15].In standard approach to deriving exponentially fitted Runge-Kutta(-Nyström) methods, the effect of the error in the internal stages to the error of the final stage is completely neglected. D’Ambrosio et al. revisited the EFRK(N) method with two stages by considering the error in the internal stages in [1, 2]. Ixaru [16] introduced A-stable explicit fourth-order Runge-Kutta methods with three stages by evaluating the errors in each internal stage.For more accurate request, inspired by the work of Chan and Tsai [17], we improve the EFRK methods of D’Ambrosio et al. [1] by introducing the second derivative into the scheme. The paper is organized as follows. In Section 2, we construct two versions of exponentially fitted TDRK methods. The second version, revised version, has coefficients depending on the equation to be solved. In Section 3, we present the local truncations errors and analyze the linear stability of the new EFTDRK methods. Numerical experiments are reported in Section 4. Finally, in Section 5, we give some conclusive remarks.


## 2. Construction of the New Methods

For the numerical integration of (1), we consider the explicit TDRK methods of the form (2)Y1=yn,Yi=yn+cihfxn,yn+h2∑j=1i-1aijgxn+cjh,Yj,i=2,…,s,yn+1=yn+hfxn,yn+h2∑i=1sbigxn+cih,Yi,where g(x,y)≔∂/∂xf(x,y)+∂/∂yf(x,y)·f(x,y). Method (2) can be expressed briefly by the Butcher tableau: (3)
c
A
b
T
=
0
c
2
a
21
⋮
⋮
⋱
c
s
a
s
1
⋯
a
s
s
-
1
b
1
⋯
b
s
-
1
b
sAs scheme (2) indicates, at each step, this special explicit TDRK method involves only one evaluation of the function f and s evaluations of the function g. The order conditions for the general-purpose TDRK method (2) are given by Chan and Tsai [17]. In this section, we will consider the explicit two-stage TDRK method of the form (4)Y1=yn,Y2=yn+c2hfxn,yn+h2a21gxn,Y1,yn+1=yn+hfxn,yn+h2b1gxn,yn+b2gxn+c2h,Y2.In order to derive the exponentially fitted version of (4), we require that the internal stages exactly integrate the linear combination of the functions (5)1,x,expνx,and the final stage exactly integrates functions (6)1,x,exp⁡νx,xexp⁡νx.For the propose of constructing new methods, we introduce the operator (7)L2h,ayx=yx+c2h-yxn-hc2y′x-h2a21y′′x.With the local assumption that there is no error in Y1, we require that the linear operator (7) related to method (4) vanishes for the reference set (5) and obtain (8)L2h,a1=0,L2h,ax=0,L2h,aeνx=eνx+c2h-eνx-c2νheνx-ν2h2a21eνx=0.Solving the last equation, we have (9)a21=ec2z-1-c2zz2,where z=νh. As z approaches zero, a21 has the following Taylor series expansion: (10)a21=c222+c236z+c2424z2+c25120z3+c26720z4+c275040z5+⋯.When a21 is chosen, the second internal stage in (4) is exact for the linear combination of the functions in (5). Note that these three functions are solutions of the differential equation: (11)y′′′x-νy′′x=0.Then the leading term of the error in the computation of the second stage Y2 can be given by (12)err2yx=h3szy′′′x-νy′′x,where s(z) is the error constant that related to Y2. To determine s(z), following the ideas in [1, 2], we evaluate (7) and (12) for y=x2 yielding (13)L2h,ax2=h2c22-2a21,err2x2=-2szzh2.Letting L2[h,a]x2=err2(x2) gives (14)sz=2a21-c222z.On the other hand, for the standard exponentially fitted Runge-Kutta-type methods, the contribution of the error in the internal stages to the error of the overall numerical scheme is neglected. In order to construct a standard EFTDRK method, we introduce(15)LSh,byx=yx+h-yx-hfx,yx-h2b1Sy′′x+b2Sy′′x+c2h,Y2,where biS, i=1,2, are the weight coefficients of a standard EFTDRK method. Under the assumption that Y2=y(xn+c2h), we require the operator LS in (15) to vanish for all functions in (6) and obtain (16)b1S=-2+ez2+zc2-1-z1+c2+c2zc2z3,b2S=e-c2z2+ezz-2+zc2z3.We denote this method as SEFTDRK.SinceY2 approximates y(xn+c2h) with error (12), which is zero when y(x)=1,x,eνx, Y2 makes no contribution to the error in the final stage yn+1 for y(x)=1,x,eνx. However, in other cases, Y2≠y(xn+c2h). The main problem we are concerned with is how the error in Y2 affects the error in yn+1. So we consider (17)y′′xn+c2h=gxn+c2h,yxn+c2h=gxn+c2h,Y2+err2,and we have (18)y′′xn+c2h=gxn+c2h,Y2+err2gyxn+c2h,Y2+Oerr22.This leads to the following operator: (19)LRh,byx=yx+h-yx-hy′x-h2b1Ry′′x+b2Ry′′x+c2h-h3gyszy′′′x-νy′′x,where biR,i=1,2, are the weight coefficients of the revised EFTDRK method and gy=∂/∂yg(x+c2h,y(x)+c2hf(x,y(x))+a21h2g(x,y)). We compute (19) for y(x)=1,x,eνx,xeνx as follows: (20)LRh,b1x=xn=0,LRh,bxx=xn=0,LRh,beνxx=xn=ez-1-z-z2b1+b2ec2z,LRh,bxeνxx=xn=ez-1-2zb1R+2b2Rec2z-b2Rz2c2+m2sz,where m2=h2gy|x=xn. Solving the last two equations yields (21)b1R=ez-z-1z2+2ec2z2+z+ezz-2N,b2R=2z+2+ezz-2N,with N=2ec2z(c2z3-m2)+m2(2+c2z(2+c2z)). The obtained method is denoted as REFTDRK.Now we have obtained two versions (standard version and revised version) of exponentially fitted TDRK methods. In the next section, we will give some properties of the new methods.


## 3. Error and Stability Properties

### 3.1. Error Analysis

The local truncation errors (LTE) for the methods SEFTDRK and REFTDRK are given by(22)LTES=h44!y4xn-2c2gxxxn,yn-4c2gxyxn,ynfxn,yn-2c2gyyxn,ynfxn,yn2-2c2gyxn,yngxn,yn+2c2-1ν2gxxn,yn+2gyxn,ynfxn,yn-gxn,yn+Oh5,(23)LTER=h44!2c2-1ννy′′xn-2y3xn+y4xn+Oh5,respectively. For c2=2/3, the two methods are of algebraic order three. We denote them as SEFTDRK3 and REFTDRK3. If we take c2=1/2, then LTES=LTER=O(h5) and the two methods are of algebraic order four. We denote these two methods as SEFTDRK4 and REFTDRK4.


### 3.2. Stability Analysis

The stability of the TDRK methods can be analyzed by the linear test equation:(24)y′=λy,Reλ<0.For the test equation (24), we have m2=λ2h2. Applying the TDRK method (4) to the test equation (24), we obtain the recursive relation: (25)yn+1=Ru,zyn,where (26)Ru,z=1+u+b1V+b2Vu2+b2Vc2u3+b2Va21u4,u=λh,is called thestability function of the method, V=S for the standard version SEFTDRK and V=R for the revised version REFTDRK.Definition 1.
A three-dimensional region(27)Ω=Reu,Im⁡u,z∣Ru,z<1is called thestability region of the EFTDRK method with R(u,z) given by (26).In Figure1, we plot the stability regions of SEFTDRK4 and REFTDRK4 projected on the planes z=-1 and z=-2. In view of Figure 1, we observe that the revised version is clearly better than the standard one for stability with the same z.Figure 1
Stability regions of the SEFTDRK4 (a) and the REFTDRK4 (b) projected on the planesz=-1 and z=-2.
(a)
(b)


## 3.1. Error Analysis

The local truncation errors (LTE) for the methods SEFTDRK and REFTDRK are given by(22)LTES=h44!y4xn-2c2gxxxn,yn-4c2gxyxn,ynfxn,yn-2c2gyyxn,ynfxn,yn2-2c2gyxn,yngxn,yn+2c2-1ν2gxxn,yn+2gyxn,ynfxn,yn-gxn,yn+Oh5,(23)LTER=h44!2c2-1ννy′′xn-2y3xn+y4xn+Oh5,respectively. For c2=2/3, the two methods are of algebraic order three. We denote them as SEFTDRK3 and REFTDRK3. If we take c2=1/2, then LTES=LTER=O(h5) and the two methods are of algebraic order four. We denote these two methods as SEFTDRK4 and REFTDRK4.


## 3.2. Stability Analysis

The stability of the TDRK methods can be analyzed by the linear test equation:(24)y′=λy,Reλ<0.For the test equation (24), we have m2=λ2h2. Applying the TDRK method (4) to the test equation (24), we obtain the recursive relation: (25)yn+1=Ru,zyn,where (26)Ru,z=1+u+b1V+b2Vu2+b2Vc2u3+b2Va21u4,u=λh,is called thestability function of the method, V=S for the standard version SEFTDRK and V=R for the revised version REFTDRK.Definition 1.
A three-dimensional region(27)Ω=Reu,Im⁡u,z∣Ru,z<1is called thestability region of the EFTDRK method with R(u,z) given by (26).In Figure1, we plot the stability regions of SEFTDRK4 and REFTDRK4 projected on the planes z=-1 and z=-2. In view of Figure 1, we observe that the revised version is clearly better than the standard one for stability with the same z.Figure 1
Stability regions of the SEFTDRK4 (a) and the REFTDRK4 (b) projected on the planesz=-1 and z=-2.
(a)
(b)


## 4. Numerical Experiments

In this section, we will compare the numerical performance of the revised exponentially fitted TDRK methods with the standard exponentially fitted TDRK methods and a revised exponentially fitted RK method in the literature. The methods for comparison are denoted as follows.(i)
REFRK3: the revised EFRK method of order three given by D’Ambrosio et al. [1].(ii)
SEFTDRK4: the standard EFTDRK method of order four (c2=1/2) derived in this paper.(iii)
SEFTDRK3: the standard EFTDRK method of order three (c2=2/3) derived in this paper.(iv)
REFTDRK4: the revised EFTDRK method of order four (c2=1/2) derived in this paper.(v)
REFTDRK3: the revised EFTDRK method of order three (c2=2/3) derived in this paper.These five methods will be applied to the three test problems used by D’Ambrosio et al. [1]. Each problem will be integrated with differential step sizes and the comparison is based on the global errors at the right end point of the interval versus the step size.Problem 1.
Consider the linear initial-value problem:(28)y′x=λyx+kxk-1eλx,y0=0.The analytic solution is given by y(x)=xkeλx. We choose k=2 and λ=-2 and integrate the problem on the interval [0,100]. The numerical results are displayed in Table 1.Table 1
The global errors in the final step points obtained by several methods for Problem1.
Step size
REFRK3
SEFTDRK4
SEFTDRK3
REFTDRK3
REFTDRK4
1/8
6.41e − 05
1.33e − 05
6.60e − 05
4.84e − 05
3.29e − 06
1/16
7.42e − 06
7.44e − 07
6.46e − 06
5.50e − 06
1.97e − 07
1/32
8.96e − 07
4.42e − 08
7.11e − 07
6.55e − 07
1.19e − 08
1/64
1.10e − 07
2.69e − 09
8.33e − 08
7.99e − 08
7.40e − 10
1/128
1.36e − 08
1.66e − 10
1.00e − 08
9.87e − 09
4.60e − 11
1/256
1.69e − 09
1.03e − 11
1.24e − 09
1.22e − 09
2.87e − 12Problem 2.
Consider the nonlinear initial-value problem:(29)y′x=λy2x+2x3e2λxyx,y1=eλ.The analytic solution is given by y(x)=x2eλx. We choose λ=-4 and integrate the problem on the interval [1,100]. The numerical results are displayed in Table 2.Table 2
The global errors in the final step points obtained by several methods for Problem2.
Step size
REFRK3
SEFTDRK4
SEFTDRK3
REFTDRK3
REFTDRK4
1/8
1.89e − 05
2.21e − 06
2.60e − 06
2.58e − 06
4.94e − 07
1/16
2.27e − 06
1.12e − 07
4.49e − 08
2.07e − 07
1.91e − 08
1/32
2.77e − 07
6.32e − 09
8.91e − 09
2.18e − 08
9.31e − 10
1/64
3.42e − 08
3.75e − 10
1.77e − 09
2.58e − 09
5.08e − 11
1/128
4.25e − 09
2.28e − 11
2.59e − 10
3.04e − 10
2.96e − 12
1/256
5.29e − 10
1.41e − 12
3.46e − 11
3.74e − 11
1.79e − 13Problem 3.
Consider the nonlinear system:(30)y1′x=3y2x-x+λy12xx3eλx,y2′x=y2xx2+2y1x+λx2y2x-λx3x31+xeλx,y1=eλ,y21=1+eλ.The analytic solution is given by y1(x)=x3eλx and y2(x)=x(1+xeλx). We choose λ=-2 and integrate the problem on the interval [1,10]. The numerical results are displayed in Table 3.Table 3
The global errors in the final step points obtained by several methods for Problem3.
Step size
REFRK3
SEFTDRK4
SEFTDRK3
REFTDRK3
REFTDRK4
1/8
Inf
4.49e − 06
1.39e − 05
1.04e − 05
6.39e − 07
1/16
4.51e − 05
2.36e − 07
1.34e − 06
1.23e − 06
3.81e − 08
1/32
5.61e − 06
1.34e − 08
1.54e − 07
1.49e − 07
2.33e − 09
1/64
6.99e − 07
8.02e − 10
1.86e − 08
1.84e − 08
1.44e − 10
1/128
8.73e − 08
4.90e − 11
2.29e − 09
2.28e − 09
8.99e − 12
1/256
1.09e − 08
3.04e − 12
2.85e − 10
2.84e − 10
5.61e − 13In Problem1, the new third-order methods SEFTDRK3 and RERTDRK3 are comparable to the method SEFRK3 of D’Ambrosio et al., while in Problems 2 and 3 the former two are more accurate than the latter. The two new fourth-order methods SEFTDRK4 and RERTDRK4 are more accurate than the other three methods in all the three problems. Notice that, in Problem 3, with the step size h=1/8, the SEFRK3 produces enormous error while the other four new methods work very well.


## 5. Conclusions

In this paper, for integration of problems with solutions of the form of exponential functions, exponentially fitted two-derivative Runge-Kutta methods are considered. In contrast to the standard exponential fitting of RK methods, we take into account the contribution of the local error of internal stages to the update. The obtained coefficients thus depend on the derivatives of the function in the differential equation. Two standard two-stage EFTDRK methods of order three and order four as well as two standard two-stage EFTDRK methods of order three and order four are constructed. The results of numerical experiments show the effectiveness of the new methods.



---

*Source: 9827952-2016-04-20.xml*