Comparison of Numerical Solvers for Diﬀerential Equationsfor Holonomic Gradient Method in Statistics
N.Takayama, T.Yaguchi, Y.Zhang
November 23, 2021

1

Introduction

Definite integrals with parameters of holonomic functions satisfy holonomic systems of linear
partial diﬀerential equations. When we restrict parameters to a one dimensional curve, then thesystem becomes a linear ordinary diﬀerential equation (ODE). For example, the definite integral
with a parameter t∫+∞exp(ty − y 3 )dy

Z(t) =
0dsatisfies the ODE (3∂t2 −t)Z(t) = 1 where ∂t = dt(see [13, Chap 6] as an introductory exposition).
We can evaluate the integral by solving the linear ODE numerically. The method of solving thelinear ODE is not a high precision method, but it gives values of the integral at several evaluation
points at once. This approach to evaluate numerically definite integrals is called the holonomic
gradient method (HGM) and it is useful to evaluate several normalizing constants in statistics. Seereferences in [12] or in ,e.g., [11], [26].
The HGM consists of 3 steps. The first step is to find a linear ODE for the definite integral
Z(t) with a parameter t. The 2nd step is to find value of Z(t) and its derivatives at some pointsof t. The 3rd step is to solve the ODE. We can apply several results in computer algebra and in
the theory of special functions of several variables for the first and the 2nd step. The HGM is
successful for a broad class of distributions and the 3rd step is relatively easy in most of them.
When the function Z(t) is the dominant solution (the maximal growth rate solution) for the
ODE when t → +∞ as in the case [11, Theorem 2], the 3rd step can be performed by a simple
application of standard numerical solver for ODEs. However, there are some HGM problems forwhich we encounter diﬃculties in the 3rd step. The diﬃculty is that when the function Z(t) is not
dominant and we solve the linear ODE in a long interval, the numerical solution might become afalse one. We have this diﬃculty in our analysis of the outage probability of MIMO WiFi systems(see the section 5 and [6]) and in solving an ODE numerically of the expectation of the Euler
characteristic of a random manifold (see the section 6 and [26]). In principle, if we give an initialvalue vector exactly and the time step is suﬃciently small, we do not get a false solution. However,these assumptions are not always satisfied in practice.
In this paper, we will discuss and compare methods to solve linear ODE’s in the 3rd step of
the HGM. We intend to provide references to the HGM users for choosing numerical solvers of
ODE’s for the HGM. We will present methods that work and methods that does not work. Of
course, it depends on a given problem.

1

2

Methods to test

We will compare some methods mainly from a point of view of using for the 3rd step of the HGM1 .
In the HGM, the function Z(t) to evaluate has an integral representation and satisfies an ODE
Lf = b,L=r∑ck (t)∂tk , ∂t =

k=0

d
dt

(1)

where b(t) is an inhomogeneous term and f (t) is an unknown function. See, e.g., [6], [11], andpapers in [12] as to examples. The diﬀerential equation (1) might have a solution u which dominates
Z at +∞. In other words, there might be a solution u such that |u/Z| → ∞ when t → +∞. Insuch cases, it is not easy to obtain numerical values of Z(t) globally by solving the ODE evenwhen the almost accurate initial values of Z(t), Z ′ (t), . . . are known. See Examples 1 and 2. We
call such case the instable HGM problem. We want to obtain correct values of Z as large intervalof t as possible.
Here is a list of methods we try.
1. The Runge-Kutta method (see, e.g., [10, Chapter 2]).
2. The implicit Runge-Kutta method (see, e.g., [10, Chapters 2 and 4]).
3. Multi-step methods (see, e.g., [10, Chapters 3 and 5]).
4. The spectral method in the approximation theory (see, e.g., [2], [19], [27]).
5. Defusing method: restricting the initial value vector to a sub linear space (see, e.g., Section
3).
6. Sparse interpolations/extrapolations by ODE: solving a generalized boundary value problem
by solving a linear equation or by solving an optimization problem (see, e.g., [3], [8] and theSection 4).
We compare these methods by diﬀerent implementations on diﬀerent languages: [2], [7], [9], [16],[17], [22], [23], [24], [25]. The robustness for input errors is an important point in our comparison,because initial values or boundary values might be evaluated by Monte-Carlo integrators.
The first 4 approaches are standard. In the HGM, we may be able to evaluate the targetfunction Z(t) at any point by numerical integration methods. Of course, the value obtained is notalways accurate. For example, the accuracy of Monte Carlo integrators are not high. The last
two methods of the defusing method and sparse interpolation/extrapolation method by ODE are
expected to be used mainly under this situation of the HGM.
The program codes of this paper are obtainable from http://www.math.kobe-u.ac.jp/OpenXM/
Math/defusing/ref.html.

3

Defusing method

In the HGM, we do not always have an exact initial value vector. We often have the case thatthe target function (the target integral) have a moderate growth, but the ODE has rapid growthsolutions, too. Then, we need to correct the initial value vector so that the rapid growth solution
causes a blow-up of the numerical solution by an error of the initial value vector. We will proposea heuristic method, which we will call defusing method, to avoid a blow-up of a solution under thissituation. Any solution of the ODE can be expressed as a linear combination of basis functionsof diﬀerent growth order. The defusing method removes some basis functions of specified growth
1 Programs of this paper can be obtainable fromhttp://www.math.kobe-u.ac.jp/OpenXM/Math/defusing/ref.html. The folder/program names are given in blue
letters. (Ignore the folder names if the folder does not exist.) The letter “-” in the file name may be “ ” for somecases.

2

order to obtain the correct solution. The defusing method is an analogy of filtering some specified
frequencies from wave data. This method may be well-known, but we do not find a relevant
literature and we describe this method in this section.
In the HGM, we want to evaluate numerically a definite integral with a parameter t on someinterval of t. We put this function f (t). The HGM gives an algebraic method to find an ordinary
diﬀerential equation (ODE) for F (t) = (f (t), f ′ (t), f ′′ (t), . . . , f (r−1) )T wher r is the rank of the
ODE. We evaluate the integral by Monte-Carlo integration methods or by the saddle point approximation for large t’s. Since F (t) satisfies the diﬀerential equation, we can apply algorithms(see, e.g., [4], [14] ) to find a basis of series solutions at t = ∞. Let F∗j , j = 1, . . . , r be the basisof series solutions around t = ∞ of the ODE. We denote by Fij is the ith component of F∗j and
Fi is the ith component of F . Any solution of the ODE is written as a linear combination of F∗j ’s
over C. We assume
F11 > F12 > · · · > F1m > · · · > F1r > 0
with diﬀerent growth orders when t → +∞. We suppose thatnumerical approximation of F1 (t)→ 0,
F1j (t)t → +∞(2)

for 1 ≤ j < m. See (25) as an example. Under these assumptions, we can conclude that F (t) isapproximated as a linear combination of F∗m , . . . , F∗r . The coeﬃcients of the linear combination
can be estimated by numerical approximation values of f (t) for large t’s. Thus, we can interpolate
and extrapolate numerical values of F (t) around t = +∞ by a set of fundamental solutions of
the ODE. We call this method the defusing method by series solutions around t = +∞. See the
Example 7.
We want to find a numerical solution of the initial value problem of the ordinary diﬀerential
equation (ODE)dF
dt
F (t0 )= P (t)F
F0true ∈ Rr

=

(3)
(4)

where P (t) is an r × r matrix, F (t) is a column vector function of size r, and F0true is the initialvalue of F at t = t0 . Let us explain a defusing method for the initial value problem. Solving thisproblem is the final step of the holonomic gradient method (HGM) [12]. We often encounter thefollowing situation in the final step.
Situation 1
F0 .

1. The initial value has at most 3 digits of accuracy. We denote this initial value

2. The property |F | → 0 when t → +∞ is known, e.g., from a background of the statistics.
3. There exists a solution F̃ of (3) such that |F̃ | → +∞ or non-zero finite value when t → +∞.
Under this situation, the HGM works only for a very short interval of t because the error of theinitial value vector makes the fake solution F̃ dominant and it hides the true solution F (t). We
call this bad behavior of the HGM the instability of the HGM .
Example 1−1dF = 0dt
0
1 0−1 1  F
0 0

The solution space is spanned by F 1 = (exp(−t), 0, 0)T , F 2 = (0, exp(−t), 0)T , F 3 = (1, 1, 1)T .
The initial value (1, 0, 0)T at t = 0 yields the solution F1 . Add some errors (1, 10−30 , 10−30 )T to
the initial value. Then, we have

3

t
50
60
70
80

value F1 by RK
1.92827e-22
8.75556e-27
1.39737e-30
1.00002e-30

diﬀerence F1 − F11
9.99959e-31
1.00000e-30
1.00000e-30
1.00000e-30

We can see the instability.
Example 2(
P (t) =

0 1
t 0

)
.

This diﬀerential equation is obtained from the Airy diﬀerential equationy ′′ (t) − ty(t) = 0
by putting F = (y(t), y ′ (t))T . It is well-known that the Airy function
( 3
)
∫ b
s
1
lim
cos
+ ts ds
Ai(t) =
π b→+∞ 0
3
is a solution of the Airy diﬀerential equation and1= 0.355028053887817 · · ·32/3 Γ(2/3)
1Ai′ (0) == −0.258819403792807 · · ·31/3 Γ(1/3)
lim Ai(t) = 0
Ai(0) =

t→+∞lim Ai′ (t)

t→+∞

= 0

Figure 1 is a graph of Airy Ai function and Airy Bi function. The function F (t) = (Ai(t), Ai′ (t))Tsatisfies the condition 2 of the Situation 1 of the instability problem.

Figure 1: Airy Ai(x) and Bi(x) drawn by Mathematica
We can also see that the condition 3 of the Situation 1 holds by applying the theory of singularity of ordinary diﬀerential equations (see, e.g., the manual DEtools/formal sol of Maple [15]and its references on the theory, which has a long history). In fact, the general solution of the
Airy diﬀerential equation is expressed as
)
(
2 3/2
−1/4
(1 + O(t−3/2 ))C1 t
exp − t3()
2 3/2
+ C2 t−1/4 expt(1 + O(t−3/2 ))3
4

when t → +∞ where Ci ’s are arbitrary constants.
We note that the high precision evaluation of the Airy function is studied by several methods(see, e.g., [5] and its references). Some mathematical software systems have evaluation functionsof the Airy function. For example, N[AiryAi[10]] gives the value of Ai(10) on Mathematica. By
utilizing these advanced evaluation methods, we use the Airy diﬀerential equation for our test caseto check the validity of our heuristic algorithm.
We are going to propose some heuristic methods to avoid the instability problem of the HGM.
Numerical schemes such as the Runge-Kutta method obtain a numerical solution by the recurrence
Fk+1 = Q(k, h)Fk

(5)from F0 where Q(k, h)2 is an r × r matrix determined by a numerical scheme and h is a small
number The vector Fk is an approximate value of F (t) at t = tk = t0 + hk.
Example 3 The Euler method assumes dF/dt(t) is approximated by (F (t + h) − F (t))/h and thescheme of this method is
Fk+1 = (E + hP (tk ))Fk
where E is the r × r identity matrix.
In case that the initial value vector F0 contains an error, the error may generate a blow-up
solution F̃ under the Situation 1 and we cannot obtain the true solution.
Let N be a suitable natural number and put
Q = Q(N − 1, h)Q(N − 2, h) · · · Q(1, h)Q(0, h)(6)

We call Q the matrix factorial of Q(k, h). The matrix Q approximates the fundamental solution
matrix of the ODE. We assume the eigenvalues of Q are positive real and distinct to simplify ourpresentation. The following heuristic algorithm avoids to get the blow-up solution.
Algorithm 1
1. Obtain eigenvalues λ1 > λ2 > · · · > λr > 0 of Q and the corresponding
eigenvectors v1 , . . . , vr .
2. Let λm be the eigenvalue which is almost equal to 0.
3. Express the initial value vector F0 containing errors in terms of vi ’s as
F 0 = f 1 v1 + · · · + f r vr ,fi ∈ R(7)

4. Choose a constant c such that F0′ := c(fm vm + · · · + fr vr ) approximates F0 .
5. Determine FN by FN = QF0′ with the new initial value vector F0′ .
We call this algorithm the defusing method for initial value problem. This is a heuristic algorithm
and the vector F0′ gives a better approximation of the initial value vector than F0 in several
examples and we can avoid the blow-up of the numerical solution.
Example 4 We set t0 = 0, h = 10−3 , N = 10 × 103 and use the 4-th order Runge-Kutta
scheme. We have λ1 = 9.708 × 109 , v1 = (−5.097, −159.919)T and λ2 = 3.247 × 10−7 , v2 =(−5.09798, 37.164813649680576037539971418209465086)T = (a, b). Then, m = 2. We assume the
3 digits accuracy of the value Ai(0) as 0.355 and set F0′ = (0.355, 0.355b/a). Then, the obtained
value F5000 at t = 5 is (0.000108088745179140, −0.000246853220440734). We have the following
accurate value by Mathematica
2 It was denoted by Q(t + kh, h) in the previous section. We denote Q(t + kh, h) by Q(k, h) as long as no00
confusion arises.

5

In[1]:= N[AiryAi[5]]Out[1]= 0.000108344
In[2]:= N[D[AiryAi[x],{x}] /. {x->5}]
Out[2]= -0.000247414
Note that 3 digits accuracy has been kept for the value Ai(5). On the other hand, we appy the
4th order Runge-Kutta method with h = 10−3 for F0 = (0.355, −0.259)T , which has the accuracyof 3 digits. It gives the value at t = 5 as (−0.147395, −0.322215), which is a completely wrong
value, and the value at t = 10 as (−102173, −320491), which is a blow-up solution.
This heuristic algorithm avoids the blow-up of the numerical solution. Moreover, when the
numerical scheme gives a good approximate solution for the exact initial value, we can give anerror estimate of the solution by our algorithm. Let | · | be the Eucledian norm.
Lemma 1 Let F (t) be the solution. When |QF0true − F (N h)| < δ holds, we have
|QF0′ − F (N h)| < |QF0′ | + |F (N h)| + 2δ(8)for any F0′ ∈ Rn .
Proof . It is a consequence of the triangular inequality. In fact, we have
|QF0′ − F (N h)|= |QF0′ − QF0true + QF0true − F (N h)|≤ |QF0′ − QF0true | + |QF0true − F (N h)|≤ |QF0′ | + |QF0true | + δ≤ |QF0′ | + |F (N h)| + 2δ//Under the Situation 1, |F (N h)| is small enough. Then, it follows from the Lemma that |QF0′ |should be small. In this context, we can give an error estimate of our algorithm. However, ournumerical experiments present that the algorithm shows a better behavior than this theoreticalerror estimate. Then, we would like to classify our defusing method as a heuristic method.

4

Sparse interpolation/extrapolation method

We solve the ODE (1) of rank r when (approximate) values of f (t) = Z(t) at t = p1 , p2 , . . . , pr are
known. We call the points (pi , qi ), qi = f (pi ) data points. In other words, we want to interpolate
or extrapolate values of f from these r values. This is a generalization of the boundary value
problems. The number of data points may be more than the rank or less than the rank of the
ODE in the methods B and C in this section. A standard numerical method to find values off on an interval [ts , te ] is as follows. Devide [ts , te ] into N intervals. Let ti be ts + hi whereh = (te − ts )/N . We assume that the set of all pj ’s are a subset of the set {ti }. We denote bykfi the value of f (ti ). We introduce the backward shift operator ∇fi = fi − fi−1 . Then ∇hkfi isapproximately equal to f (k) (ti ) or f (k) (ti−1 ) or . . . or f (k) (ti−k ). For example, when k = 1, we havef (ti ) − f (ti−1 ) = f ′ (ti )h + O(h2 ) if we make the Taylor expansion of f (ti−1 ) = f (ti − h) at t = ti
and f (ti ) − f (ti−1 ) = f ′ (ti−1 )h + O(h2 ) if we make the Taylor expansion of f (ti ) = f (ti−1 + h) at
t = ti−1 . We haver∑∇k fi+sk∼ b(ti ), 0 ≤ sk ≤ k.
(9)
ck (ti )hk
k=0Here si is an integer to choose an approximation of f (k) (ti ). By assuming the left hand side andthe right hand side are equal and giving values fj for t = pj , 1 ≤ j ≤ r, we have a system of linear
6equations of the formf0
 f1 


 · 

=B
A

 · 
 · 
fN

(10)

where A is a (N +1)×(N +1) matrix and B is a column vector of length N +1. Solving this equation,
we obtain approximate values of f (ti ). We call this method the sparse interpolation/extrapolation
method A for HGM.
Example 5 We solve Lf = 0, (intro/2021-06-11-sparse-interp.ipynb) L = (∂t − 1)(∂t2 − t) =∂t3 − ∂22 − t∂t + t − 1 on t ∈ [−9, 0] with f (−9) = Ai(−9) ∼ −0.0221337215473414, f (−4) =Ai(−4) ∼ −0.0702655329492895, f (0) = Ai(0) ∼ 0.355028053887817 and N = 100. Then, we
obtain approximate values of the Airy function Ai(t). The numerical result shows that we do nothave a false solution in spite of the factor ∂t − 1.
An alternative method to solve (10) when b = 0 is to construct r linearly indedent solutions(fundamental solutions) of (9) and to find coeﬃcients to express the solution of the generalizedboundary conditions as the linear combination of the fundamental solutions.
Let us introduce the sparse interpolation/extrapolation method B for HGM. This is also astandard method and using it in the HGM is very useful. A numerical integration for a function
g can be expressed asN∑
IN (g) =Tj g(tj )(11)
j=0

where t0 = t2 < t1 < · · · < tN −1 < tN = te and Tj ∈ R≥0 . We fix a numerical integration−ts
method. For example, the trapezoidal method can be expressed as h = teN, ti = ts + hi, Tj = h
for 1 ≤ j < N and T0 = TN = h/2. We approximately expand the solution f by a given basis
functions {ek (t)}, k = 0, 1, . . . , M asf (t) =

M∑

fk ek (t),fk ∈ R.

(12)

k=0

We put M + 1 = m to be compatible with the variable name in our program. Put this∫ expressiontinto Lf = b. We minimize the loss function, which is the numerical integration value of tse |Lf (t)−b(t)|2 dtN∑ℓ({fk }) =|(Lf )(tj ) − b(tj )|2 Tj

(13)

j=0
N
M
∑
√
√ ∑
fk (Lek )(tj ) − Tj b(tj )
=
Tj
j=0

2

k=0

under the constraints at data points
M∑

fk · ek (pi ) = qi ,i = 1, 2, . . . , r.

(14)

k=0

This is a least mean square problem under constraints (see, e.g., [18]). Since the loss function (13)is defined by the numerical integration formula (11), we have the following estimate.

7∫tLemma 2 The norm tse |Lf (t)−b(t)|2 dt is bounded by ℓ({fk })+eN (|Lf −b|2 ) where eN (|Lf −b|2 )is the error of the numerical integration method.
Solvers of the least mean square problem output the value of the loss function, then we canestimate the L2 norm of Lf − b by this Lemma.
We can also consider a least mean square problem with no constraints by the loss functionℓ̃({fk }) = αℓ({fk }) + β

(M
r
∑
∑
i=1

)2
fk · ek (pi ) − qi

+γ

N
∑

fi2 .

(15)

i=0

k=0

Here α, β, γ are hyperparameters for the optimization. In particular, γ > 0 is used to avoid an
overfitting.
Example 6 We solve the diﬀerential equation of Example 5. The method B does not work
well on the interval [−9, 0] because of the oscillation of the solution, but it works on a smaller
interval [−4, 0]. We use Chebyshef polynomials on [−4, 0] upto degree 9 as basis functions
and h = 0.01. Data points are [[−4, −0.0702655329492895], [−3, −0.37881429], [−2, 0.22740743]].
(intro/sib-intro.py) (asir-tmp/sib-intro.rr, trybintro3())
This method works well for some problems which are hard by standard methods. See the
Example 12 as to an application of this method to a 4th order ODE. An application of thismethod to a problem on random matrices is given in Example 14.
Let us introduce the sparse interpolation/extrapolation method C for HGM. We use the basis
functions ej (t) as in the method C. We construct an m × m symmetric matrix S of which (i, j)element is∫te(Lei (t))(Lej (t))dt

(16)

ts

Let F be the column vector of length m = M + 1: F = (f0 , f1 , . . . , fM )T . Then, we have
∫ te
ts


L

M
∑

2
fj ej (t) dt = F T SF.

(17)

j=0

∑M
It follows from this idenity that when j=0 fj ej (t) is an approximate solution of Lf = 0, we
may expect the value of the quadratic form F T SF is small. Let (pk , qk )’s be given data pointsof (t, f (t)). The method C finds the vector F by solving the quadratic programming problem ofminimizing
F T SF +

r
∑




k=1

= F T SF +

r
∑
(

M
∑

2
fj ej (pk ) − qk 

j=0

F T Sk2 F − 2qk Sk F + qk2

)

k=1
T

= F T SF + F (PeT Pe )F − 2PeT Q + QT Q
where Sk = diag(e0 (pk ), . . . , eM (pk )) ande0 (p1 ) · · · eM (p1 ) e0 (p2 ) · · · eM (p2 ) 
Pe = ,·····e0 (pr ) · · · eM (pr )8Q = (q1 , . . . , qr )T

(18)
(19)

(20)

Or, we may solve the quadratic programming problem of minimizing F T SF under the constraintsM∑

fj ej (pk ) − qk = 0,

k = 1, . . . , r.

(21)

j=0

The method C can also solve the problem of Example 6. (intro/sic-intro.py) (asir-tmp/sic-intro.rr,
trycintro3())

5

Tests of the methods to the function Hnk

Let n and k be positive integers. We define the function Hnk (x, y) by∫ x
Hnk (x, y) =tk exp(−t)0 F1 (; n; yt)dt
0
∫
√
Γ(n)
tk (1 − s2 )n−3/2 exp(−t − 2s yt)dtds= √πΓ(n − 1/2) D(x)

(22)
(23)

where D(x) = {(t, s) ∈ [0, x] × [−1, 1]}
This function appears in studies of the outage probability of MIMO WiFi systems (see, e.g., [6]).
We will compare the methods in the section 2 for this function on several systems. The function
satisfies the following system of linear partial diﬀerential equations.
Proposition 1 ([6])
1. The function u = Hnk (x, y) satisfies{θy (θy + n − 1) + y(θx − θy − k − 1)} • u= 0,(θx − θy − k − 1 + x) θx • u

= 0.

∂
∂
where θx = x ∂x,θy = y ∂y
. The holonomic rank of this system is 4.

2. The function u is a solution of the following ODE with respect to y.
y 2 ∂y4 +(−y +2n+2)y∂y3 +(−yx+(−k−n−3)y +n(n+1))∂y2 +((y −n)x−n(k+2))∂y +(k+1)x(24)(Prog paper/19-a19-n-pf.rr)
When y → +∞, solutions of the system has the following asymptotic behavior. It is shown by
the DEtools[formal sol] function of Maple [15]. (Hkn/a4-f2.ml)h1h2= (xy)−1/2(1/2+n) exp(−2(xy)1/2 )(1 + O(1/y 1/2 )),= y −k−1 (1 + O(1/y)),h3h4= (xy)−1/2(1/2+n) exp(2(xy)1/2 )(1 + O(1/y 1/2 )),= y 1−n+k exp(y)(1 + O(1/y)),Which is the asymptotic behavior of the function Hnk (x, y) when x is fixed? We compare the valueof h4 and the value by a numerical integration in Mathematica3 . The integration by Mathematica
is done as follows.
--> hh[k_,n_,x_,y_]:=NIntegrate[t^k*Exp[-t]*HypergeometricPFQ[{},{n},t*y],{t,0,x}];--> hh[10,1,1,1000]
3 The quality of HypergeometricPFQ of mathematica is extremely high. However, the method to evaluate hypergeometric functions in Mathematica is still a black box. It is not easy to give a numerical evaluator of hypergeometricfunctions which matches to Mathematica in all ranges of parameters and independent variables.

9

y
1000
2000
3000

Ratio
7.36595030875893e-452
2.64621603289928e-881
2.67723893601667e-1311

whereRatio = (H110 (1/2, y))/(y 1−n+k exp(y)).

(25)

This computational experiments suggest that Hnk is expressed by h1 , h2 , h3 without the dominantcomponent h4 as explained in (2).
Example 7 We approximate H110 (1, y) by h3 of the asymptotic series terms truncated by O(y −3 ).
(Hkn/2021-06-06-asymp.rr) Let ys be a suﬃciently large number. We determine the constant Cs
h (y)−H110 (1,y)
by H110 (1, ys )/h3 (ys ). Then we have the following data of relative errors re (y) = 3 H 10 (1,y)
.
1

ys
102
104

re (ys + 10)
0.0763
5.7 × 10−10re (ys + 90)
0.259
5.09 × 10−9

The data tells that the defusing by series around y = +∞ works well near y = 104 , but we need
to use other methods near y = 102 . The sparse interpolation/extrapolation method B in the
Example 12 is a more refined variation of this method.
Let us apply methods presented in the Section 2 to evaluate H110 (1, y).
Example 8 We apply the implicit Runge-Kutta method to H110 (1, y) which is a solution of the
ODE (24). We translate the ODE into a system of ODE for U by U = (u, u′ , u′′ , u′′′ )T and derive
an ODE for F = U exp(−y)y −(1−n+k) This is the Gauge transformation by the exponential partof h4 so that the solutions keep bounded values. Then, the column vector function F satisfies
F ′ = P F where −y−k+n−1

1
0
0
y
−y−k+n−1


1
0
0


y
P =
(26)
.
−y−k+n−1
0
0
1


y
(−k−1)xy2(−y+n)x+nk+2ny2

yx+(k+n+3)y−n2 −ny2−k−n−3y

The following data is obtained by the implicit Runge-Kutta method (Gauss method of order 10,[10, IV.5]) implemented in Mathematica [16]. The Figure 2 presents the value of log H110 (1, y) (thealmost straight graph), the value by the implicit Runge-Kutta method, and the value by NDSolvewith the ExplicitRungeKutta option (the graph of yellow green). (Hkn/2021-05-17-i-RK.nb) The
initial value vector is evaluated at y = 1 by the numerical integration with the accuracy about
16 digits and the initial step size is 10−3 . The implicit Runge-Kutta method works upto abouty = 25. This interval [1, 25] is larger than the valid interval [1, 4.5] of the ExplicitRungeKutta
method of NDSolve. Why the implicit method works well? By virtue of the Gauge transformation,
the maximal eigenvalue of the constant part of the coeﬃcient matrix is close to 0. Then, it seems
that we have relatively larger valid interval by the implicit Runge-Kutta method (see, e.g., [10,IV.3]).
We note that the implicit Runge-Kutta method is used in a large scale computation withthe spectral diferred correction (SDC), see, e.g., [28], [20]. It is a method to accelerate the implicit Runge-Kutta method by making a preconditioning to iteration schemes for solving algebraic
equations of the implicit method. We do not need this method for our relatively small rank
ODE’s in this paper. However, it will be useful for solving high rank ODE’s in HGM. (numht3/pySDC/airy playground.py, Airy implicit.py)
When we can obtain exact initial values by the numerical integration (e.g., of (22)), it will bea good strategy to extrapolate the values in a short interval by solving the ODE by the implicit
10

60
50
40
30
20
10

5

10

15

20

25

30

35

-10

Figure 2: The Runge-Kutta(yellow green), the implicit Runge-Kutta method(blue), the exact
value(orcher)
18
105

16

110

115

120

125

130

-1. × 10-6

-2. × 10-6

14

-3. × 10-6

12
-4. × 10-6

110

120

130

140

150

Figure 3: The implicit Runge-Kutta method vs the exact value. The right graph is the relative
error (Hr − H)/H where H is the exact value and Hr is the value by the implicit Runge-Kutta
method.
Runge-Kutta method. The Figure 3 describes the value of log H110 (1, y) and the value by the
implicit Runge-Kutta method on the interval [100, 150]. The initial value vector is evaluated at
y = 100 by the numerical integration with the accuracy about 16 digits and the inital step size is10−3 . It works upto about y = 130.
Example 9 The spectral method in the approximation theory is a powerful method and we applied the ApproxFun package [2] implemented in Julia for solving the ODE (24). The outputs are
remarkable. We solve the ODE on [1, 40] with the initial value at y = 1 as (u, u′ , u′′ , u′′′ )(1) =(0.07810139136088563, 0.05096276584900834, 0.02050273784371611, 0.005887855153702640) whereu(y) = H110 (1, y). (Julia/hkn4i.jl) Then, the output value at y = 40 is 815.0105773595695,
which agrees with the value of numerical integration by Mathematica 815.0105773587113 upto
11 digits. It is also powerful to solve boundary value problems. We give the boundary values(u(1), u(1.1), u(39.9), u(40)) = (0.0781014, 0.0833012, 803.121, 815.011), which is less accurate thanthe previous example. (Julia/hkn4b.jl) Then, the output value at y = 20 is 27.021711397385513,
which agrees with the value of numerical integration by Mathematica 27.021701160033859079 upto
6 digits. Note that we only give the boundary values in 6 digits accuracy. Although it works remarkably with a proper setting, we have had some troubles. For example, when we input the
ODE (24)×y 2 , the program returns NaN. Another example is that it did not work for a first ordersystem obtained from (24). (Julia/hkn3b.jl)
Example 10 We can try several solvers of ODE by solve ivp on scipy/python [23]. (Hkn/202106-02-hkn.ipynb) We solve the initial value problem of (26) on y ∈ [1, 30] with the absolute
tolerance 0 and relative tolerance 10−3 4 . Note that the default values of the absolute and relative
tolerances does not work well. Here is a table of y’s for which the relative error becomes morethan 0.3. Refer to [23] on details on the following methods. (Hkn/2021-06-08-hkn-rerr2.ipynb)
4 atol=0,

rtol=1e-3

11Figure 4: Solving by the sparse interpolation
A

Figure 5: Solving the boundary value problem with random errors

Method
RK45
Radau(Implicit Runge-Kutta)
BDF
LSODA(Adams/BDF on ODEPACK)

Failure y (relative error > 0.3)
21.5
26.1
16.8
20.0

We can also try the deSolve package in R [7]. Some codes are common (e.g., FOTRAN package
ODEPACK) with the python package and results are analogous with above. (Hkn/hkn.r)
Example 11 We apply the sparse interpolation/extrapolation method A. In other words, wesolve a generalized boundary value problem of the ODE (24) on y ∈ [104 , 104 + 40]. (Hkn/hkn10-1-10000si-random.py) We give generalized boundary values of u(y) = Hnk (1, y) × 10−80 at
p1 = 104 , p2 = 104 + 1333h, p3 = 104 + 2000h, p4 = 104 + 4000h = 104 + 40 where h = 0.01. We
approximate derivatives asu(1) (y)=u(2) (y)=u(3) (y)=u(4) (y)=1(u(y) − u(y − h))h
1(u(y + h) − 2u(y) + u(y − h))h2
1(u(y + h) − 3u(y) + 3u(y − h) − u(y − 2h))h3
1(u(y + 2h) − 4u(y + h) + 6u(y) − 4u(y − h) + u(y − 2h))
h4

(27)

and solve the linear equation (10). The linear equation is solved by the function linsolv in the
scipy package [22]. The condition number of the matrix A becomes about 2.6 × 1013 . We give
random errors by assuming the significant digits of the given boundary values are 3. The Figure 4
is a graph of the solutions of 30 tries of random errors. It works well. We note that the robustness
with respect to random errors depends on a choice of pi ’s. (Hkn/hkn-10-1-10000bv-random.py)
For example, we change p2 to p1 + h and p3 to p4 − h. This stands for solving the boundary value
problem with the boundary values of f and f ′ at the boundary p1 and p4 . We solve the linear
equation with random errors. We can see that errors are magnified in middle from the Figure 5.
Although the boundary value problem with exact boundary values of Hnk (1, y) and its first
derivatives can be solved nicely on [1, 40], application of this method on [1, 40] with internalreference points is not good. (Hkn/hkn-10-1-1-40-si.py, Hkn/2021-06-11-sparse-interp.ipynb) We
take the four points as 1, 13.99675, 20.5, 40 and h = 0.00975 and apply the method. It has an
acceptable relative error 33.23 in the interval [1, 6] and the relative error becomes acceptablesmaller errors out of this interval.
(Hkn/2021-06-09-riccati-Hkn-rtol.ipynb) A clever method to solve boundary value problems is
to use the Riccati transformation [8, p.1059] (see also [1, p.149]). We implemented this methodwith solve ivp on scipy/python. It works for a small interval, e.g., y ∈ [1, 2]. But, it fails
12

on y ∈ [1, 10]. The numerical solution of the associated Riccati equation increases from 0 to
−1 × 1065 by the Runge-Kutta method (method=RK4 with rtol=1e-13, atol=1e-10) and the
backward equation cannot be solved because of the overflow. We tried on a smaller interval [1, 3]
with the Adams/BDF method (method=LSODA). Evaluation of the Riccati equation did not stopin 8 minutes on the google colaboratory.
Example 12 We apply the sparse interpolation/extrapolation method B explained in Section 4.
By utilizing the local solution expansion at the infinity of the ODE (24), we use the set of the 4functionsej (t) = t−3/4 exp(2t1/2 )t−j/2 , j = 0, 1, 2, 3, t = y(28)
as the basis for (12). We will approximate the solution on [ts , te ]. The pi ’s in the constraint (14)are
p0 = ts , p1 = ts + 5, p2 = ts + 10, . . . , pk = te − 1.
We do not use the constraint and add the diﬀerence of the approximate value and the true value
to the the loss function ask∑()2(approximate value at pj ) − H110 (1, pj )j=0

We expect that this basis will give a good approximation when ts and te are large. In fact, when[ts , te ] = [1, 40], the maximum of the relative error is 5.47 by our implementation on least squaresin scipy/python. On the other hand, this basis gives a good approximation on [20, 20 + 40] and[104 , 104 + 40]. We give random errors to the value of qj = H110 (1, pj ) as qj × (1 + O(10−3 )). It isalso robust to these random errors. (Hkn/hkn-10-1-10000sib-random.py Hkn/hkn-10-1-20-60sibrandom.py)[ts , te ]max of relative errors with exact qj
max of relative errors with 30 random errors[20, 60]
6.21 × 10−3
1.39 × 10−2[104 , 104 + 40]
2.67 × 10−12
4.07 × 10−3

In summary, we should use the implicit Runge-Kutta method or the sparse interpolation/extrapolation
method A when ts and te are small and use the sparse interpolation/extrapolation method B when
they become large for the problem Hnk . An alternative choice when ts and te are relatively small
will be the following defusing method.
Example 13 We implement the defusing method in tk ode by mpfr.rr 5 for the Risa/Asir [21].
It generates C codes utilizing the MPFR [17] for bigfloat and the GSL [9] for eigenvalues andeigenvectors. We apply the defusing method for initial value problem to H110 (1, y) which is asolution of the ODE (24). We use the step size h = 10−3 and the bigfloat of 30 digits of accuracy. (Hkn/defusing/tmp-test.c, tmp-proj.c) (Generated by asir-tmp/tk-ode-assert.rr, tk-odeassert.hkn1(), tk-ode-assert.hkn2()) The Figure 6 shows that the adaptive Runge-Kutta methodof GSL [9] fails before y becomes 30. (Hkn/defusing/a26-y.c) The Figure 7 presents the relative
error of values by the defusing method and exact values. It shows that the defusing method works
even when y = 103 .

6

Tests of some methods to the function E[χ(Mt )]We consider the integral studied in [26]E[χ(Mt )]==F (s1 , s2 , m11 , m21 , m22 ; t)∫∫ ∞∫ 2π ∫ 2π
{ 1 }
s1 s2
1 ∞
dσ
db
dθ
dϕ(σ 2 − b2 )exp − R ,
2
2 t
(2π)
2
−∞
0
0

5 http://www.math.kobe-u.ac.jp/OpenXM/Current/doc/asir-contrib/ja/tk_ode_by_mpfr-html/tk_ode_by_mpfr- ja.html. Todo, English manual.

13

(29)

Figure 6: log H110 (1, y). Exact value (by numerical integration) and the value by our defusing method agree. The adaptive RungeKutta method with the initial relative error
10−20 (upper curve) does not agree with theexact value when y is larger than about 25.

Figure 7: The relative error of H110 (1, y) of
our defusing method. The relative error isdefined as (Hd −H)/H where Hd is the value
by the defusing method and H is the exact
value.

whereR2= s1 (b sin θ sin ϕ + σ cos θ cos ϕ − m11 ) + s2 (σ sin θ cos ϕ − b cos θ sin ϕ − m21 )222+s1 (σ cos θ sin ϕ − b sin θ cos ϕ) + s2 (b cos θ cos ϕ + σ sin θ sin ϕ − m22 ) .
and
m11 = 1, m21 = 2, m22 = 3, s1 = 103 , s2 = 102 .
By virtue of the Euler characteristic method, the expectation of the Euler characteristic of arandom manifold Mt E[χ(Mt )] approximates the probability that the maximal eigenvalue of 2 × 2
random matrices is less than t. See [26, Sec 3] for details.
The function E[χ(Mt )] is a solution of the rank 11 ODE Lf = 0 discussed in [26, Example 5].
The operator L is of the form()
(30)(−4.72 × 10−52 t29 + · · ·)∂t10 + · · · + (−7.78 × 10−22 t35 + · · ·) ∂t
by multiplying a constant 10−31 /8.66 6 and ∂t from the right to the operator given in https://yzhang1616.github.io/ec1/ec1.html. (or ec/tryb6/ann3.txt, variable ODE) Note that thisODE is of a form of the singular perturbation problem and it seems to be hard to solve initialvalue problems like the Runge-Kutta method. In [26], since the exact numerical integration of(29) is not easy, they use values of the Monte-Carlo simulation on t ∈ [3.8, 3.81] to solve the ODE
by power series. This series approximates the expectation upto t = 3.8633, but it does not when t
is larger than it. See [26] for details. The sparse interpolation/extrapolation method B overcomes
this diﬃculty.
Example 14 (ec/tryb6/2021-07-09-try6b-tmpb.py) (It is generated by asir-tmp/sib-yi5b.rr, tryb6())
The sparse extrapolation method B gives an approximate solution which has larger valid domain
than the solution by the power series method. We use the basis functions
ej (t) = (t − 3.8055)j .

(31)

Let pi = 3.8 + i/1000, i = 0, 1, . . . , 9. The values qi at pi are[0.067160, 0.065485, 0.064732, 0.063315, 0.061814,
0.060477, 0.059611, 0.058257, 0.057520, 0.055971](32)

6 This constant is chosen so that the maximal absolute value of the coeﬃcients of f ’s in (13) of the sparsekinterpolation/extrapolation method B is 1 in Example 14.

14

Figure 8: The graph of F29 (t) and simulation values in the left and relative errors in the right.
The data points are marked with ’x’.

Figure 9: The graph of F29 (t)’s with random errors for data points and simulation values andrelative errors in the right.
by a Monte-Carlo simulation. (ec/tryb6/yiex5b.r) (pi , qi )’s are data points for the sparse extrapolation method B. We divide the interval [ts , te ] = [3.8, 4.0] into 200 segments and use thetrapezoidal formula for ℓ(f ). The optimal solution[f0 , f1 , . . . , f29 ]= [0.060145405402867516, −1.20074804549872, 9.438660716835336,−29.022737131194667, −46.606911486598264, 587.9594544508735, . . .]can be found in about 18 seconds 7 by least squares of the scipy package [24] on python. The
value of the loss function ℓ̃ is equal to 3.85 × 10−7 . The integral∫ te2(LF29 ) dt where F29 (t) =
ts

29
∑

fj ej (t)

(33)

j=0

is equal to 2.58 × 10−4 and |F29 (pi ) − qi | < 3.18 × 10−4 for all i. The left graph of Figure 8 is
the graph of F29 (t) and dots are approximate values of E(χ(Mt )) by the Monte-Carlo simulation.
The data points are marked with ’x’. The right graph of Figure 8 is the relative error of the valuesof F29 and the values by the simulation. Although the relative error is rather large, we disagree
that the method is useless. For example, we have F29 (3.935) = 0.00241, then we can conclude thevalue of t satisfying f (t) = 0.001 may be found in the domain t > 3.935 assuming that the relative
error is less than 1.4.
Example 15 (ec/tryb5/yi5b-random.py) We give relative errors less than 10−3 for qk of the data
points (pk , qk ). We use the same scheme as Example 14. We execute 30 tries and obtain the resultsin Figure 9. Relative errors may be more than 40 at far points from the data points. We conclude
that we should give accurate values of data points to extrapolate by the sparse extrapolation
method B for this problem.
We expected that the sparse interpolation/extrapolation method C works for this problem,
too. However, it does not seem to work well in real implementations on the double arithmetic.
7 Debian 10.2 with Intel(R) Xeon(R) CPU E5-4650 2.70GHz

15

In fact, the cvxopt package [25] in python returns ArithmeticError:5 in our implementation of
the method C. (ec/tryb6/2021-07-14-tryc6-tmpb-modifed.py) We expect that an implementationof the quadratic programming with the bigfloat will solve this problem.
Acknowledgments. The first and the second authors are supported by the JST CREST Grant
Number JP19209317. The first author is supported by the JSPS KAKENHI Grant Number
21K03270. The thrid author is supported by XJTLU Research Development Funding RDF-2001-12, NSFC Young Scientist Fund No. 12101506, and Natural Science Foundation of the Jiangsu
Higher Education Institutions of China Programme-General Programme No. 21KJB110032.

References[1] F. S. Acton, Numerical Methods that Work, 1990, Mathematical Association of America.
[2] ApproxFun.jl, Julia package for function approximation,
https://github.com/JuliaApproximation/ApproxFun.jl,
https://juliaapproximation.github.io/ApproxFun.jl/latest/.
[3] U.M.Ascher , R.M.M.Mattheij, R.D.Russel, Numerical Solution of Boundary Value Problems
for Ordinary Diﬀerential Equations, 1987, SIAM.
[4] M. Barkatou, An algorithm to Compute the Exponential Part of a Formal Fundamental
Matrix Solution of a Linear Diﬀerential System, Applicable Algebra in Engineering, Communication, and Computing 8 (1997), 1–23.
[5] S. Chevillard, M. Mezzarobba, Multiple-precision evaluation of the Airy Ai function with
reduced cancellation, arxiv:1212.4731.
[6] F.H.Danufane, K.Ohara, N.Takayama, C.Siriteanu, Holonomic Gradient Method-Based CDF
Evaluation for the Largest Eigenvalue of a Complex Noncentral Wishart Matrix.
https://arxiv.org/abs/1707.02564[7] DeSolve (R-package), Solvers for Initial Value Problems of Diﬀerential Equations,
http://desolve.r-forge.r-project.org/[8] L. Dieci, M. Osborne, R. Russell, A Riccati Transformation Method for Solving Linear BFPs
I, SIAM Journal on Numerical Analysis 25 (1988), 1055–1073.
[9] GSL, GNU scientific library,https://www.gnu.org/software/gsl/[10] E.Hailer, S.P.Norsett, G.Wanner, Solving Ordinary Diﬀerential Equations I, II. Springer,
1993.
[11] H.Hashiguchi, Y.Numata, N.Takayama, A.Takemura, Holonomic gradient method for the
distribution function of the largest root of a Wishart matrix, Journal of Multivariate Analysis,
117, (2013) 296-312.
[12] References for HGM, http://www.math.kobe-u.ac.jp/OpenXM/Math/hgm/ref-hgm.html[13] T. Hibi and et al, Gröbner Bases: Statistics and Software Systems, 2013, Springer[14] M. van Hoeij, Formal Solutions and Factorization of Diﬀerential Operators with Power Series
Coeﬃcients, Journal of Symbolic Computation 24 (1997), 1–30.
[15] DEtools in Maple,
https://www.maplesoft.com/support/help/maple/view.aspx?path=DEtools16[16] NDSolve in Mathematica,
https://reference.wolfram.com/language/tutorial/NDSolveImplicitRungeKutta.
html.
[17] The GNU MPFR library, https://www.mpfr.org/[18] J.Nocedal, S.Wright, Numerical Optimization, 2006, Springer.
[19] S. Olver, A. Townsend, A fast and well-conditioned spectral method, SIAM Review 55 (2013),
462489.
[20] pySDC project, a Python implementation of the spectral diferred correction,
https://parallel-in-time.org/pySDC/[21] Risa/Asir, a Computer algebra system, http://www.openxm.org[22] scipy.linalg.solve
https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.solve.html#scipy.linalg.solve[23] scipy.integrate.solve ivp,
https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html[24] scipy.optimize.least squares,
https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.least_squares.html[25] cvxopt.solver.qp,
https://cvxopt.org/userguide/coneprog.html#quadratic-programming[26] N.Takayama, L.Jiu, S.Kuriki, Y.Zhang, Computations of the Expected Euler Characteristic
for the Largest Eigenvalue of a Real Wishart Matrix, Journal of Multivariate Analysis 179(2020), 104642.
[27] L. N. Trefethen, Approximation Theory and Approximation Practice, 2020, SIAM.
[28] M.Winkel, R.Speck, D.Ruprecht, A high-order Boris integrator, Journal of Computational
Physics 295 (2015), 456–474.

17