A new way of deriving implicit Runge-Kutta methods based on repeated integrals
Hana Mizerová Katarína Tvrdá












A new way of deriving implicit Runge-Kutta methods
 based on repeated integrals
    Hana MizerováCorresponding author: hana.mizerova(at)fmph.uniba.sk. The work of H.M. was supported by the grants VEGA 1/0084/23 and VEGA 1/0709/24. 
The work of K.T. was supported by the grants KEGA 030STU-4/2023 and VEGA 1/0155/23.
 
		Katarína Tvrdá ^**


    
====================================================================================================================================================================================================================================================================







^* Department of Mathematical Analysis and Numerical Mathematics, Comenius University
Mlynská dolina, 842 48 Bratislava, Slovakia

^** Faculty of Civil Engineering, Slovak University of Technology in Bratislava
Radlinského 11, 810 05 Bratislava, Slovakia




Runge-Kutta methods have  an irreplaceable position among numerical methods designed to solve ordinary differential equations. Especially, implicit ones are suitable for approximating solutions of stiff initial value problems. We propose a new way of deriving  coefficients of implicit Runge-Kutta methods. This approach based on repeated integrals yields both new and well-known  Butcher's tableaux. We discuss the properties of newly derived methods and compare them with standard collocation implicit Runge-Kutta methods in a series of numerical experiments. In particular, we observe higher accuracy  and higher experimental order of convergence of some newly derived methods.




Keywords: implicit Runge-Kutta method, Cauchy's repeated integration, modified Newton-Cotes quadratures, stiff  initial value problem 






§ INTRODUCTION


Stiff problems form an important class of ordinary differential equations (ODEs). They appear in  many physical and chemical processes, for instance in fluid dynamics, chemical reaction kinetics, medicine, plasticity, neutron kinetics, porous media, gas transmission networks, transient magnetodynamics, in the study of spring and damping system,   and others, see <cit.> and the references therein.
 It is difficult to define stiffness, but it is even more challenging to solve stiff problems numerically - standard numerical methods are unstable, requiring unexpectedly small discrete step in regions of smooth solution curve.  
Runge-Kutta methods represent an important  class of well-established numerical  methods for solving ODEs. The explicit ones with small stability regions turned out to be unsuitable for stiff problems - discrete step size has to meet stability rather than accuracy requirements <cit.>. However, implicit  Runge-Kutta methods with large stability regions are frequently used to solve stiff problems despite being more demanding from the computational point of view. 

 


Let us consider an initial value problem (IVP) of the form


    y'(x)   =f(x,y(x)),   x∈ [a,b],
    
    y(a)   =y_a
 

with a<b, and f:[a,b]×ℝ→ℝ,  y_a ∈ℝ given.  The independent variable x not always denotes time, but it is usually referred to as a time variable. Quantity y represents the dependent variable.  Knowing the value y_a of a quantity y at time x=a we want to predict what happens to y as time evolves. Initial value problem is an essential tool in scientific computations and mathematical modelling. Therefore it is of a great interest to design efficient and reliable numerical schemes to approximate the exact solutions with a desired precision, especially in  the case of stiff systems which are pertinent for many real-world problems described by differential equations.

The famous Euler method published in 1768-1770 is based on a simple idea of a moving particle given by a differential equation: in a short period of time in which the velocity has not changed significantly, the change
in position will be approximately equal to the change in time multiplied by the initial velocity.
It was at the end of 19^ th and the beginning of 20^ th century when Runge, Heun and Kutta generalized the Euler method by allowing several evaluations of the derivative in one computational step. Nowadays, the Runge-Kutta methods (RK) form an important family of numerical schemes used to solve IVPs. See, e.g, <cit.> and the references therein. 

As already mentioned, stiff ODEs exhibit distressing behaviour when solved by classical numerical methods. Since explicit RK applied to stiff equations are usually unstable, and implementation of fully implicit RK is costly, there has been a growing interest in designing implicit methods with reduced computational costs. For instance, the so-called DIRK (diagonal implicit Runge-Kutta) and SDIRK (singly diagonal implicit Runge-Kutta) methods, see, e.g., <cit.> and the references therein. It is impossible to provide an exhausting literature overview on the design of new implicit Runge-Kutta(-type) methods. Let us mention, for instance, new optimized implicit-explicit Runge-Kutta methods <cit.>, new collocation methods based on weighted integrals  <cit.>, a two step fifth order RK for differential-algebraic equations <cit.>, usage of interval analysis tools to compute coefficients of Runge-Kutta methods <cit.>, and  implicit seven stage tenth order Runge-Kutta methods based on Gauss-Kronrod-Lobatto quadrature formula <cit.>.

In the present paper we provide a new way of deriving coefficients of implicit Runge-Kutta methods. Our approach yields both new and known schemes depending on the quadrature formula we choose.  In particular, we derive a stiffly accurate implicit four stage Runge-Kutta method with explicit first line of order 4. In addition, its experimental order of convergence in the case of linear IVPs is 6.   To the best of our knowledge there is no  result in the literature for computing  coefficients of Runge-Kutta methods based on repeated integration. 


The paper is organized as follows:  we provide preliminary material in Section 2. Section 3 contains the main result - new strategy of designing implicit Runge-Kutta methods. Several numerical experiments are presented in Section 4. Conclusion is followed by Appendix containing the list of newly derived and some well-known implicit Runge-Kutta methods.







§ PRELIMINARIES


We provide a necessary mathematical apparatus, notation and known results used in the paper.  




 §.§ A general Runge-Kutta method


Let y(x) be the exact solution of (<ref>). Let N≥ 1. We seek approximations  of exact values y() at points 

    x_n=a+nh, h=b-a/N, n=0,…,N.

Here h denotes the discrete step size.
The  initial condition (<ref>) yields y_0=y(x_0)=y(a)=y_a.

A general Runge-Kutta method is a single step method  employing  s stages K_i in one step. To determine value  from known value  we compute 

    = + h∑_i=1^s b_iK_i 
    
    K_i   =f(+h c_i,  +h∑_j=1^s a_ijK_j),  i=1,…, s.

The choice of coefficients a_i,j, b_i and c_i  defines the method itself. 
Every Runge-Kutta method can be  represented by a practical Butcher's tableau <cit.>, see Table <ref>.


Obviously, the structure of the matrix 𝐀 decides whether the corresponding Runge-Kutta method is explicit (lower triangular matrix with zeros on the diagonal),  DIRK (lower triangular matrix), SDIRK (lower triangular matrix with the same diagonal elements) or implicit (full matrix).
Typical examples are the explicit  and implicit Euler  methods,  Heun's method - explicit trapezoidal rule, implicit trapezoidal rule, and the prominent explicit fourth-order Runge–Kutta method (RK4). See Table <ref>.





  §.§.§ A-stability


The concept of A-stability introduced by Dahlquist, originally used for multistep methods, suggests to   solve IVP 

    y'(x)   =λ y(x)
    
    y(0)   =1
 
with λ∈ℂ, Re(λ)<0 numerically. 
 Desirably our numerical solution has the same behaviour as the exact solution  y(x)=e^λ x, namely

    → 0   n →∞.

Applying a general RK on (<ref>) we get, after one computational step, =R(hλ), and by induction, 
    =R(hλ)^n y_0.
 Thus the stability condition (<ref>) is equivalent to |R(hλ)|<1.

The function R(z) is called the stability function and the set S={z∈ℂ: |R(z)|≤ 1} is called the stability domain. 
If S ⊃ℂ^- (left half of complex plain), we say the method is A-stable.
The stability function of a general RK   reads 

    R(z)= det(𝐈 - z𝐀 +z𝐞𝐛^T)/ det(𝐈 - z𝐀),

where 𝐈, 𝐞 are the identity matrix and the vector of ones, respectively, see <cit.>. Figure <ref> depicts the stability domains of well-known RK from Table <ref>.






  §.§.§ Order conditions


The order of a Runge-Kutta method is p if and only if the local truncation error is 𝒪(h^p+1). There is no explicit Runge-Kutta method of order p with s=p stages for p>4.
The highest possible order of an implicit method with s stages is p=2s and is only attained by Gauss-Legendre methods. See <cit.> for more details on the so-called Butcher's barriers.

One of the great results of Butcher's theory for Runge-Kutta methods are the necessary and sufficient conditions to derive a new method of order p. These order conditions or Butcher's rules were firstly presented  and connected with the rooted tree theory by Butcher, see <cit.>. These conditions yield a system of equations with the unknowns a_i,j, b_i and c_i. 
The number of constraints for each order increases exponentially and is different for one- and high-dimensional problem. Table <ref> contains 17 order conditions <cit.> to be satisfied by a Runge-Kutta method of order p=5. Orders 6, 7 and 8 impose 37, 85 and 200 constraints, respectively.

For linear ODEs the number of order conditions can be reduced, see <cit.>. For instance, to obtain a sixth order method  it is enough to satisfy  16 conditions collected in Table <ref>.

A Runge-Kutta method is of order p if and only if all order conditions up to order p are satisfied.  

There are further constraints to be imposed in order to obtain a desired structure of the method, for instance:

  0em

  - explicit: a_i,j=0, ∀ j≥ 1

  - explicit first line: a_1,1=…=a_1,s=0

  - diagonal implicit: a_i,j=0, ∀ j> 1

  - singly diagonal: a_1,1=…=a_s,s

  - stiffly accurate: a_s,i=b_i, ∀ i=1, …, s

  - fully  implicit: a_i,j≠0, ∀ i,j=1, …, s.
  
Note that an additional consistency constraint (row-sum condition), c_i=∑_j=1^s a_i,j,  is typically taken into account, but it is not necessary, see, e.g., <cit.>.
We refer the interested reader to <cit.> for more details. 






 §.§ Standard implicit Runge-Kutta methods


In what follows we give a brief note on the derivation of standard implicit Runge-Kutta methods (sIRK) which shall be compared to our newly derived methods in terms of order and accuracy. 



  §.§.§ Collocation methods

Let us integrate the differential equation  (<ref>) over the interval [,x] to get 

    y(x)-y()=xf(x,y(x)).

We approximate the exact value y() by  and replace the integrand f(x,y(x)) by its unique Lagrange interpolation polynomial of degree at most s-1 corresponding to s points z_n,i=+τ_i h, 0≤τ_1<…τ_i < …<τ_s≤ 1, i.e.

    y(x)≈ +xL_s-1(x)=+∑_i=1^sf(z_n,i,y(z_n,i))x_i(x),

where

    L_s-1(x)=∑_i=1^s f(z_n,i,y(z_n,i))_i(x),   _i(x)=∏_i≠ j=1^s(x-z_n,j/z_n,i-z_n,j).

The principle of collocation method requires the above identity to hold at any point y(z_n,i). Thus, the approximations y_n,i≈ y(z_n,i) can be computed by a nonlinear system of equations

    y_n,i=+∑_j=1^sf(z_n,j,y(z_n,j))z_n,i_j(x),    i=1,…,s.

In case τ_s=1 we set =y_n,s, otherwise we set 

    =+∑_j=1^sf(z_n,j,y(z_n,j))_j(x).

Denoting K_i=f(z_n,i,y(z_n,i)),  c_i=τ_i, i=1,…, s, and

    b_i   =z_n,s_i(x),    
    a_i,j=z_n,i_j(x)

we get a standard implicit Runge-Kutta method  with s stages. It shall be referred to as sIRKs. For reader's convenience we list the corresponding Butcher's tableaux  for s=2,3,4,5 in Appendix, see Tables <ref> and <ref>.


Choosing equally spaced points z_n,i=+hi-1/s-1, i=1,…,s, h=-, means that coefficients b_i are computed using a (closed) Newton-Cotes quadrature formula on the interval [,]. 




  §.§.§ General implicit Runge-Kutta methods


Every collocation method is a Runge-Kutta method but not every Runge-Kutta method is a collocation method. A general IRK can be derived using  the so-called simplifying order conditions introduced by Butcher <cit.>. In particular, the derivation of the coefficients b_i and a_i,j relies on three of them:

    B(p)   =∑_i=1^s b_i c_i^k-1=1/k,   k=1, …, p
    
    C(q)   =∑_j=1^s a_i,j c_j^k-1=c_i^k/k,   k=1, …, q,  i=1…,s
    
    D(r)   =∑_i=1^s  b_i c_i^k-1a_i,j=b_j/k(1-c_j^k),   k=1, …, r,  j=1…,s.

The coefficients c_i are chosen according to the quadrature formula as specified below. The following theorem was proved in <cit.>.

If a Runge-Kutta method satisfies B(p), C(q) and D(r) with p≤ q+r+1 and p≤ 2q+2, then its order of convergence is p. 



 
There are three groups of general IRK stemming from the choice of quadrature formula:


  i) Gauss-Legendre methods which are collocation methods based on the roots of the Legendre polynomial P^*_s(x) shifted to [0,1]. The coefficients b_i and a_i,j are obtained from B(s) and C(s), respectively. It can be shown that they are of maximal order  2s, where s is the number of stages. 

  ii) Radau methods are divided to two subgroups: Radau IA based on the Radau left quadrature formula: c_1=0 and remaining c_i are the roots of P_s^*(x)+P_s-1^*(x),  and Radau IIA  based on the Radau right quadrature formula:  c_s=1 and remaining c_i are the roots of P_s^*(x)-P_s-1^*(x). Moreover, Radau IA impose conditions B(s) and D(s), while Radau IIA require conditions B(s) and C(s) to be satisfied. The well-known example of Radau IIA for s=1 is the implicit Euler method. The Radau methods are not collocation methods and are all of order  2s-1.

  iii) Lobatto methods, also called Lobatto III, are based on the choice c_1=0, c_s=1 and remaining c_i being the roots of the polynomial P_s^*(x)-P_s-2^*(x). Nowadays there are three established subgroups of Lobatto methods based on the choice of the coefficients a_i,j:Lobatto IIIA impose condition C(s) and for s=2 yield the implicit trapezoidal method, Lobatto IIIB impose D(s), and Lobatto IIIC impose conditions a_i,1=b_i for i=1,…, s and C(s-1).  They are all of order 2s-2.

We refer the reader to <cit.> for more details.  In Appendix we give Butcher's tableaux for selected Gauss-Legendre, Radau and Lobatto methods, see Tables <ref> and <ref>.




 §.§ Repeated integrals


We aim to derive coefficients for RK in order to approximate the  solution to IVP (<ref>).  Our new strategy  requires the moments of y'(x) and f(x,y(x)) to be equal. It means that besides the standard integral identity typical for methods based on numerical integration, 

    y'(x)   = f(x,y(x))

we want the following identities to be satisfied:

    (-x)y'(x)   = (-x)f(x,y(x))
    …   ……
    (-x)^γ-1 y'(x)   = (-x)^γ-1 f(x,y(x))

for arbitrary γ≥ 1 to be chosen later. As we shall see, it is useful to rewrite these moments equations as repeated integrals. This is possible thanks to 
Cauchy's formula for repeated integration <cit.> which allows us to compress γ antiderivatives of one function to a single integral. See <cit.> for the proof.
[Cauchy's repeated integration formula]
Let g be a continuous function on the real line. Then the γ^ th repeated integral of g, denoted by g^(-γ), is given by a single integral,

    g^(-γ):=∫_a^b∫_a^x_1…∫_a^x_γ-1 g(x_γ) dx_γ…dx_1 = 1/(γ-1)!∫_a^b (b-t)^γ-1g(t)  dt.


 From now on, to avoid confusion with the discretization points x_n, n=0,…,N, we use x for each integration variable on the left-hand side of (<ref>).

By the Newton-Leibniz formula (the fundamental theorem of calculus) we have

    y(x)-y()=xy'(x).

Due to (<ref>) identities (<ref>) for γ≥ 2 become 

    x…xy'(x)_(γ-1)-   =x…xf(x,y(x))_(γ-1)-.

In view of (<ref>) the integral on the left-hand side is nothing but 

    x…xy'(x)_(γ-1)-   =
    x…xy(x)_(γ-2)--y()x…x1_(γ-2)-
       =x…xy(x)_(γ-2)--y()h^γ-1/(γ-1)!, h=-.
  
Combining the above with (<ref>) we get 

    h y()+xf(x,y(x))   =y(x)
    h^2/2y()+xxf(x,y(x))   =xy(x)
       …
    h^S/S!y()+x…xf(x,y(x))_S-   =x…xy(x)_(S-1)-

 for S=γ-1≥ 1. 



  §.§.§ Numerical quadrature for repeated integrals


Another step towards deriving coefficients of new implicit RK is a suitable numerical approximation of repeated integrals in (<ref>). We refer to  <cit.>, where the detailed derivations of such quadrature formulas were presented. Here we briefly recall the results. 

Gauss quadrature for  repeated integral.
Approximation of a repeated integral g^(-S) by means of a general Gauss quadrature formula of order 2m-1  reads

    g^(-S)≈(b-a)^S/S!∑_j=1^m w_S,j g(t_j),
 
where  

    t_j=(b-a)(x_j+1)+2a/2,      w_S,j=S w_j (1-x_j)^S-1/2^S

are the nodal Gaussian points shifted from [-1,1] to [a,b], and the weights depending on the number of repeated integrals S, respectively. 
 Here x_j ∈ [-1,1] are the original nodal points and w_j the corresponding weights of Gauss quadrature formula on [-1,1]. The reader might have noticed that (<ref>) with (<ref>) stems from Cauchy's formula (<ref>). See <cit.> for more details. 

 
Newton-Cotes quadrature for  repeated integral. 
 Another possibility to approximate g^(-S) is by the so-called modified Newton-Cotes formulas introduced in <cit.>. In this case, we firstly approximate function g by its unique Lagrange interpolation  polynomial L_m-1(x) corresponding to equidistant nodal points and then compute the weights directly, i.e. without applying Cauchy's formula (<ref>).
 
  For the case of closed Newton-Cotes formula let H=b-a/m-1 and T_j=a+(j-1)H, j=1,…,m. Then we write 

    g(x)≈ L_m-1(x)=∑_j=1^m g(T_j)ℓ_j(x),   ℓ_j(x)=∏_j≠ i =1^m (x-T_i/T_j-T_i),
 
where ℓ_j(x) are the basis Lagrange interpolation polynomials.
Consequently, for the repeated integral we obtain

 
    g^(-S)≈∑_j=1^m g(T_j)∫_a^b∫_a^x…∫_a^xℓ_j(x) dx…dx _(S-1)-= H^S ∑_j=1^m W_S,j g(T_j),   H^S W_S,j=ℓ_j^(-S).
 
An analogous formula holds for open Newton-Cotes formulas which do not include the endpoints a and b. Indeed, 
 
    g^(-S)≈ H^S ∑_j=1^m-2 W_S,j g(T_j),   H^S W_S,j=ℓ_j^(-S),    T_j=a+jH, j=1,…,m-2
 

with the interval [a,b] being divided into m-1 subintervals of equal length.
 







§ NEW FAMILIES OF IMPLICIT RUNGE-KUTTA METHODS


Our motivation is to solve (stiff) initial value problem (<ref>) by an implicit Runge-Kutta method (IRK). We propose a new way of deriving the entries of 𝐀, 𝐛 and 𝐜, cf. Table <ref>. Instead of focusing on simplifying order conditions (<ref>) we employ the moment approach and repeated integrals as described in Subsection <ref> to determine the coefficients a_i,j, b_i and c_i. 

We proceed in four steps:

0em


  (Step 1) choose a numerical quadrature  and determine c_i accordingly

  (Step 2) compute b_i by the chosen quadrature formula 

  (Step 3) for the number of unknown stages write the corresponding number of identities (<ref>)  and apply chosen quadrature formula on repeated integrals

  (Step 4) solve the resulting  system of algebraic equations for unknowns y_n,i and get a_i,j  
 
We shall follow these steps and detail the procedure for three different cases: the cases of closed and open Newton-Cotes  quadrature formulas (NC) as well as the case of Gaussian quadratures of Legendre-, Radau- and  Lobatto-type. 



 §.§ New way based on closed Newton-Cotes formulas


(Step 1) For a closed NC with the nodal pints z_n,i=+hi-1/s-1, i=1,…,s, we have 
    c_i=i-1/s-1,  i=1,…,s.


(Step 2) Integrating (<ref>) over [,] and applying the closed NC from Step 1 to approximate  the integral of f(x,y(x)) over [,]  we get 

    =+h∑_i=1^s b_i K_i,    b_i=ℓ_i(x),  i=1,…,s.

Recall that  K_i=f(z_n,i,y(z_n,i)). 
 The coefficients b_i are simply the integrals of the corresponding Lagrange basis interpolation polynomials, analogously as in the case of standard collocation IRK, cf. Remark <ref>. Indeed, replacing f(x,y(x)) by its unique Lagrange interpolation polynomial corresponding to the nodal points from  Step 1 we get

    y()=y()+ f(x,y(x))≈= + ∑_i=1^s f(z_n,i,y(z_n,i)) ℓ_i(x).

The difference compared to collocation methods from Subsubsection <ref> lies in deriving the entries of the matrix 𝐀. 

(Step 3) Let s≥ 3. Since z_n,1= and z_n,s=, we  have S=(s-2) unknown stages. Thus we consider the first S equations of (<ref>) to which we apply the modified Newton-Cotes quadrature (<ref>) with H=h/s-1. 
This results in a system of S equations with the unknowns y(z_n,2), …, y(z_n,s-1)), 

    h+ H^2∑_i=1^s W_2,iK_i   =H∑_i=1^s W_1,i y(z_n,i)
    h^2/2+H^3∑_i=1^s W_3,iK_i   =H^2∑_i=1^s W_2,iy(z_n,i)
       …
    h^S/S!+H^S+1∑_i=1^s W_S+1,iK_i   =H^S∑_i=1^s W_S,iy(z_n,i),

where K_i=f(z_n,i,y(z_n,i)) are parameters.
 Note  that y(z_n,1)=y()≈ is known from the previous computational step and y(z_n,s)=y()≈ is given by (<ref>) in Step 2.
 
 
(Step 4) Each y(z_n,i), i=2,…,s-1 as a solution to  (<ref>) taking into account (<ref>) can be expressed in the form

    y(z_n,i)=+h∑_j=1^s a_i,jK_j,  i=2, …, s-1,

from which we get the coefficients a_i,j, i=2,…,s-1, j=1,…,s.
Due to y(z_n,1)=y()≈ we have a_1,j=0, j=1,…,s, and from (<ref>) simply a_s,j=b_j, j=1,…,s.



For 𝐬=2 we obviously get the Lobatto IIIA with 2 stages and order 2, the so-called implicit trapezoidal method. Since K_1=f(,y()) and K_2=f(,y()), there are no unknown stages and thus no need to use repeated integrals (<ref>). 

For 𝐬=3 there is only one unknown stage K_2=f(z_n,2,y(z_n,2)) with z_n,2=+1/2h. Thus, from  the first equation of (<ref>) we get 

    h + 2H^2/3(K_1+2K_2)=H/3(+2y(z_n,2)+),  H=h/2.

In addition, (<ref>) results in

    = + h/6(K_1+4K_2+K_3 ).

Therefore, 

    y(z_n,2)    = (3/2-1/8 -1/4)+ h/4(K_1+2K_2)- h/24(K_1+4K_2+K_3 ) 
        = +h/24(5K_1+8K_2-K_3 ).

We conclude

    𝐀=([     0     0     0;  5/24  8/42 -1/24;   1/6   2/3   1/6;       ]),  𝐛^T=(1/6,2/3,1/6),
       𝐜=(
     [   0; 1/2;   1;     ]).

This method is again known as Lobatto IIIA.

For 𝐬=4 we deal with the well-known 3/ 8-Simpson rule.
It means we divide the interval [,] into three subintervals of equal length H=h/3. Then 

    f(x,y(x))≈3/8 H (K_1+3K_2+3K_3+ K_4),
 
where z_n,i=x_n+(i-1)H, i=1,2,3,4. Thus c_i=i-1/3, i=1,2,3,4, and
Step 2 yields 

    = +  3/8 H (K_1+3K_2+3K_3+ K_4).

It means b_1=1/8=b_4 and b_2=3/8=b_3. 
Since the stages K_2 and K_3 including the values y(z_n,2) and y(z_n,3) are at this point unknown, we need two additional equations to close the system. Therefore we consider the first two identities of (<ref>) to get

    3H   +3/40 H^2 (13K_1+36K_2+9K_3+ 2K_4)=3/8 H(+3y(z_n,2)+3y(z_n,3)+)
    9H^2/2   +9/80 H^3 (12K_1+27K_2+0 K_3+ K_4 )=3/40 H^2 (13+36y(z_n,2)+9y(z_n,3)+ 2).

Let us denote α_j∈{1,3,3,1}, β_j∈{13,36,9,2}, γ_j∈{12,27,0,1}.
 Using (<ref>) we further rewrite the latter equations to get 


    y(z_n,2)   = 2 + 2/3 H ∑_i=1^4(β_i/10-3α_i/16)K_i-y(z_n,3)
    
     y(z_n,3)   = 5 + 1/6 H ∑_i=1^4(γ_i-α_i/2)K_i-y(z_n,2).


Now we substitute (<ref>) into (<ref>) to obtain

    y(z_n,2)   =  + H/18∑_i=1^4(α_i/4-2β_i/5+γ_i)K_i,

and analogously, (<ref>) into (<ref>) to get

    y(z_n,3)   =  + H/18∑_i=1^4(-5α_i/2+8β_i/5-γ_i)K_i.

Now, let 
    A_i=α_i/4-2β_i/5+γ_i∈{141/20, 267/20, -57/20, 9/20},     B_i=-5α_i/2+8β_i/5-γ_i∈{63/10, 231/10, 69/10, -3/10}.

Hence (<ref>), (<ref>) and (<ref>) become 

    =+3H/8∑_i=1^4α_iK_i,    
    y(z_n,2)=  + H/18∑_i=1^4A_iK_i,    
    y(z_n,3)=  + H/18∑_i=1^4B_iK_i.

Recall H=h/3 and
K_i=f(z_n,i,y(z_n,i)).  We have derived four stage IRK given by

    =+ h ∑_i=1^4α_i/8K_i
    
    K_1   =f(,)
    
    K_2   =f(+h/3, + h  ∑_i=1^4 A_i/54K_i) 
    
      K_3   =f(+2/3h, + h∑_i=1^4B_i/54K_i)
    
       K_4   =f(+h,  + h ∑_i=1^4α_i/8K_i) .

We shall refer to it as nIRK4. Note that Lobatto IIIA with 4 stages differs in the choice of c_i and sIRK4 differs in the coefficients a_2,j, a_3,j, j=1,2,3,4. 

[properties of  nIRK4]
Newly derived four stage implicit Runge-Kutta method  given by Butcher's tableau 

    [       0       0       0       0       0;     1/3  47/360  89/360 -19/360   3/360;     2/3  21/180  77/180  23/180  -1/180;       1     1/8     3/8     3/8     1/8;             1/8     3/8     3/8     1/8;         ]

 is  stiffly accurate with explicit first line. It is  A-stable with the order of convergence p=4.    


Obviously, the first line is explicit and since a_s,j=b_j, j=1,…,s the method is stiffly accurate according to <cit.>. 
The method is indeed A-stable. The stability function  reads

    R(z)= -z^3+ 12z^2+ 60 z +120/z^3- 12z^2+ 60 z-120,

and direct calculation for z=u+iv ∈ℂ,  yields

    |R(z)| ≤ 1 ⇔   u(u^4+v^4+2u^2v^2+70u^2+30v^2+600) ≤ 0

which is only true for u≤ 0. Thus the stability domain is  S={z∈ℂ:  Re(z)≤ 0}.
It is easy to see that the method satisfies the simplifying order conditions B(p), C(q) and D(r) with p=4, q=3 and r=0, which by Theorem <ref> means the order is p=4.
Another possibility is to check that all 8 order conditions from Table <ref> up to order p=4 are satisfied. 


For 𝐬≥ 5 one can analogously as above derive IRK  based on closed Newton-Cotes formulas. We give  Butcher's tableau and the corresponding order of nIRK5 in Table <ref>.  







 §.§ New way based on open Newton-Cotes formulas


The only difference in deriving the coefficients of IRK by an open Newton-Cotes formula compared to the closed one is that the end points of the interval [,] are not included. It means the resulting method neither has the explicit first line nor is stiffly accurate according to <cit.>.

(Step 1) In order to derive an s stage Runge-Kutta method, we need to consider s+1 subintervals of [,] with the corresponding nodal points z_n,i=+hi/s+1, i=1, …, s. Obviously, 
    c_i=i/s+1, i=1,…, s.


(Step 2) By the same reason as before, it holds that 

    =+h∑_i=1^s b_i K_i,    b_i=ℓ_i(x),  i=1,…,s.

Note the difference compared to (<ref>) lies in ℓ_i(x) being polynomials defined for different nodal points and   of different degree.

(Step 3)  For s≥ 2 we have S=s unknown stages and the corresponding first s equations of (<ref>) to be solved. Note that y(z_n,1)≠ y() and y(z_n,s)≠ y(). We approximate the integrals in (<ref>) by the modified open Newton-Cotes quadrature formula (<ref>) which results in 

    h+ H^2∑_i=1^s W_2,iK_i   =H∑_i=1^s W_1,i y(z_n,i)
    h^2/2+H^3∑_i=1^s W_3,iK_i   =H^2∑_i=1^s W_2,iy(z_n,i)
       …
    h^S/S!+H^S+1∑_i=1^s W_S+1,iK_i   =H^S∑_i=1^s W_S,iy(z_n,i),

where K_i=f(z_n,i,y(z_n,i)) are parameters, and y(z_n,1),…,y(z_n,s) are the unknowns.
Note that y()≈ is known from the previous computational step and H=h/s+1. 

(Step 4) We solve system (<ref>) and as a result we obtain the coefficients a_i,j, i,j=1,…,s from 
  
    y(z_n,i)=+h∑_j=1^s a_i,jK_j,  i=1, …, s.


For s=3,4 the newly derived IRK are listed in Table <ref> and are referred to as nIRK3o and nIRK4o. They are both A-stable and of order p=4.  



Dividing the interval [,] into s subintervals of equal length H=h/s yields either (s+1) or (s-1) stage method for closed or open NC, respectively. 
 



  §.§.§ New way based on Newton-Cotes formulas combined with Cauchy's formula


As we have seen, weights in modified Newton-Cotes quadrature formulas (<ref>) are obtained as repeated integrals of the basis Lagrange interpolation polynomials. However, we can also rewrite them using Cauchy's repeated  integration formula (<ref>) and then compute the weights.  
When we use this approach in (<ref>) in Step 3 we derive new yet different implicit Runge-Kutta methods. Their Butcher's tableaux and corresponding orders for s=2,3,4,5,6 in the case of closed formulas and for s=3,4 in the case of open formulas are given in Tables <ref>,<ref> and <ref>. They are refereed to as nIRKsc and nIRKsoc, respectively. Here “c” stands for “Cauchy”. 





 §.§ New way based on Gauss quadratures



It turns out that our approach based on repeated integrals combined with Gauss quadratures yields the well-known implicit Runge-Kutta methods. For completeness, let us briefly summarize the results. 



  §.§.§ Gauss-Legendre quadrature

 
 (Step 1)
 The choice of c_i is obvious:  the roots of the Legendre polynomial P_s^*(x) shifted to [0,1]. 

 (Step 2) 
 The coefficients b_i are exactly the weights w_i of the corresponding Gauss-Legendre quadrature. Indeed,  

    f(x,y(x))≈ h∑_i=1^s f(z_n,i,y(z_n,i)) w_i =h∑_i=1^s b_i K_i.

Recall K_i=f(z_n,i,y(z_n,i)) with z_n,i=+c_i h, i=1,…,s. 

 (Step 3)
Let s≥ 2. We realize that z_n,1≠ and z_n,s≠. Thus, we have S=s unknown stages to be determined in order to compute  given . To this end, we consider the first s identities of (<ref>). The repeated integrals of f(x,y(x)) and y(x) are now approximated by means of (<ref>). The system reads

    h+ h^2/2!∑_i=1^s w_2,iK_i   =h∑_i=1^s w_1,i y(z_n,i)
    h^2/2+h^3/3!∑_i=1^s w_3,iK_i   =h^2/2!∑_i=1^s w_2,iy(z_n,i)
       …
    h^s/s!+h^s+1/(s+1)!∑_i=1^s w_s+1,iK_i   =h^s/s!∑_i=1^s w_s,iy(z_n,i).


 (Step 4) Solving the underlying system yields y(z_n,1), …, y(z_n,s), in the form 
  
    y(z_n,i)=+h∑_j=1^s a_i,jK_j,  i=1, …, s,

 and finally we have the coefficients a_i,j, i,j=1,…,s.




  §.§.§ Gauss-Radau quadrature


(Step 1)
 Set c_1=0 and c_i stemming  from left Radau quadratures or set c_s=1 and c_i stemming  from right Radau quadratures as described above in Subsubsection <ref>.

(Step 2) Analogously as for the Gauss-Legendre quadrature, the coefficients b_i coincide with the weights of Radau left quadrature formula or Radau right quadrature formula, respectively.

 (Step 3) In both cases we have S=s-1 unknown stages. We consider the first s-1 equations of (<ref>) where the repeated integrals are approximated by (<ref>) using the corresponding nodal points and weights of left or right Radau quadratures, respectively. 

(Step 4)
Solving the system of s-1 first equations of (<ref>) with unknowns y(z_n,2), …, y(z_n,s) in the case of left Radau quadrature, or with unknowns y(z_n,1), …, y(z_n,s-1) in the case of right Radau quadrature, yields the coefficients a_i,j, i,j=1,…,s. Note that  a_1,j=0 for left Radau quadrature
and a_s,j=b_j  for right  Radau quadrature formulas.



  §.§.§ Gauss-Lobatto quadrature

 
 (Step 1)
We set c_1=0,  c_s=1 and c_i, i=2,…,s-2, to be the roots  of  P_s^*(x)-P_s-2^*(x). 

 (Step 2)
The coefficients b_i coincide with the weights of chosen Lobatto quadrature formula, and again we have

    =+h∑_i=1^s b_i K_i.


 (Step 3)
We need to consider S=s-2 equations of (<ref>) together with (<ref>) on noting that  is known from the previous computational step. We approximate the repeated integrals by means of (<ref>) using the nodal points and weights corresponding to chosen Gauss-Lobatto quadrature formula. 

 (Step 4)
Coefficients a_i,j, i,j=2,…,s-2, are obtained from 
 
    y(z_n,i)=+h∑_j=1^s a_i,jK_j,  i=2, …, s-2

which solve the system from Step 3  taking into account (<ref>).
 Since z_n,1= we have the explicit first line,  a_1,j=0, j=1,…,s. Finally, a_s,j=b_j, j=1,…,s because z_n,s=. 

The resulting IRK of Gauss-type (nIRK-Gs,  nIRK-RIs, nIRK-RIIs, nIRK-Ls according to type of Gauss quadrature) are listed in Tables <ref> and <ref> and are classified according to existing families of general Gauss-type IRK.


Note that the advantage of our proposed approach compared to the derivation of general IRK based on Gauss quadratures as described in Subsubsection <ref> lies in the computation of a_i,j and b_i. For instance, in the case of Gauss-Legendre quadrature, our strategy directly yields b_i and requires solving system of  s equations of (<ref>) to determine a_i,j. The general approach requires solving s^2+s equations stemming from the simplifying order conditions B(s) and C(s). 


       




§ NUMERICAL EXPERIMENTS


We present numerical experiments for selected newly derived implicit Runge-Kutta methods and compare relative errors and experimental order of convergence (EOC) with standard collocation implicit Runge-Kutta methods  with the same number of stages. 

We compute  relative error

    e:=max_n=1,…,N|y(x_n)-y_n|/max_n=1,…,N|y(x_n)|

for different values of N, and  EOC given by 

    EOC_N,2N:=log_2(e^N/e^2N),

where e^N is the relative error computed with the discrete step size h=1/N and e^2N is the relative error computed with half the discrete step size. 




 §.§ Experiment 1: linear stability


Firstly, we consider standard IVP to test linear stability of Runge-Kutta methods,

    y'(x)   =-15y(x),  x ∈ [0,1]
    
    y(0)   =1

with the exact solution y(x)=e^-15x. 

In Table <ref> we present relative errors and EOC for the following four stage stiffly accurate implicit Runge-Kutta methods with the explicit first line:  

nIRK4 - new IRK  based on closed Newton-Cotes formula,
 nIRK4c - new IRK based on closed Newton-Cotes formula using Cauchy's formula,
 sIRK4 - standard collocation IRK with τ_i=i/4, i=0,…,4.


We term them “closed” because they include the approximation of function values at the end points,  and . 




 Table <ref>  contains relative errors and EOC for the following four stage fully implicit Runge-Kutta methods:

nIRK4o - new IRK  based on open Newton-Cotes formula,
 nIRK4oc  - new IRK  based on open Newton-Cotes formula using Cauchy's formula,
 sIRK4o - standard collocation IRK with τ_i=i/5, i=1,…,4.

Obviously, we refer to these methods as “open” because, as opposed to the “closed” ones, their stages do not include the end points  and . 

The graphs of relative errors in logarithmic scale for the four stage implicit Runge-Kutta methods from Tables <ref> and <ref> are depicted in Figure <ref>.
 



In Table <ref> we present the results obtained  by selected five stage stiffly accurate implicit Runge-Kutta methods with the explicit first line:  

nIRK5 - new IRK  based on closed Newton-Cotes formula,
 nIRK5c - new IRK based on closed Newton-Cotes formula using Cauchy's formula,
 
sIRK5 - standard collocation IRK with τ_i=i/5, i=0,…,5. See the left plot in Figure <ref>. 
   

Finally, Table <ref> contains the relative errors and EOC for selected three stage fully implicit Runge-Kutta methods:

nIRK3o - new IRK  based on open Newton-Cotes formula,
 nIRK3oc  - new IRK  based on open Newton-Cotes formula using Cauchy's formula,
 
sIRK3o - standard collocation IRK with τ_i=i/4, i=1,…,3.  See  the right plot in Figure <ref>.




Note, in particular, the EOC of our newly derived methods nIRK4, nIRK5 and nIRK3o is 6, 8 and 6, respectively, whereas the theoretical orders of convergence are 4,6 and 4, respectively.  The reason in the case of nIRK4 might be that this method satisfies 13 out of 16 order conditions for p=6 corresponding to linear ODEs from Table <ref>. The first sum  for p=5 and the first two for p=6 only yield 
   values 11/54, 19/108 and 1/540 instead of 1/5, 1/6 and 1/30, respectively. Similarly, nIRK3o satisfies 11 out of 16 order conditions from Table <ref>.
   



 §.§ Experiment 2: stiff equation


Solving stiff equations is difficult for many numerical methods. We test our newly derived implicit Runge-Kutta methods on the stiff initial value problem presented in <cit.> and compare the results with the standard implicit Runge-Kutta methods with the same number of stages.

Let us consider the IVP

    y'(x)   =-100y(x)+99e^2x,  x ∈ [0,0.5]
    
    y(0)   =0

with the exact solution 
    y(x)=33/34(e^2x-e^-100x).

We compute the errors given in (<ref>)
for the values N=8,16,32,64,128, and the corresponding EOC as defined in (<ref>).

Table <ref> and Figure <ref> indicates better performance of our newly derived implicit Runge-Kutta method with four stages nIRK4  compared to the standard implicit collocation method sIRK4 when applied to the stiff initial value problem (<ref>).   EOC of our newly derived method nIRK4 is 6 compared to 4  of sIRK4.
Moreover, we can observe that our “open” method nIRK4o has the same EOC=4 as sIRK4o. In particular, both newly derived methods approximate the exact solution with smaller errors for each N.



   
    





 §.§ Experiment 3: another stiff equation

In the third experiment we consider the IVP

    y'(x)   =-50(y(x)-cosx),  x ∈ [0,1]
    
    y(0)   =0

with exact solution 
    y(x)=50(sinx-50cosx)/2501-2500e^-50x/2501.

We perform computations for the values N=8,16,32,64,128.

Relative errors and EOC presented in Table <ref> confirm the same behaviour of methods as observed 
in Table <ref> for the stiff IVP considered in Experiment 2. Namely, newly derived “closed” methods nIRK4 and nIRK5 approximate the exact solution with smaller relative errors and higher EOC (6 and 8) as compared to the standard collocation methods sIRK4 and sIRK5 (with EOC 4 and 6). See also Figure <ref>.
  



  
   


  


 §.§ Experiment 4: nonlinear stiff equation

In this experiment we consider a simple model of flame propagation, 

    y'(x)   =y^2(x)-y^3(x),  x ∈[0,2/δ] 
    
    y(0)   =δ,

where δ <<0. 
The exact solution is given by 

    y(x)=1/W(A e^A-x)+1,

where A=1/δ-1 and W(t) is the so-called Lambert W function  <cit.>. It solves the equation W(t) e^W(t)=t. 
We compute the above IVP by nIRK4, sIRK4 and nIRK5, sIRK5 for δ=0.01. Figure <ref> depicts graphs of exact and numerical solutions computed for N=16. The relative errors and EOC for N=2^k, k=3,4,5,6 are presented in Table <ref>.  See also Figure <ref>. We can observe smaller relative errors of our newly derived IRK compared to standard collocation IRK. The EOC confirms theoretical order of convergence, p=4 for selected four stage  methods and p=6 for selected five stage methods. 
   
  


   
 



§ CONCLUSION

We have proposed the new way of deriving coefficients of implicit Runge-Kutta methods used for numerical  approximation of initial value problems (<ref>). Instead of standard collocation approach or the general approach based on simplifying order conditions as described in Subsection <ref> we have imposed  identities for  the  moments of y'(x) and f(x,y(x)).  The number of resulting equations containing repeated integrals depends on the number of stages s and the numerical quadrature being used, see Section <ref>.
We have derived several Butcher's tableaux based on closed and open modified Newton-Cotes quadrature formulas with and without employing Cauchy's repeated integration formula, as well as  based on Gauss-Legendre, -Radau and -Lobatto quadrature formulas. Various numerical experiments performed for selected newly derived implicit Runge-Kutta methods have confirmed their theoretical order of convergence indicated in Tables <ref> and <ref>. Moreover, the latter methods are A-stable and yield approximation of solutions  to considered linear and  nonlinear stiff initial value problems with higher accuracy than the standard collocation implicit Runge-Kutta methods with the same number of stages. Especially, recall nIRKs with EOC 2s-2 for linear IVPs which is typical for Lobatto methods. While nIRK2 and nIRK3 indeed coincide with Lobatto IIIA methods, nIRK4 and nIRK5 do not belong to family of Lobatto methods.

10
ahmad
R. R. Ahmad, N. Yaacob, A.-H. Mohd Murid.
Explicit methods in solving stiff ordinary differential equations.
Int. J. Comput. Math. 81(11): 1407–1415, 2004.

BuFa
R. L. Burden, J. D. Faires.
Numerical Analysis, 5th edition.
PWS-Kent Publishing Company, 1993.

butcher_history
J. C. Butcher.
A history of Runge-Kutta methods.
Appl. Numer. Math. 20: 247–260, 1996. 

butcher_oc
J. C. Butcher.
Coefficients for the Study of Runge-Kutta Integration Processes.
J. Aust. Math. Soc 3: 185–201, 1963.

butcher_radau
J. C. Butcher.
Integration Processes Based on Radau Quadrature Formulas. 
Math. Comput. 18: 233–244, 1964.

butcher_soc
J. C. Butcher.
Implicit Runge-Kutta Processes. 
Math. Comput. 18: 50–64, 1964.

butcher_oc2
J. C. Butcher.
Numerical Methods for Ordinary Differential Equations. 
Wiley, 2003.

butcher_book
J. C. Butcher.
Numerical Methods for Ordinary Differential Equations. 2nd edition.
Wiley, 2008.

butcher_bar
J. C. Butcher.
The non-existence of ten stage eighth order explicit Runge-Kutta methods.
BIT 25: 521–540, 1985. 

cauchy
A. L. Cauchy. 
Trente-Cinquiéme Leo̧n. In: Résum
é des leo̧ns données á l’Ecole
royale polytechnique sur le calcul infinitésimal. 
Imprimerie Royale, Paris 1823.

W
R. M. Corless, G. H. Gonnet, D. E. G. Hare, D. J. Jeffrey, D. E. Knuth. 
On the Lambert W Function.
Adv. Comput. Math. 5:  329–359, 1996.

dor
J. R. Dormand. 
Numerical Methods for
Differential Equations.
Boca Raton CRC Press, 1996.

FeSp
L. Ferracina, M.N. Spijker.
Strong stability of singly-diagonally-implicit
Runge-Kutta methods.
Appl. Numer. Math. 58(11): 1675–1686, 2008.

Fo
G. B. Folland. 
Advanced Calculus.
Prentice Hall, 2002.

gupta
G. K. Gupta, R. Sacks-Davis, T. E. Tischer.
A review of recent development in solving ODEs.
Comput. Surv. 17(1): 5–47, 1985.
 
HaNoWa
E. Hairer, S. P. Nørsett, G. Wanner.
Solving ordinary differential equations I.: Nonstiff problems.
Springer Series in Computational Mathematics 8,
Springer, 1993.

HaWa
E. Hairer, G. Wanner.
Solving ordinary differential equations II.: Stiff and Differential-Algebraic Problems.
Springer Series in Computational Mathematics 14,
Springer, 1996. 

nasa
C. Kennedy, M. Carpenter.
Diagonally Implicit Runge-Kutta Methods for Ordinary Differential Equations. A Review.
NASA/TM-2016-219173.
NASA Langley Research Center, 2016.

imex
P. K. Maurya, V. S. Yadav, B. Mahato, N. Ganta, M. K. Rajpoot, Y. G. Bhumkar.
New optimized implicit-explicit Runge-Kutta methods with applications to the hyperbolic conservation laws.
Journal of Computational Physics 446, 2021.

ramos
H. Ramos. 
A non-standard explicit integration scheme
for initial-value problems.
Appl. Math. Comput. 189(1): 710–718, 2007.

interval
J. A. dit Sandretto. 
Runge–Kutta Theory and Constraint Programming.
arXiv preprint No. 1804.04847, 2018.

skvortsov
L. M. Skvortsov. A fifth order implicit method for the numerical solution of differential algebraic equations. 
Comput. Math. and Math. Phys. 55: 962–968, 2015.

MT
K. Tvrdá, M. Minárov a.
Computation of definite integral over repeated integral.
Tatra Mt. Math. Publ. 72: 141–154, 2018.
 
NT
K. Tvrdá, P. Novotný. 
Modifications of Newton-Cotes formulas for computation of repeated integrals and derivatives. 
Czech. Math. J. 73(148): 1175–1188, 2023.

col_weight
J. Urevc, B. Starman, A. Maček, M. Halilovič. 
A novel class of collocation methods based on the weighted integral form of ODEs.
Comp. Appl. Math. 40(135), 2021.

7-10
Y. Ying, N. Yaacob.
Implicit 7-stage Tenth Order Runge-Kutta Methods Based on Gauss-Kronrod-Lobatto Quadrature Formula. 
MATEMATIKA: MJIAM  31(1): 93–109, 2015.

ZiChi
D. W. Zingg, T.T. Chisholm. 
Runge-Kutta methods for linear ordinary
differential equations.
Appl. Numer. Math. 31: 227–238,  1999.

Zl
Z. Zlatev.
Modified Diagonally Implicit Runge–Kutta Methods.
SIAM J. Sci. Stat. Comp. 2(3): 321–334, 1981.



 










§ APPENDIX

We provide the list of newly derived  as well as of some well-known implicit Runge-Kutta methods. The tables have the following columns: name/type of method, Butcher's tableau, theoretical order of convergence and values of p, q, r corresponding to simplifying order conditions (<ref>)  satisfied by the underlying methods.  




 §.§ List of new IRK  based on closed Newton-Cotes formulas


Table <ref> contains new IRK  derived by modified closed Newton-Cotes formulas, cf. Subsection <ref> for derivation of nIRKs and nIRKsc. 
There is an additional column with the corresponding stability functions R(z) in Table <ref>.


   


 §.§ List of new IRK based on open Newton-Cotes formulas


Table <ref> contains new IRK derived by modified open Newton-Cotes formulas, cf. Subsection <ref> for derivation of nIRKso and nIRKsoc. 



 §.§ List of Lobatto IIIA methods


Table <ref> contains Lobatto IIIA  methods with s=2,3,4,5 stages. The first column indicates which of the proposed methods to derive coefficients result in the well-known Lobatto IIIA for given s. See Subsection <ref>. 



 §.§ List of Gauss-Legendre and Radau methods


Table <ref> contains Radau and Gauss-Legendre methods with s=2,3 stages. Again, the first column contains the well-known name and the acronym of newly proposed approach which results in the same Butcher's tableau. See Subsection <ref>.



 §.§ List of standard collocation IRK


Table <ref>, for completeness, contains the list of standard collocation IRK which were used in our numerical experiments. For the derivation of these methods see Subsubsection <ref>. 



   



  












