# Comparison of Numerical Solvers for Differential Equations for Holonomic
Gradient Method in Statistics

N.Takayama, T.Yaguchi, Y.Zhang

(August 13, 2023)

##  1 Introduction

Definite integrals with parameters of holonomic functions satisfy holonomic
systems of linear partial differential equations. When we restrict parameters
to a one dimensional curve, then the system becomes a linear ordinary
differential equation (ODE). For example, the definite integral with a
parameter $t$

$Z(t)=\int_{0}^{+\infty}\exp(ty-y^{3})dy$

satisfies the ODE $(3\partial_{t}^{2}-t)Z(t)=1$ where
$\partial_{t}=\frac{d}{dt}$ (see [16, Chap 6] as an introductory exposition).
We can evaluate the integral by solving the linear ODE numerically. The method
of solving the linear 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. See references in [15] or in ,e.g., [14], [30].

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 points of $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 2nd steps. Results in
numerical integration can be applied to 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\rightarrow+\infty$ as in the case [14,
Theorem 2], the 3rd step can be performed by a simple application of standard
numerical solver for ODEs. However, there are some HGM problems for which we
encounter difficulties in the 3rd step. The difficulty 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 a false one. We have this
difficulty in our analysis of the outage probability of MIMO WiFi systems (see
the section 6 and [8]) and in solving an ODE numerically of the expectation of
the Euler characteristic of a random manifold (see the section 7 and [30]). In
principle, if we give an initial value vector exactly and the time step is
sufficiently 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.

##  2 Methods to test

We will compare some methods mainly from a point of view of using for the 3rd
step of the HGM111Programs of this paper can be obtainable from
http://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 some cases. . In
the HGM, the function $Z(t)$ to evaluate has an integral representation and
satisfies an ODE

$Lf=b,\quad L=\sum_{k=0}^{r}c_{k}(t)\partial_{t}^{k},\
\partial_{t}=\frac{d}{dt}$ (1)

where $b(t)$ is an inhomogeneous term and $f(t)$ is an unknown function. See,
e.g., [8], [14], and papers in [15] as to examples. The differential equation
(1) might have a solution $u$ which dominates $Z$ at $+\infty$. In other
words, there might be a solution $u$ such that $|u/Z|\rightarrow\infty$ when
$t\rightarrow+\infty$. In such cases, it is not easy to obtain numerical
values of $Z(t)$ globally by solving the ODE even when the almost accurate
initial values of $Z(t),Z^{\prime}(t),\ldots$ are known. See Examples 1 and 2.
We call such case the unstable HGM problem. We want to obtain correct values
of $Z$ as large interval of $t$ as possible.

Here is a list of methods we try.

  1. 1.

The Runge-Kutta method (see, e.g., [13, Chapter 2]).

  2. 2.

The implicit Runge-Kutta method (see, e.g., [13, Chapters 2 and 4]).

  3. 3.

Multi-step methods (see, e.g., [13, Chapters 3 and 5]).

  4. 4.

The spectral method in the approximation theory (see, e.g., [2], [23], [31],
Section 5).

  5. 5.

Defusing method (filter method): restricting the initial value vector to a sub
linear space (see, e.g., Section 3).

  6. 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], [10] and the Section 4).

We compare these methods by different implementations on different languages:
[2], [6], [9], [12], [20], [21], [26], [27], [28], [29]. 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 target function $Z(t)$ at any point by numerical integration methods. Of
course, the value obtained is not always 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 (filter method)

In the HGM, we do not always have an exact initial value vector. We often have
the case that the target function (the target integral) have a moderate
growth, but the ODE has rapid growth solutions, 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
propose a heuristic method, which we will call defusing method, to avoid a
blow-up of a solution under this situation. Any solution of the ODE can be
expressed as a linear combination of basis functions of different growth
order. The defusing method removes some basis functions of specified growth
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 some interval of $t$. We put this function $f(t)$. The HGM
gives an algebraic method to find an ordinary differential equation (ODE) for
$F(t)=(f(t),f^{\prime}(t),f^{\prime\prime}(t),\ldots,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 differential equation, we can apply algorithms (see, e.g., [5],
[17] ) to find a basis of series solutions at $t=\infty$. Let $F_{*j}$,
$j=1,\ldots,r$ be the basis of series solutions around $t=\infty$ of the ODE.
We denote by $F_{ij}$ is the $i$th component of $F_{*j}$ and $F_{i}$ is the
$i$th component of $F$. Any solution of the ODE is written as a linear
combination of $F_{*j}$’s over ${\bf C}$. We assume

$F_{11}>F_{12}>\cdots>F_{1m}>\cdots>F_{1r}>0$

with different growth orders when $t\rightarrow+\infty$. We suppose that

$\left|\frac{\mbox{numerical approximation of
$F_{1}(t)$}}{F_{1j}(t)}\right|\rightarrow 0,\quad t\rightarrow+\infty$ (2)

for $1\leq j<m$. See (37) as an example. Under these assumptions, we can
conclude that $F(t)$ is approximated as a linear combination of
$F_{*m},\ldots,F_{*r}$. The coefficients 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=+\infty$ by a set of fundamental solutions of the ODE. We call this method
the defusing method by series solutions around $t=+\infty$. See the Example
11.

We want to find a numerical solution of the initial value problem of the
ordinary differential equation (ODE)

$\displaystyle\frac{dF}{dt}$ $\displaystyle=$ $\displaystyle P(t)F$ (3)
$\displaystyle F(t_{0})$ $\displaystyle=$ $\displaystyle F_{0}^{\rm
true}\in{\bf R}^{r}$ (4)

where $P(t)$ is an $r\times r$ matrix, $F(t)$ is a column vector function of
size $r$, and $F_{0}^{\rm true}$ is the initial value of $F$ at $t=t_{0}$. Let
us explain a defusing method for the initial value problem. Solving this
problem is the final step of the holonomic gradient method (HGM) [15]. We
often encounter the following situation in the final step.

###### Situation 1

  1. 1.

The initial value has at most 3 digits of accuracy. We denote this initial
value $F_{0}$.

  2. 2.

The property $|F|\rightarrow 0$ when $t\rightarrow+\infty$ is known, e.g.,
from a background of the statistics.

  3. 3.

There exists a solution ${\tilde{F}}$ of (3) such that
$|{\tilde{F}}|\rightarrow+\infty$ or non-zero finite value when
$t\rightarrow+\infty$.

Under this situation, the HGM works only for a very short interval of $t$
because the error of the initial value vector makes the fake solution
${\tilde{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

$\frac{d}{dt}F=\left(\begin{array}[]{ccc}-1&1&0\\\ 0&-1&1\\\
0&0&0\end{array}\right)F$

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 $F_{1}$. Add some errors
$(1,10^{-30},10^{-30})^{T}$ to the initial value. Then, we have

$t$ | value $F_{1}$ by RK | difference $F_{1}-F^{1}_{1}$
---|---|---
50 | 1.92827e-22 | 9.99959e-31
60 | 8.75556e-27 | 1.00000e-30
70 | 1.39737e-30 | 1.00000e-30
80 | 1.00002e-30 | 1.00000e-30

We can see the instability.

###### Example 2

$P(t)=\left(\begin{array}[]{cc}0&1\\\ t&0\\\ \end{array}\right).$

This differential equation is obtained from the Airy differential equation

$y^{\prime\prime}(t)-ty(t)=0$

by putting $F=(y(t),y^{\prime}(t))^{T}$. It is well-known that the Airy
function

${\rm
Ai}(t)=\frac{1}{\pi}\lim_{b\rightarrow+\infty}\int_{0}^{b}\cos\left(\frac{s^{3}}{3}+ts\right)ds$

is a solution of the Airy differential equation and

$\displaystyle{\rm Ai}(0)$ $\displaystyle=$
$\displaystyle\frac{1}{3^{2/3}\Gamma(2/3)}=0.355028053887817\cdots$
$\displaystyle{\rm Ai}^{\prime}(0)$ $\displaystyle=$
$\displaystyle\frac{1}{3^{1/3}\Gamma(1/3)}=-0.258819403792807\cdots$
$\displaystyle\lim_{t\rightarrow+\infty}{\rm Ai}(t)$ $\displaystyle=$
$\displaystyle 0$ $\displaystyle\lim_{t\rightarrow+\infty}{\rm
Ai}^{\prime}(t)$ $\displaystyle=$ $\displaystyle 0$

Figure 1 is a graph of Airy Ai function and Airy Bi function. The function
$F(t)=({\rm Ai}(t),{\rm Ai}^{\prime}(t))^{T}$ satisfies the condition 2 of the
Situation 1 of the instability problem.

Figure 1: Airy ${\rm Ai}(x)$ and ${\rm 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 differential equations (see, e.g., the
manual DEtools/formal_sol of Maple [19] and its references on the theory,
which has a long history). In fact, the general solution of the Airy
differential equation is expressed as

$\displaystyle
C_{1}t^{-1/4}\exp\left(-\frac{2}{3}t^{3/2}\right)(1+O(t^{-3/2}))$
$\displaystyle+$ $\displaystyle
C_{2}t^{-1/4}\exp\left(\frac{2}{3}t^{3/2}\right)(1+O(t^{-3/2}))$

when $t\rightarrow+\infty$ where $C_{i}$’s are arbitrary constants.

We note that the high precision evaluation of the Airy function is studied by
several methods (see, e.g., [7] and its references). Some mathematical
software systems have evaluation functions of the Airy function. For example,
N[AiryAi[10]] gives the value of ${\rm Ai}(10)$ on Mathematica. By utilizing
these advanced evaluation methods, we use the Airy differential equation for
our test case to 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

$F_{k+1}=Q(k,h)F_{k}$ (5)

from $F_{0}$ where $Q(k,h)$222It was denoted by $Q(t_{0}+kh,h)$ in the
previous section. We denote $Q(t_{0}+kh,h)$ by $Q(k,h)$ as long as no
confusion arises. is an $r\times r$ matrix determined by a numerical scheme
and $h$ is a small number The vector $F_{k}$ is an approximate value of $F(t)$
at $t=t_{k}=t_{0}+hk$.

###### Example 3

The Euler method assumes $dF/dt(t)$ is approximated by $(F(t+h)-F(t))/h$ and
the scheme of this method is

$F_{k+1}=(E+hP(t_{k}))F_{k}$

where $E$ is the $r\times r$ identity matrix.

In case that the initial value vector $F_{0}$ contains an error, the error may
generate a blow-up solution ${\tilde{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)\cdots 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 our presentation. The following
heuristic algorithm avoids to get the blow-up solution.

###### Algorithm 1

  1. 1.

Obtain eigenvalues $\lambda_{1}>\lambda_{2}>\cdots>\lambda_{r}>0$ of $Q$ and
the corresponding eigenvectors $v_{1},\ldots,v_{r}$.

  2. 2.

Let $\lambda_{m}$ be the first non-positive eigenvalue.

  3. 3.

Express the initial value vector $F_{0}$ containing errors in terms of
$v_{i}$’s as

$F_{0}=f_{1}v_{1}+\cdots+f_{r}v_{r},\quad f_{i}\in{\bf R}$ (7)

  4. 4.

Choose a constant $c$ such that
$F_{0}^{\prime}:=c(f_{m}v_{m}+\cdots+f_{r}v_{r})$ approximates $F_{0}$.

  5. 5.

Determine $F_{N}$ by $F_{N}=QF_{0}^{\prime}$ with the new initial value vector
$F_{0}^{\prime}$.

We call this algorithm the defusing method for initial value problem. This is
a heuristic algorithm and the vector $F_{0}^{\prime}$ gives a better
approximation of the initial value vector than $F_{0}$ in several examples and
we can avoid the blow-up of the numerical solution.

###### Example 4

We set $t_{0}=0$, $h=10^{-3}$, $N=10\times 10^{3}$ and use the 4-th order
Runge-Kutta scheme. We have $\lambda_{1}=9.708\times 10^{9}$,
$v_{1}=(-5.097,-159.919)^{T}$ and $\lambda_{2}=3.247\times 10^{-7}$,
$v_{2}=(-5.09798,37.164813649680576037539971418209465086)^{T}=(a,b)$. Then,
$m=2$. We assume the $3$ digits accuracy of the value ${\rm Ai}(0)$ as $0.355$
and set $F_{0}^{\prime}=(0.355,0.355b/a)$. Then, the obtained value $F_{5000}$
at $t=5$ is $(0.000108088745179140,-0.000246853220440734)$. We have the
following accurate value by Mathematica

    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 ${\rm Ai}(5)$. On
the other hand, we appy the 4th order Runge-Kutta method with $h=10^{-3}$ for
$F_{0}=(0.355,-0.259)^{T}$, which has the accuracy of $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 an error estimate of the solution by our
algorithm. Let $|\cdot|$ be the Eucledian norm.

###### Lemma 1

Let $F(t)$ be the solution. When $|QF_{0}^{\rm true}-F(Nh)|<\delta$ holds, we
have

$|QF_{0}^{\prime}-F(Nh)|<|QF_{0}^{\prime}|+|F(Nh)|+2\delta$ (8)

for any $F_{0}^{\prime}\in{\bf R}^{n}$.

Proof. It is a consequence of the triangular inequality. In fact, we have

$\displaystyle|QF_{0}^{\prime}-F(Nh)|$ $\displaystyle=$
$\displaystyle|QF_{0}^{\prime}-QF_{0}^{\rm true}+QF_{0}^{\rm true}-F(Nh)|$
$\displaystyle\leq$ $\displaystyle|QF_{0}^{\prime}-QF_{0}^{\rm
true}|+|QF_{0}^{\rm true}-F(Nh)|$ $\displaystyle\leq$
$\displaystyle|QF_{0}^{\prime}|+|QF_{0}^{\rm true}|+\delta$
$\displaystyle\leq$ $\displaystyle|QF_{0}^{\prime}|+|F(Nh)|+2\delta$

$//$

Under the Situation 1, $|F(Nh)|$ is small enough. Then, it follows from the
Lemma that $|QF_{0}^{\prime}|$ should be small. In this context, we can give
an error estimate of our algorithm. However, our numerical experiments present
that the algorithm shows a better behavior than this theoretical error
estimate. Then, we would like to classify our defusing method as a heuristic
method.

###### Example 5

(Solving Airy differential equation by the defusing method.) (
defusing/intro/2023-07-21-airy.rr) Give initial values at $t=-20$ as
$F_{0}=[-0.17640612707798468959,0.89286285673647123840]$ (${\rm Ai}(-20)$ and
${\rm Ai}^{\prime}(-20)$). It keeps $3$ digits of accuracy at $t=5$ even when
we start at $t=-20$, but the accuracy is broken at $t=6$; The defusing method
gives the value $1.09\times 10^{-5}$ at $t=6$, but ${\rm Ai}(6)\sim 9.95\times
10^{-6}$. The graph in Figure 2 looks nice without a blowing-up, but we have
to be careful that it is not accurate. For example, the defusing method gives
the value $5.19\times 10^{-49}$ at $t=30$, but ${\rm Ai}(30)\sim 3.21\times
10^{-49}$. We need to correct errors at several values of $t$ by different
methods of evaluation for this problem. On the other hand, it works more
accurately for the function $H^{k}_{n}$. See Example 19. It is a future
problem to give a practical error bound of the defusing method.

Figure 2: Solving initial value problem of Airy differential equation by
defusing, $t\in[-20,30]$

##  4 Sparse interpolation/extrapolation methods

The methods in this section are standard approximation methods in functional
analysis. They are not new and are not named, but they are very useful for
unstable HGM problems. It is for this reason that we name methods of this
section sparse interpolation/extrapolation methods.

We solve the ODE (1) of rank $r$ when (approximate) values of $f(t)=Z(t)$ at
$t=p_{1},p_{2},\ldots,p_{r}$ are known. We call the points
$(p_{i},q_{i}),q_{i}=f(p_{i})$ 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 of $f$ on an
interval $[t_{s},t_{e}]$ is as follows. Devide $[t_{s},t_{e}]$ into $N$
intervals. Let $t_{i}$ be $t_{s}+hi$ where $h=(t_{e}-t_{s})/N$. We assume that
the set of all $p_{j}$’s are a subset of the set $\\{t_{i}\\}$. We denote by
$f_{i}$ the value of $f(t_{i})$. We introduce the backward shift operator
$\nabla f_{i}=f_{i}-f_{i-1}$. Then $\frac{\nabla^{k}f_{i}}{h^{k}}$ is
approximately equal to $f^{(k)}(t_{i})$ or $f^{(k)}(t_{i-1})$ or $\ldots$ or
$f^{(k)}(t_{i-k})$. For example, when $k=1$, we have
$f(t_{i})-f(t_{i-1})=f^{\prime}(t_{i})h+O(h^{2})$ if we make the Taylor
expansion of $f(t_{i-1})=f(t_{i}-h)$ at $t=t_{i}$ and
$f(t_{i})-f(t_{i-1})=f^{\prime}(t_{i-1})h+O(h^{2})$ if we make the Taylor
expansion of $f(t_{i})=f(t_{i-1}+h)$ at $t=t_{i-1}$. We have

$\sum_{k=0}^{r}c_{k}(t_{i})\frac{\nabla^{k}f_{i+s_{k}}}{h^{k}}\sim
b(t_{i}),\quad 0\leq s_{k}\leq k.$ (9)

Here $s_{i}$ is an integer to choose an approximation of $f^{(k)}(t_{i})$. By
assuming the left hand side and the right hand side are equal and giving
values $f_{j}$ for $t=p_{j},1\leq j\leq r$, we have a system of linear
equations of the form

$A\left(\begin{array}[]{c}f_{0}\\\ f_{1}\\\ \cdot\\\ \cdot\\\ \cdot\\\
f_{N}\end{array}\right)=B$ (10)

where $A$ is a $(N+1)\times(N+1)$ matrix and $B$ is a column vector of length
$N+1$. Solving this equation, we obtain approximate values of $f(t_{i})$. We
call this method the sparse interpolation/extrapolation method A for HGM in
this paper.

###### Example 6

We solve $Lf=0$, ( intro/2021-06-11-sparse-interp.ipynb)
$L=(\partial_{t}-1)(\partial_{t}^{2}-t)=\partial_{t}^{3}-\partial_{2}^{2}-t\partial_{t}+t-1$
on $t\in[-9,0]$ with $f(-9)={\rm Ai}(-9)\sim-0.0221337215473414$, $f(-4)={\rm
Ai}(-4)\sim-0.0702655329492895$, $f(0)={\rm Ai}(0)\sim 0.355028053887817$ and
$N=100$. Then, we obtain approximate values of the Airy function ${\rm
Ai}(t)$. The numerical result shows that we do not have a false solution in
spite of the factor $\partial_{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 coefficients to
express the solution of the generalized boundary conditions as the linear
combination of the fundamental solutions.

Let us introduce the sparse interpolation/extrapolation method B for HGM. This
is also a standard method and using it in the HGM is very useful. We minimize

$\int_{t_{s}}^{t_{e}}|Lf(t)-b(t)|^{2}d\mu(t)$ (11)

where $d\mu(t)$ is a measure in the $t$ space. We approximate this integral. A
numerical integration for a function $g$ can be expressed as

$I_{N}(g)=\sum_{j=0}^{N}T_{j}g(t_{j})$ (12)

where $t_{0}=t_{2}<t_{1}<\cdots<t_{N-1}<t_{N}=t_{e}$ and $T_{j}\in{\bf
R}_{\geq 0}$. We fix a numerical integration method. For example, the
trapezoidal method can be expressed as $h=\frac{t_{e}-t_{s}}{N}$,
$t_{i}=t_{s}+hi$, $T_{j}=h$ for $1\leq j<N$ and $T_{0}=T_{N}=h/2$. We
approximately expand the solution $f$ by a given basis functions
$\\{e_{k}(t)\\}$, $k=0,1,\ldots,M$ as

$f(t)=\sum_{k=0}^{M}f_{k}e_{k}(t),\quad f_{k}\in{\bf R}.$ (13)

We put $M+1=m$ to be compatible with the variable name in our program. Put
this expression into $Lf=b$. We minimize the following loss function, which is
approximately equal to (11),

$\displaystyle\ell(\\{f_{k}\\})$ $\displaystyle=$
$\displaystyle\sum_{j=0}^{N}|(Lf)(t_{j})-b(t_{j})|^{2}T_{j}$ $\displaystyle=$
$\displaystyle\sum_{j=0}^{N}\left|\sqrt{T_{j}}\sum_{k=0}^{M}f_{k}(Le_{k})(t_{j})-\sqrt{T_{j}}b(t_{j})\right|^{2}$

under the constraints333If data are given in derivative values of $f$, these
constraints become relations among $q_{i}$’s. at data points

$\sum_{k=0}^{M}f_{k}\cdot e_{k}(p_{i})=q_{i},\quad i=1,2,\ldots,r.$ (15)

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

###### Lemma 2

The norm $\int_{t_{s}}^{t_{e}}|Lf(t)-b(t)|^{2}d\mu(t)$ is bounded by
$\ell(\\{f_{k}\\})+e_{N}(|Lf-b|^{2})$ where $e_{N}(|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 can estimate the $L^{2}$ norm of $Lf-b$ by this Lemma.

We can also consider a least mean square problem with no constraints by the
loss function

${\tilde{\ell}}(\\{f_{k}\\})=\alpha\ell(\\{f_{k}\\})+\beta\sum_{i=1}^{r}\left(\sum_{k=0}^{M}f_{k}\cdot
e_{k}(p_{i})-q_{i}\right)^{2}+\gamma\sum_{i=0}^{N}f_{i}^{2}.$ (16)

Here $\alpha,\beta,\gamma$ are hyperparameters for the optimization. In
particular, $\gamma>0$ is used to avoid an overfitting.

###### Example 7

We solve the differential equation of Example 6. 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. See the Example 16 as
to an application of this method to a 4th order ODE. An application of this
method to a problem on random matrices is given in Example 20.

Let us introduce the sparse interpolation/extrapolation method C for HGM. We
use the basis functions $e_{j}(t)$ as in the method C. We construct an
$m\times m$ symmetric matrix $S$ of which $(i,j)$ element is

$\int_{t_{s}}^{t_{e}}(Le_{i}(t))(Le_{j}(t))d\mu(t)$ (17)

Let $F$ be the column vector of length $m=M+1$:
$F=(f_{0},f_{1},\ldots,f_{M})^{T}$. Then, we have

$\int_{t_{s}}^{t_{e}}\left(L\sum_{j=0}^{M}f_{j}e_{j}(t)\right)^{2}d\mu(t)=F^{T}SF.$
(18)

It follows from this idenity that when $\sum_{j=0}^{M}f_{j}e_{j}(t)$ is an
approximate solution of $Lf=0$, we may expect the value of the quadratic form
$F^{T}SF$ is small. Let $(p_{k},q_{k})$’s be given data points of $(t,f(t))$.
The method C finds the vector $F$ by solving the quadratic programming problem
of minimizing

$\displaystyle
F^{T}SF+\sum_{k=1}^{r}\left(\sum_{j=0}^{M}f_{j}e_{j}(p_{k})-q_{k}\right)^{2}$
(19) $\displaystyle=$ $\displaystyle
F^{T}SF+\sum_{k=1}^{r}\left(F^{T}S_{k}^{2}F-2q_{k}S_{k}F+q_{k}^{2}\right)$
$\displaystyle=$ $\displaystyle
F^{T}SF+F^{T}(P_{e}^{T}P_{e})F-2P_{e}^{T}Q+Q^{T}Q$ (20)

where $S_{k}={\rm diag}(e_{0}(p_{k}),\ldots,e_{M}(p_{k}))$ and

$P_{e}=\left(\begin{array}[]{ccc}e_{0}(p_{1})&\cdots&e_{M}(p_{1})\\\
e_{0}(p_{2})&\cdots&e_{M}(p_{2})\\\ \cdot&\cdots&\cdot\\\
e_{0}(p_{r})&\cdots&e_{M}(p_{r})\\\ \end{array}\right),\quad
Q=(q_{1},\ldots,q_{r})^{T}$ (21)

Or, we may solve the quadratic programming problem of minimizing $F^{T}SF$
under the constraints

$\sum_{j=0}^{M}f_{j}e_{j}(p_{k})-q_{k}=0,\quad k=1,\ldots,r.$ (22)

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

##  5 The Chebyshef function method

The approximation by Chebyshef functions is known to be very good. We briefly
overview the Chebyshef function methods for ODE’s and show that it falls in
the framework of sparse interpolation/extrapolation methods. A comprehensive
reference on the method is the book by Trefethen [31]. A nice package chebfun
[6] for Matlab has been developed, see also
https://en.wikipedia.org/wiki/Chebfun. The chebfun project was initiated in
2002 by Lloyd N. Trefethen and his student Zachary Battles.

A different solver with validation and Chebyshef functions is proposed in [4].
The advantage of the method is that matrices in the solver are banded and
validation is given. We will test this method for the HGM in a next paper.

The $n$-th Chebyshef function (polynomial) is

$T_{n}(x)=\cos(n\theta),\quad x=\cos\theta$ (23)

The extreme points of the curve $y=T_{n}(x)$ in $[-1,1]$, which we mean points
that take the values $y=1$ or $y=-1$, are called $n+1$ Chebyshef points (of
the second kind) of $T_{n}$. For example, $T_{2}(x)=2x^{2}-1$ and the
Chebyshef points are $\\{-1,0,1\\}$.

Let $f(x)$ be a function. Fix the set of Chebyshef points for $T_{n}(x)$. Let
the value of $f$ at Chebyshef point $x_{j}$ be $f_{j}$. The (degree $n$)
Chebyshev interpolant is

$p(x)={\sum_{j=0}^{n}}^{\prime}\frac{(-1)^{j}f_{j}}{x-x_{j}}/{\sum_{j=0}^{n}}^{\prime}\frac{(-1)^{j}}{x-x_{j}}$
(24)

The primes on the summation signes signify that the terms $j=0$ and $j=n$ are
multiplied by $1/2$. It takes the value $p(x_{j})=f_{j}$ and $p(x)$ agrees
with the Lagrange interpolation polynomial for $(x_{j},f_{j})$’s. See, e.g.,
[31, Th 5.1].

###### Theorem 1

(Bernstein 1911, 1912. See, e.g., Th 8.2, Th 8.3 in [31]) If $f$ is analytic
on $[-1,1]$, its Chebyshef coefficients $a_{k}$ decrease geometrically. If $f$
is analytic and $|f|\leq M$ in the Bernstein $\rho$-ellipse about $[-1,1]$,
then $|a_{k}|<2M\rho^{-k}$. The degree $n$ Chebyshev interpolant has accuracy
$O(M\rho^{-n})$ by the sup norm.

Here, $a_{k}$ is defined by
$a_{k}=\frac{2}{\pi}\int_{-1}^{1}\frac{f(x)T_{k}(x)}{\sqrt{1-x^{2}}}dx$ and
satisfies $f=\sum_{k=0}^{\infty}a_{k}T_{k}(x)$ when $f$ is Lipschitz
continuous. The Bernstein $\rho$-ellipse is defined as follows. Consider the
radius $\rho$ circle in the $z$-plane. Map it by $x=(z+z^{-1})/2$ and then we
obtain the $\rho$-ellipse. Its foci is at $-1$ and $1$. $\rho$ is the sum of
semi-major and semi-minor axes.

Chapter 16 of [31] shows that Chebyshev interpolants are often as good as the
best approximation in practice.

Let us proceed on explaining briefly a method of solving ODE’s with Chebyshef
functions and the implementation of chebfun [6]. The point is an approximation
of differential operators in terms of Chebyshef friendly matrices.

Let $X=\\{X_{0},X_{1},\ldots,X_{n-1}\\}$ be the set of the $n$ Chebyshef
points (of the second kind) for the Chebyshef function $T_{n-1}$. The $i$-th
element of $X$ $(0\leq i<n)$ is given as

$\cos\left(\pi\frac{n-i-1}{n-1}\right).$ (25)

We denote by $t$ the independent variable for ODE’s instead of $x$. We denote
by $\ell_{j}(X;t)$ the $j$-th polynomial of the Lagrange interpolation for
$X$. We have $\ell_{j}(X;X_{j})=1$ and $\ell_{j}(X;X_{i})=0$ when $i\not=j$.

When $f(t)$ is a polynomial of degree $n-1$, we have

$f(t)=\sum_{j=0}^{n-1}\ell_{j}(X;t)f_{j}$ (26)

by putting $f_{j}=f(X_{j})$.

Now, let $Y$ be the set of the $(n-m)$ Chebyshef points for $T_{n-m-1}$ where
$m\geq 0$. We approximate $f(t)$ by the values at $Y$, which is called “down-
sampling” in [11]. Since we have

$f^{(m)}(t)=\sum_{j=0}^{n-1}\ell_{j}^{(m)}(X;t)f_{j}$ (27)

and it is a polynomial of degree $n-m-1$, we have

$\displaystyle f^{(m)}(t)$ $\displaystyle=$
$\displaystyle\sum_{k=0}^{n-m-1}\ell_{k}(Y;t)\left(\sum_{j=0}^{n-1}\ell_{j}^{(m)}(X;Y_{k})f_{j}\right)$
(28) $\displaystyle=$
$\displaystyle\sum_{j=0}^{n-1}\left(\sum_{k=0}^{n-m-1}\ell_{k}(Y;t)\ell_{j}^{(m)}(X;Y_{k})\right)f_{j}$
(29)

###### Definition 1

The $(n-m)\times n$ matrix with $(i,j)$ entries

$\sum_{k=0}^{n-m-1}\ell_{k}(Y;Y_{i})\ell_{j}^{(s)}(X;Y_{k})$ (30)

is denoted by $M(n-m,n;s)$ .

For $s\leq m$, the $s$-th derivative value of $f(t)$ at $Y_{i}$ is
$(\mbox{$i$-th row of $M(n-m,n;s)$})\cdot(f_{0},\ldots,f_{n-1})^{T}$. Then,
the matrix can be used to approximate $\frac{d^{s}f}{dt^{s}}$. The ODE

$\sum_{i=0}^{r}c_{i}(t)\frac{d^{i}f}{dt^{t}}=0$

is approximately translated into the linear equation

$\left(\sum_{i=0}^{r}{\rm
diag}(c_{i}(Y_{0}),c_{i}(Y_{1}),\ldots,c_{i}(Y_{n-r-1}))M(n-r,n;i)\right)(f_{0},f_{1},\ldots,f_{n-1})^{T}=0$
(31)

in the chebfun system where $Y$ is the $n-r$ Chebyshef points for $T_{n-r-1}$.
In the chebfun system, the matrix $M(n-m,n;s)$ is obtained by the function
diffmat(n-m,n,s). The solution $f(t)$ is recovered by the Chebyshef
interpolant (24) from $f_{j}$’s.

The system of linear equations (31) is used in the Chebyshef function method.

###### Example 8

The Airy equation

$f^{\prime\prime}-tf=0$

is solved by solving a linear equation derived by $M(n-2,n;2)$ and
$M(n-2,n;0)$ with two boundary conditions. Symbolically, we solve

$\left(M(n-2,n;2)-{\rm diag}(Y)M(n-2,n;0)\right)F=0$ (32)

where $F=(f_{0},f_{1},\ldots,f_{n-1})^{T}$ with giving, e.g., values of
$f_{0}$ and $f_{n-1}$.
See https://www.chebfun.org/examples/ode-linear/SpectralDisc.html.

Figure 3: Solving the Airy differential equation by chebfun

###### Example 9

Initial value problem for Airy ${\rm Ai}(t)$. (
defusing/intro/y2023_07_16_airy_initial_value.m) We give the initial value at
$t=-20$ as ${\rm Ai}(-20)=-0.176406127077984689590192292219$ and ${\rm
Ai}^{\prime}(-20)=0.892862856736471238398409934114$. These values are
evaluated by Mathematica. Chebfun gives reasonable values upto $t=9$, but
divergent values appear when $t$ is larger than $9$. See the left graph of
Figure 3.

###### Example 10

Boundary value problem for Airy ${\rm Ai}(t)$. (
defusing/intro/y2023_07_16_airy_boundary_value.m) We give the boundary value
at $t=-20$ and $t=11$ as ${\rm Ai}(-20)=-0.176406127077984689590192292219$ and
${\rm Ai}(-11)=4.22627586496035959129883545080\times 10^{-12}$. Divergent
values do not appear. See the right graph of Figure 3.

As we have seen in two examples, the Chebyshef function method (solving
boundary value problem) proposed and implemented in [11] works well, but
solving the initial value problem is affected by numerical error.

Finally, we note that the Chebyshef function method can be regarded as a
special case of the sparse interpolation/extrapolation method B. In fact, the
numerical integration scheme of the Chebyshef quadrature is

$\int_{-1}^{1}\sqrt{1-t^{2}}g(t)dt\sim\sum_{i=1}^{k-1}T_{i}g(Y_{i})$ (33)

where $Y$ is the set of the Chebyshef points (25) for $T_{k}$ and the weight
$T_{i}$ is

$T_{i}=\frac{\pi}{k}\sin^{2}\left(\frac{i}{k}\pi\right)$

Put $g(t)=|Lf|^{2}$ and $d\mu(t)=\sqrt{1-t^{2}}dt$. Since the left hand side
of (31) are values at the set of Chebyshef points $Y$, assuming it is equal to
the zero vector is equivalent to that the integral by the Chebyshef quadrature
over $Y$ is equal to zero.

##  6 Tests of the methods to the function $H^{k}_{n}$

Let $n$ and $k$ be positive integers. We define the function $H^{k}_{n}(x,y)$
by

$\displaystyle H^{k}_{n}(x,y)$ $\displaystyle=$
$\displaystyle\int_{0}^{x}t^{k}\exp(-t){}_{0}F_{1}(;n;yt)dt$ (34)
$\displaystyle=$
$\displaystyle\frac{\Gamma(n)}{\sqrt{\pi}\Gamma(n-1/2)}\int_{D(x)}t^{k}(1-s^{2})^{n-3/2}\exp(-t-2s\sqrt{yt})dtds$
$\displaystyle\ \mbox{where}\ D(x)=\\{(t,s)\in[0,x]\times[-1,1]\\}$

This function appears in studies of the outage probability of MIMO WiFi
systems (see, e.g., [8]). We will compare the methods in the section 2 for
this function on several systems. The function satisfies the following system
of linear partial differential equations.

###### Proposition 1

([8])

  1. 1.

The function $u=H_{n}^{k}(x,y)$ satisfies

$\displaystyle\\{\theta_{y}(\theta_{y}+n-1)+y(\theta_{x}-\theta_{y}-k-1)\\}\bullet
u$ $\displaystyle=$ $\displaystyle 0,$
$\displaystyle(\theta_{x}-\theta_{y}-k-1+x)\,\theta_{x}\bullet u$
$\displaystyle=$ $\displaystyle 0.$

where $\theta_{x}=x\frac{\partial}{\partial
x}$,$\theta_{y}=y\frac{\partial}{\partial y}$. The holonomic rank of this
system is $4$.

  2. 2.

The function $u$ is a solution of the following ODE with respect to $y$.

$y^{2}\partial_{y}^{4}+(-y+2n+2)y\partial_{y}^{3}+(-yx+(-k-n-3)y+n(n+1))\partial_{y}^{2}+((y-n)x-n(k+2))\partial_{y}+(k+1)x$
(36)

( Prog_paper/19-a19-n-pf.rr)

When $y\rightarrow+\infty$, solutions of the system has the following
asymptotic behavior. It is shown by the DEtools[formal_sol] function of Maple
[19]. ( Hkn/a4-f2.ml)

$\displaystyle h_{1}$ $\displaystyle=$
$\displaystyle(xy)^{-1/2(1/2+n)}\exp(-2(xy)^{1/2})(1+O(1/y^{1/2})),$
$\displaystyle h_{2}$ $\displaystyle=$ $\displaystyle y^{-k-1}(1+O(1/y)),$
$\displaystyle h_{3}$ $\displaystyle=$
$\displaystyle(xy)^{-1/2(1/2+n)}\exp(2(xy)^{1/2})(1+O(1/y^{1/2})),$
$\displaystyle h_{4}$ $\displaystyle=$ $\displaystyle
y^{1-n+k}\exp(y)(1+O(1/y)),$

Which is the asymptotic behavior of the function $H^{k}_{n}(x,y)$ when $x$ is
fixed? We compare the value of $h_{4}$ and the value by a numerical
integration in Mathematica444The 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 hypergeometric functions which matches to Mathematica in all ranges of
parameters and independent variables. . 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]

$y$ | Ratio
---|---
1000 | 7.36595030875893e-452
2000 | 2.64621603289928e-881
3000 | 2.67723893601667e-1311

where

$\mbox{Ratio}=(\mbox{$H^{10}_{1}(1/2,y)$})/(y^{1-n+k}\exp(y)).$ (37)

This computational experiments suggest that $H^{k}_{n}$ is expressed by
$h_{1},h_{2},h_{3}$ without the dominant component $h_{4}$ as explained in
(2).

###### Example 11

We approximate $H^{10}_{1}(1,y)$ by $h_{3}$ of the asymptotic series terms
truncated by $O(y^{-3})$. ( Hkn/2021-06-06-asymp.rr) Let $y_{s}$ be a
sufficiently large number. We determine the constant $C_{s}$ by setting
$H^{10}_{1}(1,y_{s})/h_{3}(y_{s})=C_{s}$ and use $C_{s}h_{3}(y)$ as an
approximation of $H^{10}_{1}(1,y)$ (a naive method). Then we have the
following data of relative errors
$r_{e}(y)=\frac{C_{s}h_{3}(y)-H^{10}_{1}(1,y)}{H^{10}_{1}(1,y)}$.

$y_{s}$ | $r_{e}(y_{s}+10)$ | $r_{e}(y_{s}+90)$
---|---|---
$10^{2}$ | $0.0763$ | $0.259$
$10^{4}$ | $5.7\times 10^{-10}$ | $5.09\times 10^{-9}$

The data tells that the naive method works well near $y=10^{4}$, but we need
to use other methods near $y=10^{2}$. The sparse interpolation/extrapolation
method B in the Example 16 is a more refined variation of this method.

Let us apply methods presented in the Section 2 to evaluate $H^{10}_{1}(1,y)$.

###### Example 12

We apply the implicit Runge-Kutta method to $H^{10}_{1}(1,y)$ which is a
solution of the ODE (36). We translate the ODE into a system of ODE for $U$ by
$U=(u,u^{\prime},u^{\prime\prime},u^{\prime\prime\prime})^{T}$ and derive an
ODE for $F=U\exp(-y)y^{-(1-n+k)}$. This is the Gauge transformation by the
exponential part of $h_{4}$ so that the solutions keep bounded values. Then,
the column vector function $F$ satisfies $F^{\prime}=PF$ where

$P=\left(\begin{array}[]{cccc}\frac{-{y}-{k}+{n}-1}{{y}}&1&0&0\\\
0&\frac{-{y}-{k}+{n}-1}{{y}}&1&0\\\ 0&0&\frac{-{y}-{k}+{n}-1}{{y}}&1\\\
\frac{(-{k}-1){x}}{{y}^{2}}&\frac{(-{y}+{n}){x}+{n}{k}+2{n}}{{y}^{2}}&\frac{{y}{x}+({k}+{n}+3){y}-{n}^{2}-{n}}{{y}^{2}}&\frac{-{k}-{n}-3}{{y}}\\\
\end{array}\right).$ (38)

The following data is obtained by the implicit Runge-Kutta method (Gauss
method of order 10, [13, IV.5]) implemented in Mathematica [20]. The Figure 4
presents the value of $\log H^{10}_{1}(1,y)$ (the almost straight graph), the
value by the implicit Runge-Kutta method, and the value by NDSolve with 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 about $y=25$.
This interval $[1,25]$ is larger than the valid interval $[1,4.5]$ of the
ExplicitRungeKutta method of NDSolve.

Figure 4: The Runge-Kutta(yellow green), the implicit Runge-Kutta method(blue)
starting at $y=1$, the exact value(orcher)

When we can obtain exact initial values by the numerical integration (e.g., of
(34)), it will be a good strategy to extrapolate the values in a short
interval by solving the ODE by the implicit Runge-Kutta method. The Figure 5
describes the value of $\log H^{10}_{1}(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 is $10^{-3}$. It works upto about $y=130$.

Figure 5: The implicit Runge-Kutta method starting at $y=100$ vs the exact
value. The right graph is the relative error $(H_{r}-H)/H$ where $H$ is the
exact value and $H_{r}$ is the value by the implicit Runge-Kutta method.

We note that the implicit Runge-Kutta method is used in a large scale
computation with the spectral diferred correction (SDC), see, e.g., [32],
[24]. 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. ( num-ht3/pySDC/airy_playground.py, Airy_implicit.py)

###### Example 13

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
(36). The outputs are remarkable. We solve the ODE on $[1,40]$ with the
initial value at $y=1$ as
$(u,u^{\prime},u^{\prime\prime},u^{\prime\prime\prime})(1)=(0.07810139136088563,0.05096276584900834,0.02050273784371611,0.005887855153702640)$
where $u(y)=H^{10}_{1}(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 than the 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 (36)$\times y^{2}$, the program returns NaN.
Another example is that it did not work for a first order system obtained from
(36). ( Julia/hkn3b.jl)

###### Example 14

We can try several solvers of ODE by solve_ivp on scipy/python [27]. (
Hkn/2021-06-02-hkn.ipynb) We solve the initial value problem of (38) on
$y\in[1,30]$ with the absolute tolerance $0$ and relative tolerance $10^{-3}$
555atol=0, rtol=1e-3. 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 more than $0.3$. Refer to [27] on details on the
following methods. ( Hkn/2021-06-08-hkn-rerr2.ipynb)

Method | Failure $y$ ($\mbox{relative error}>0.3$)
---|---
RK45 | 21.5
Radau(Implicit Runge-Kutta) | 26.1
BDF | 16.8
LSODA(Adams/BDF on ODEPACK) | 20.0

We can also try the deSolve package in R [9]. Some codes are common (e.g.,
FOTRAN package ODEPACK) with the python package and results are analogous with
above. ( Hkn/hkn.r)

###### Example 15

We apply the sparse interpolation/extrapolation method A. In other words, we
solve a generalized boundary value problem of the ODE (36) on
$y\in[10^{4},10^{4}+40]$. ( Hkn/hkn-10-1-10000si-random.py) We give
generalized boundary values of $u(y)=H^{k}_{n}(1,y)\times 10^{-80}$ at
$p_{1}=10^{4},p_{2}=10^{4}+1333h,p_{3}=10^{4}+2000h,p_{4}=10^{4}+4000h=10^{4}+40$
where $h=0.01$. We approximate derivatives as

$\displaystyle u^{(1)}(y)$ $\displaystyle=$
$\displaystyle\frac{1}{h}\left(u(y)-u(y-h)\right)$ (39) $\displaystyle
u^{(2)}(y)$ $\displaystyle=$
$\displaystyle\frac{1}{h^{2}}\left(u(y+h)-2u(y)+u(y-h)\right)$ $\displaystyle
u^{(3)}(y)$ $\displaystyle=$
$\displaystyle\frac{1}{h^{3}}\left(u(y+h)-3u(y)+3u(y-h)-u(y-2h)\right)$
$\displaystyle u^{(4)}(y)$ $\displaystyle=$
$\displaystyle\frac{1}{h^{4}}\left(u(y+2h)-4u(y+h)+6u(y)-4u(y-h)+u(y-2h)\right)$

and solve the linear equation (10). The linear equation is solved by the
function linsolv in the scipy package [26]. The condition number of the matrix
$A$ becomes about $2.6\times 10^{13}$. We give random errors by assuming the
significant digits of the given boundary values are $3$. The Figure 7 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
$p_{i}$’s. ( Hkn/hkn-10-1-10000bv-random.py) For example, we change $p_{2}$ to
$p_{1}+h$ and $p_{3}$ to $p_{4}-h$. This stands for solving the boundary value
problem with the boundary values of $f$ and $f^{\prime}$ at the boundary
$p_{1}$ and $p_{4}$. We solve the linear equation with random errors. We can
see that errors are magnified in middle from the Figure 7.

Figure 6: Solving by the sparse interpolation A |  Figure 7: Solving the boundary value problem with random errors
---|---

Although the boundary value problem with exact boundary values of
$H^{k}_{n}(1,y)$ and its first derivatives can be solved nicely on $[1,40]$,
application of this method on $[1,40]$ with internal reference 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 unacceptable relative error $33.23$ in the interval $[1,6]$, but the
relative error becomes acceptable smaller 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 [10, p.1059] (see also [1,
p.149]). We implemented this method with solve_ivp on scipy/python. It works
for a small interval, e.g., $y\in[1,2]$. But, it fails on $y\in[1,10]$. The
numerical solution of the associated Riccati equation increases from $0$ to
$-1\times 10^{65}$ 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 stop in 8 minutes
on the google colaboratory.

###### Example 16

We apply the sparse interpolation/extrapolation method B explained in Section
4. By utilizing the local asymptotic solution expansion at the infinity of the
ODE (36), we use the set of the $4$ functions

$e_{j}(t)=t^{-3/4}\exp(2t^{1/2})t^{-j/2},\quad j=0,1,2,3,\ t=y$ (40)

as the basis for (13). We will approximate the solution on $[t_{s},t_{e}]$.
The $p_{i}$’s in the constraint (15) are

$p_{0}=t_{s},p_{1}=t_{s}+5,p_{2}=t_{s}+10,\ldots,p_{k}=t_{e}-1.$

We do not use the constraint and add the difference of the approximate value
and the true value to the the loss function as

$\sum_{j=0}^{k}\left((\mbox{approximate value at
$p_{j}$})-H^{10}_{1}(1,p_{j})\right)^{2}$

We expect that this basis will give a good approximation when $t_{s}$ and
$t_{e}$ are large. In fact, when $[t_{s},t_{e}]=[1,40]$, the maximum of the
relative error is $5.47$ by our implementation on least_squares in
scipy/python. On the other hand, this basis gives a good approximation on
$[20,20+40]$ and $[10^{4},10^{4}+40]$. We give random errors to the value of
$q_{j}=H^{10}_{1}(1,p_{j})$ as $q_{j}\times(1+O(10^{-3}))$. It is also robust
to these random errors. ( Hkn/hkn-10-1-10000sib-random.py
Hkn/hkn-10-1-20-60sib-random.py)

$[t_{s},t_{e}]$ | $[20,60]$ | $[10^{4},10^{4}+40]$
---|---|---
max of relative errors with exact $q_{j}$ | $6.21\times 10^{-3}$ | $2.67\times 10^{-12}$
max of relative errors with $30$ random errors | $1.39\times 10^{-2}$ | $4.07\times 10^{-3}$

Figure 8: Solving a boundary value problem by chebfun. Left:
$H^{10}_{1}(1,y)$. Right: $\log H^{10}_{1}(1,y)$. Values should be magnified
by $10^{8678}$.

###### Example 17

(Boundary value problem for $H^{k}_{n}(x,y)$ for $x=1$ and
$y\in[10^{8},10^{8}+2\times 10^{5}]$.)

We give the boundary values of $H^{10}_{1}(1,y)$ and $\frac{\partial
H^{10}_{1}}{\partial y}(1,y)$ at $y=10^{8}$ and $y=10^{8}+2\times 10^{5}$. We
apply the chebfun package for this boundary value problem. (
defusing/Hkn/y2023_07_25_hkn_valid10power8.m) See Figure 8. To check the
accuracy, we compare the values by the chebfun package and by the numerical
integration by Mathematica at $y=10^{8}+200$. The chebfun package keeps $4$
digits accuracy at the point and the ODE is solved in 1.66s666Apple M1, 2020,
Matlab 2022b. On the otherhand, the numerical integration by Mathematica
(2022) ( defusing/Hkn/2023-07-09-hkn-int.m) took 23.58s777AMD EPYC 7552
48-Core Processor, 1499.534MHz.

If we want to use values by the Monte-Carlo integration, the robustness of the
sparse interpolation for boundary values with numerical errors is necessary.

Figure 9: Robustness of chebfun. $H^{10}_{1}(1,y)$ on $[10^{4},10^{4}+40]$
with randomly perturbed boundary values. The values should be magnified by
$10^{82}$.

###### Example 18

(Robustness of the Chebyshev function method.) We solve the boundary value
problem for $H^{10}_{1}(1,y)$ on $[10^{4},10^{4}+40]$. We perturb the boundary
values so that they keep about 3 digits accuracy. (
defusing/Hkn/y2023_07_27_hkn_10power4_robustness.m) Figure 9 is a graph of
$10$ such solutions.

We use a gauge transformed equation to test methods other than Chebyshef
function method. The chebfun system does not work well for the transformed
equation. Multiplication by the transformation seems to give some numerical
errors.

In summary, we should use the implicit Runge-Kutta method or the sparse
interpolation/extrapolation method A when $t_{s}$ and $t_{e}$ are small and
use the sparse interpolation/extrapolation method B or the Chebyshef function
method when they become large for the problem $H^{k}_{n}$. An alternative
choice when $t_{s}$ and $t_{e}$ are relatively small will be the following
defusing method.

###### Example 19

We implement the defusing method in tk_ode_by_mpfr.rr
888http://www.math.kobe-u.ac.jp/OpenXM/Current/doc/asir-
contrib/ja/tk_ode_by_mpfr-html/tk_ode_by_mpfr-ja.html (in Japanese). for the
Risa/Asir [25]. It generates C codes utilizing the MPFR [21] for bigfloat and
the GSL [12] for eigenvalues and eigenvectors. We apply the defusing method
for initial value problem to $H^{10}_{1}(1,y)$ which is a solution of the ODE
(36). 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-ode-assert.hkn1(), tk-ode-assert.hkn2()) The Figure 11 shows
that the adaptive Runge-Kutta method of GSL [12] fails before $y$ becomes
$30$. ( Hkn/defusing/a26-y.c) The Figure 11 presents the relative error of
values by the defusing method and exact values. It shows that the defusing
method works even when $y=10^{3}$.

Figure 10: $\log H^{10}_{1}(1,y)$. Exact value (by numerical integration) and the value by the defusing method agree. The adaptive Runge-Kutta method with the initial relative error $10^{-20}$ (upper curve) does not agree with the exact value when $y$ is larger than about $25$. |  Figure 11: The relative error of $H^{10}_{1}(1,y)$ of the defusing method. The relative error is defined as $(H_{d}-H)/H$ where $H_{d}$ is the value by the defusing method and $H$ is the exact value.
---|---

##  7 Tests of some methods to the function $E[\chi(M_{t})]$

We consider the integral studied in [30]

$\displaystyle{\rm E}[\chi(M_{t})]$ $\displaystyle=$ $\displaystyle
F(s_{1},s_{2},m_{11},m_{21},m_{22};t)$ (41) $\displaystyle=$
$\displaystyle\frac{1}{2}\int_{t}^{\infty}d\sigma\int_{-\infty}^{\infty}db\int_{0}^{2\pi}d\theta\int_{0}^{2\pi}d\phi(\sigma^{2}-b^{2})\frac{s_{1}s_{2}}{(2\pi)^{2}}\exp\Bigl{\\{}-\frac{1}{2}R\Bigr{\\}},$

where

$\displaystyle R$ $\displaystyle=$ $\displaystyle
s_{1}\left(b\sin\theta\sin\phi+\sigma\cos\theta\cos\phi-
m_{11}\right)^{2}+s_{2}\left(\sigma\sin\theta\cos\phi-b\cos\theta\sin\phi-
m_{21}\right)^{2}$
$\displaystyle+s_{1}\left(\sigma\cos\theta\sin\phi-b\sin\theta\cos\phi\right)^{2}+s_{2}\left(b\cos\theta\cos\phi+\sigma\sin\theta\sin\phi-
m_{22}\right)^{2}.$

and

$m_{11}=1,m_{21}=2,m_{22}=3,s_{1}=10^{3},s_{2}=10^{2}.$

By virtue of the Euler characteristic method, the expectation of the Euler
characteristic of a random manifold $M_{t}$ $E[\chi(M_{t})]$ approximates the
probability that the maximal eigenvalue of $2\times 2$ random matrices is less
than $t$. See [30, Sec 3] for details.

The function $E[\chi(M_{t})]$ is a solution of the rank $11$ ODE $Lf=0$
discussed in [30, Example 5]. The operator $L$ is of the form

$\left((-4.72\times
10^{-52}t^{29}+\cdots)\partial_{t}^{10}+\cdots+(-7.78\times
10^{-22}t^{35}+\cdots)\right)\partial_{t}$ (42)

by multiplying a constant $10^{-31}/8.66$ 999This constant is chosen so that
the maximal absolute value of the coefficients of $f_{k}$’s in (4) of the
sparse interpolation/extrapolation method B is $1$ in Example 20. and
$\partial_{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 this ODE is of a form of the singular perturbation problem and
it seems to be hard to solve initial value problems like the Runge-Kutta
method. In [30], since the exact numerical integration of (41) is not easy,
they use values of the Monte-Carlo simulation on $t\in[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 [30] for details. The sparse
interpolation/extrapolation method B overcomes this difficulty.

###### Example 20

( 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

$e_{j}(t)=(t-3.8055)^{j}.$ (43)

Let $p_{i}=3.8+i/1000$, $i=0,1,\ldots,9$. The values $q_{i}$ at $p_{i}$ are

$\displaystyle[0.067160,0.065485,0.064732,0.063315,0.061814,$ $\displaystyle\
\ 0.060477,0.059611,0.058257,0.057520,0.055971]$ (44)

by a Monte-Carlo simulation. ( ec/tryb6/yiex5b.r) $(p_{i},q_{i})$’s are data
points for the sparse extrapolation method B. We divide the interval
$[t_{s},t_{e}]=[3.8,4.0]$ into $200$ segments and use the trapezoidal formula
for $\ell(f)$. The optimal solution

$\displaystyle[f_{0},f_{1},\ldots,f_{29}]$ $\displaystyle=$
$\displaystyle[0.060145405402867516,-1.20074804549872,9.438660716835336,$
$\displaystyle\
-29.022737131194667,-46.606911486598264,587.9594544508735,\ldots]$

can be found in about 18 seconds 101010Debian 10.2 with Intel(R) Xeon(R) CPU
E5-4650 2.70GHz by least_squares of the scipy package [28] on python. The
value of the loss function $\tilde{\ell}$ is equal to $3.85\times 10^{-7}$.
The integral

$\int_{t_{s}}^{t_{e}}\left(LF_{29}\right)^{2}dt\mbox{ where
}F_{29}(t)=\sum_{j=0}^{29}f_{j}e_{j}(t)$ (45)

is equal to $2.58\times 10^{-4}$ and $|F_{29}(p_{i})-q_{i}|<3.18\times
10^{-4}$ for all $i$. The left graph of Figure 12 is the graph of $F_{29}(t)$
and dots are approximate values of $E(\chi(M_{t}))$ by the Monte-Carlo
simulation. The data points are marked with ’x’. The right graph of Figure 12
is the relative error of the values of $F_{29}$ and the values by the
simulation. Although the relative error is rather large, we disagree that the
method is useless. For example, we have $F_{29}(3.935)=0.00241$, then we can
conclude the value 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$.

Figure 12:  The graph of $F_{29}(t)$ and simulation values in the left and
relative errors in the right. The data points are marked with ’x’.

###### Example 21

( ec/tryb5/yi5b-random.py) We give relative errors less than $10^{-3}$ for
$q_{k}$ of the data points $(p_{k},q_{k})$. We use the same scheme as Example
20. We execute $30$ tries and obtain the results in Figure 13. 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.

Figure 13:  The graph of $F_{29}(t)$’s with random errors for data points and
simulation values and relative errors in the right.

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. In fact, the cvxopt package [29] 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 implementation of
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-20-01-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 Differential Equations, 1987, SIAM.
  * [4] F.Bréhard, N.Brisebarre, M.Joldes, Validated and numerically efficient Chebyshev spectral methods for linear ordinary differential equations, ACM Transactions on Mathematical Software, 44 (2018), 1-42. https://doi.org/10.1145/3208103
  * [5] M. Barkatou, An algorithm to Compute the Exponential Part of a Formal Fundamental Matrix Solution of a Linear Differential System, Applicable Algebra in Engineering, Communication, and Computing 8 (1997), 1–23.
  * [6] chebfun for Matlab, https://chebfun.org
  * [7] S. Chevillard, M. Mezzarobba, Multiple-precision evaluation of the Airy Ai function with reduced cancellation, arxiv:1212.4731.
  * [8] 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

  * [9] DeSolve (R-package), Solvers for Initial Value Problems of Differential Equations,
http://desolve.r-forge.r-project.org/

  * [10] L. Dieci, M. Osborne, R. Russell, A Riccati Transformation Method for Solving Linear BFPs I, SIAM Journal on Numerical Analysis 25 (1988), 1055–1073.
  * [11] T.A.Driscoll, N.Hale, Rectangular spectral collocation, IMA Journal of Numerical Analysis, 36 (2016), 108–132 https://doi.org/10.1093/imanum/dru062(2016)
  * [12] GSL, GNU scientific library,
https://www.gnu.org/software/gsl/

  * [13] E.Hailer, S.P.Norsett, G.Wanner, Solving Ordinary Differential Equations I, II. Springer, 1993.
  * [14] 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.
  * [15] References for HGM, http://www.math.kobe-u.ac.jp/OpenXM/Math/hgm/ref-hgm.html
  * [16] T. Hibi and et al, Gröbner Bases: Statistics and Software Systems, 2013, Springer
  * [17] M. van Hoeij, Formal Solutions and Factorization of Differential Operators with Power Series Coefficients, Journal of Symbolic Computation 24 (1997), 1–30.
  * [18] M.Joldes, Validated Numerics: Algorithms and Practical Applications in Aerospace, Proceedings of the 2022 International Symposium on Symbolic and Algebraic Computation, 2022, 1–2. https://doi.org/10.1145/3476446.3535505
  * [19] DEtools in Maple,
https://www.maplesoft.com/support/help/maple/view.aspx?path=DEtools

  * [20] NDSolve in Mathematica,
https://reference.wolfram.com/language/tutorial/NDSolveImplicitRungeKutta.html.

  * [21] The GNU MPFR library, https://www.mpfr.org/
  * [22] J.Nocedal, S.Wright, Numerical Optimization, 2006, Springer.
  * [23] S. Olver, A. Townsend, A fast and well-conditioned spectral method, SIAM Review 55 (2013), 462–489.
  * [24] pySDC project, a Python implementation of the spectral diferred correction,
https://parallel-in-time.org/pySDC/

  * [25] Risa/Asir, a Computer algebra system, http://www.openxm.org
  * [26] scipy.linalg.solve
https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.solve.html#scipy.linalg.solve

  * [27] scipy.integrate.solve_ivp,
https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html

  * [28] scipy.optimize.least_squares,
https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.least_squares.html

  * [29] cvxopt.solver.qp,
https://cvxopt.org/userguide/coneprog.html#quadratic-programming

  * [30] 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.
  * [31] L. N. Trefethen, Approximation Theory and Approximation Practice, 2020, SIAM.
  * [32] M.Winkel, R.Speck, D.Ruprecht, A high-order Boris integrator, Journal of Computational Physics 295 (2015), 456–474.