A Novel Third Order Numerical Method for Solving Volterra Integro-Differential
Equations
Sachin Bhalekar, Jayvant Patade111Corresponding author
Department of Mathematics, Shivaji University, Kolhapur - 416004, India.
Email<EMAIL_ADDRESS><EMAIL_ADDRESS>(Sachin Bhalekar),
<EMAIL_ADDRESS>(Jayvant Patade)

###### Abstract

In this paper we introduce a numerical method for solving nonlinear Volterra
integro-differential equations. In the first step, we apply implicit trapezium
rule to discretize the integral in given equation. Further, the Daftardar-
Gejji and Jafari technique (DJM) is used to find the unknown term on the right
side. We derive existence-uniqueness theorem for such equations by using
Lipschitz condition. We further present the error, convergence, stability and
bifurcation analysis of the proposed method. We solve various types of
equations using this method and compare the error with other numerical
methods. It is observed that our method is more efficient than other numerical
methods.

Keywords: Volterra Integro-differential equations, Trapezium rule, Daftardar-
Gejji and Jafari method, numerical solution, error, convergence, stability,
bifurcation.

##  1 Introduction

The equations involving derivative of a dependent variable are frequently used
to model the natural phenomena. The derivative of such function at a point
$t_{0}$ can be approximated by using its values in a small neighborhood of
$t_{0}$. Thus the derivative is a local operator which cannot model the memory
and hereditary properties involved in the real world problems. To overcome
this drawback, one has to include some nonlocal operator in the given
equation. Volterra [1, 2] suggested the use of an integral operator (which is
nonlocal) to model the problems in population dynamics. In his monograph [1],
he discussed the theory of integral, integro-differential and functional
equations. Existence-uniqueness, stability and applications of integro-
differential equations (IDE) is presented in a book by Lakshmikantham and Rao
[3]. Existence theory of nonlinear IDEs is also discussed in [4]. Recently,
Kostic [5] discussed the theory of abstract Volterra integro-differential
equations.

A-stable linear multi-step methods to solve Volterra IDEs (VIDE) are proposed
by Matthys in [6]. Brunner presented various numerical methods to solve VIDEs
in [7]. Day [8] used trapezoidal rule to devise a numerical method to solve
nonlinear VIDEs. Linz [9] derived fourth order numerical methods for such
equations. Runge-Kutta method, predictor-corrector method and explicit
multistep method are derived for IDEs by Wolfe and Phillips in [10]. Dehghan
and Salehi [11] developed a numerical scheme based on the moving least square
method for these equations. Singular IDEs in generalized Holder spaces are
solved by using collocation and mechanical quadrature methods in [12]. Saeedi
et al. [13] described the operational Tau method for solving nonlinear VIDEs
of the second kind. The stability of numerical methods for VIDEs is studied in
[14, 15]. A survey of numerical treatments of VIDEs is taken in [16].

The applications of integro-differential equations can be found in many fields
such as mechanics and electromagnetic theory [17], nuclear reactor [18, 3],
visco-elasticity [19, 20], heat conduction in materials with memory [21, 22]
and man-environment epidemics [23].

In the present work we utilize a decomposition method proposed by Daftardar-
Gejji and Jafari (DJM) [24] to generate a more accurate and faster numerical
method for solving nonlinear VIDEs.The paper is organized as follows:

The preliminaries are given in section 2. A new numerical method is presented
in Section 3. Analysis of this numerical method is given in Section 4. Section
5 deals with different types of illustrative examples. In Section 6 we provide
software package and conclusions are summarized in Section 7.

##  2 Preliminaries

###  2.1 Basic Definitions and Results:

In this section, we discuss some basic definitions and results[25, 3, 9].

###### Definition 2.1

Let $y_{j}$ be the approximation to the exact value $y(x_{j})$ obtained by a
given method with step-size h. Then a method is said to be convergent if and
only if

$\lim\limits_{h\rightarrow 0}\mid y(x_{j})-y_{j}\mid\rightarrow 0,\quad
j=1,2\cdots N.$ (2.1)

###### Definition 2.2

A method is said to be of order p if p is the largest number for which there
exists a finite constant C such that

$\mid y(x_{j})-y_{j}\mid\leq\textsc{C}h^{p},\quad j=1,2\cdots N.$ (2.2)

###### Theorem 2.1

Consider the integral equation

$y(t)=f(t)+\int_{0}^{x}K(x,t,y(t))dt,\quad y(x_{0})=y_{0}.$ (2.3)

Assume that
(i) f is continuous in $0\leq x\leq a$ ,
(ii) K(x, t, y) is continuous function for $0\leq t\leq x\leq a$ and
$\parallel y\parallel<\infty$,
(iii) K(x, t, y) satisfies a Lipschitz condition

$\parallel K(x,t,y_{1})-K(x,t,y_{2})\parallel\leq L\parallel
y_{1}-y_{2}\parallel$ (2.4)

for all $0\leq t\leq x\leq a$. Then the Eq. (2.3) has a unique solution.

###### Theorem 2.2

Assume that $f\in C[I\times\mathbb{R}^{n},\mathbb{R}^{n}]$, $K\in C[I\times
I\times\mathbb{R}^{n},\mathbb{R}^{n}]$ and $\int_{s}^{x}\mid K(t,s,y(s))\mid
dt\leq N$, for $x_{0}\leq s\leq x\leq x_{0}+a$, $y\in\Omega=\\{\phi\in
C[I,\mathbb{R}^{n}]:\phi(x_{0})=x_{0}\ and\ \mid\phi(x)-y_{0}\mid\leq b\\}$.
Then IVP (2.3) possesses at least one solution.

###  2.2 Existence and Uniqueness Theorem

In this subsection we prove existence and uniqueness theorem, which is
generalization of Theorem 2.1.

###### Theorem 2.3

Consider the Volterra integro-differential equation

$\displaystyle y^{\prime}(x)$ $\displaystyle=$ $\displaystyle
f(x,y(x))+\int_{x_{0}}^{x}K(x,t,y(t))dt,$ (2.5) $\displaystyle\textrm{with
initial condition}\quad y(x_{0})$ $\displaystyle=$ $\displaystyle y_{0}.$

Assume that $f$ and $K$ are continuous and satisfy Lipschitz condition

$\displaystyle\parallel f(x,y_{1})-f(x,y_{2})\parallel$ $\displaystyle\leq$
$\displaystyle L_{1}\parallel y_{1}-y_{2}\parallel$ (2.6)
$\displaystyle\parallel K(x,t,y_{1})-K(x,t,y_{2})\parallel$
$\displaystyle\leq$ $\displaystyle L_{2}\parallel y_{1}-y_{2}\parallel$ (2.7)

for every $\mid x-x_{0}\mid\leq a,\mid t-x_{0}\mid\leq a,\parallel
y_{1}\parallel<\infty,\parallel y_{2}\parallel<\infty\quad\textrm{and}\quad
a>0$. Then the Eq. (2.5) has unique solution.

Proof: Integrating Eq.(2.5) and using $y(x_{0})=y_{0}$, we get

$\displaystyle y(x)$ $\displaystyle=$ $\displaystyle
y_{0}+\int_{x_{0}}^{x}f(x,y(x))dx+\int_{x_{0}}^{x}\left(\int_{x_{0}}^{x}K(x,t,y(t))dt\right)dx$
$\displaystyle\Rightarrow y(x)$ $\displaystyle=$ $\displaystyle
y_{0}+\int_{x_{0}}^{x}\left(f(x,y(x))+\int_{x_{0}}^{x}K(x,t,y(t))dt\right)dx$
$\displaystyle\Rightarrow y(x)$ $\displaystyle=$ $\displaystyle
y_{0}+\int_{x_{0}}^{x}G(x,y(x))dx,$ (2.8) $\displaystyle\textrm{where}\quad
G(x,y(x))$ $\displaystyle=$ $\displaystyle
f(x,y(x))+\int_{x_{0}}^{x}K(x,t,y(t))dt.$

The Eq. (2.8) is of the form (2.3).
We have following observations
(i) $y_{0}$ is continuous because it is constant,
(ii) Kernel $G$ is continuous for $0\leq x\leq a$, because $f$ and $K$ are
continuous in the same domain,

$\displaystyle(iii)\parallel G(x,y_{1})-G(x,y_{2})\parallel$ $\displaystyle=$
$\displaystyle\parallel
f(x,y_{1}(x))+\int_{x_{0}}^{x}K(x,t,y_{1}(t))dt-f(x,y_{2}(x))$
$\displaystyle-\int_{x_{0}}^{x}K(x,t,y_{2}(t))dt\parallel$ $\displaystyle\leq$
$\displaystyle\parallel
f(x,y_{1}(x))-f(x,y_{2}(x))\parallel+\parallel\int_{x_{0}}^{x}K(x,t,y_{1}(t))dt$
$\displaystyle-\int_{x_{0}}^{x}K(x,t,y_{2}(t))dt\parallel$ $\displaystyle\leq$
$\displaystyle L_{1}\parallel y_{1}-y_{2}\parallel+aL_{2}\parallel
y_{1}-y_{2}\parallel\quad(\because\mid x-x_{0}\mid\leq a)$ $\displaystyle\leq$
$\displaystyle(L_{1}+aL_{2})\parallel y_{1}-y_{2}\parallel$

$\therefore G$ satisfy a Lipschitz condition.
$\Rightarrow$ All the conditions of the Theorem 2.1 are satisfies.
Hence the Eq.(2.5) has unique solution.

###  2.3 Daftardar-Gejji and Jafari Method

A new iterative method was introduced by Daftardar-Gejji and Jafari (DJM) [24]
in 2006 for solving nonlinear functional equations. The DJM has been used to
solve a variety of equations such as fractional differential equations [26],
partial differential equations [27], boundary value problems [28, 29],
evolution equations [30] and system of nonlinear functional equations [31].
The method is successfully employed to solve Newell-Whitehead-Segel equation
[32], fractional-order logistic equation [33] and some nonlinear dynamical
systems [34] also. Recently DJM has been used to generate new numerical
methods [35, 36, 37] for solving differential equations. In this section we
describe DJM which is very useful for solving the equations of the form

$u=g+L(u)+N(u),$ (2.9)

where $g$ is a given function, $L$ and $N$ are linear and nonlinear operators
respectively. DJM provides the solution to Eq.(2.9) in the form of series

$u=\sum_{i=0}^{\infty}u_{i}.$ (2.10)

where

$\displaystyle u_{0}$ $\displaystyle=$ $\displaystyle g,$ $\displaystyle
u_{m+1}$ $\displaystyle=$ $\displaystyle L(u_{m})+G_{m},\quad m=0,1,2,\cdots,$
(2.11) $\displaystyle\textrm{and}\quad G_{m}$ $\displaystyle=$ $\displaystyle
N\left(\sum_{j=0}^{m}u_{m}\right)-N\left(\sum_{j=0}^{m-1}u_{m}\right),m\geq
1.$

The $k$-term approximate solution is given by

$\displaystyle u=\sum_{i=0}^{k-1}u_{i}$ (2.12)

for suitable integer $k$.

The following convergence results for DJM are described in [38].

###### Theorem 2.4

If $N$ is $C^{(\infty)}$ in a neighborhood of $u_{0}$ and
$\left\|N^{(n)}(u_{0})\right\|\leq L$, for any $n$ and for some real $L>0$ and
$\left\|u_{i}\right\|\leq M<\frac{1}{e}$, $i=1,2,\cdots,$ then the series
$\sum_{n=0}^{\infty}G_{n}$ is absolutely convergent to $N$ and moreover,

$\left\|G_{n}\right\|\leq LM^{n}e^{n-1}(e-1),\quad n=1,2,\cdots.$

###### Theorem 2.5

If $N$ is $C^{(\infty)}$ and $\left\|N^{(n)}(u_{0})\right\|\leq M\leq e^{-1}$,
$\forall n$, then the series $\sum_{n=0}^{\infty}G_{n}$ is absolutely
convergent to $N$.

##  3 Numerical Method

In this section we present a numerical method based on DJM to solve Volterra
integro-differential equation.

Integrating Eq. (2.5) from $x=x_{j}$ to $x=x_{j}+h$, we get

$y({x_{j}+h})=y(x_{j})+\int_{x_{j}}^{x_{j}+h}f(x,y(x))dx+\int_{x_{j}}^{x_{j}+h}\int_{x_{0}}^{x}K(x,t,y(t))dtdx$
(3.1)

Applying trapezium formula [39] to evaluate integrals on right of Eq. (3.1),
we get

$\displaystyle y({x_{j}+h})$ $\displaystyle=$ $\displaystyle
y(x_{j})+\frac{h}{2}f(x_{j},y_{j})+\frac{h^{2}}{4}\left(K(x_{j},x_{0},y_{0})+K(x_{j},x_{j},y_{j})+K(x_{j+1},x_{0},y_{0})\right)$
(3.2)
$\displaystyle+\frac{h^{2}}{2}\left(\sum_{i=1}^{j-1}K(x_{j},x_{i},y_{i})+\sum_{i=1}^{j}K(x_{j+1},x_{i},y_{i})\right)+\frac{h}{2}f(x_{j+1},y_{j+1})$
$\displaystyle+\frac{h^{2}}{4}K(x_{j+1},x_{j+1},y_{j+1})+O(h^{3}).$

If $y_{j}$ is an approximation to $y(x_{j})$, then approximate solution is
given by

$\displaystyle y_{j+1}$ $\displaystyle=$ $\displaystyle
y_{j}+\frac{h}{2}f(x_{j},y_{j})+\frac{h^{2}}{4}\left(K(x_{j},x_{0},y_{0})+K(x_{j},x_{j},y_{j})+K(x_{j+1},x_{0},y_{0})\right)$
(3.3)
$\displaystyle+\frac{h^{2}}{2}\left(\sum_{i=1}^{j-1}K(x_{j},x_{i},y_{i})+\sum_{i=1}^{j}K(x_{j+1},x_{i},y_{i})\right)+\frac{h}{2}f(x_{j+1},y_{j+1})$
$\displaystyle+\frac{h^{2}}{4}K(x_{j+1},x_{j+1},y_{j+1}).$

Equation (3.3) is of the form (2.9), where

$\displaystyle u$ $\displaystyle=$ $\displaystyle y_{j+1},$ $\displaystyle g$
$\displaystyle=$ $\displaystyle
y_{j}+\frac{h}{2}f(x_{j},y_{j})+\frac{h^{2}}{4}\left(K(x_{j},x_{0},y_{0})+K(x_{j},x_{j},y_{j})+K(x_{j+1},x_{0},y_{0})\right)$
$\displaystyle+\frac{h^{2}}{2}\left(\sum_{i=1}^{j-1}K(x_{j},x_{i},y_{i})+\sum_{i=1}^{j}K(x_{j+1},x_{i},y_{i})\right)$
$\displaystyle N(u)$ $\displaystyle=$
$\displaystyle\frac{h}{2}f(x_{j+1},y_{j+1})+\frac{h^{2}}{4}K(x_{j+1},x_{j+1},y_{j+1}).$

Applying DJM to equation (3.3), we obtain 3-term solution as

$\displaystyle u$ $\displaystyle=$ $\displaystyle u_{0}+u_{1}+u_{2}$
$\displaystyle=$ $\displaystyle u_{0}+N(u_{0})+N(u_{0}+u_{1})-N(u_{0})$
$\displaystyle=$ $\displaystyle u_{0}+N(u_{0}+u_{1})$ $\displaystyle=$
$\displaystyle u_{0}+N(u_{0}+N(u_{0})).$

That is

$\displaystyle y_{j+1}$ $\displaystyle=$ $\displaystyle
y_{j}+\frac{h}{2}f(x_{j},y_{j})+\frac{h^{2}}{4}\left(K(x_{j},x_{0},y_{0})+K(x_{j},x_{j},y_{j})+K(x_{j+1},x_{0},y_{0})\right)$
(3.4)
$\displaystyle+\frac{h^{2}}{2}\left(\sum_{i=1}^{j-1}K(x_{j},x_{i},y_{i})+\sum_{i=1}^{j}K(x_{j+1},x_{i},y_{i})\right)$
$\displaystyle+N\left(y_{j}+\frac{h}{2}f(x_{j},y_{j})+\frac{h^{2}}{4}\left(K(x_{j},x_{0},y_{0})+K(x_{j},x_{j},y_{j})+K(x_{j+1},x_{0},y_{0})\right)\right.$
$\displaystyle\left.+\frac{h^{2}}{2}\left(\sum_{i=1}^{j-1}K(x_{j},x_{i},y_{i})+\sum_{i=1}^{j}K(x_{j+1},x_{i},y_{i})\right)\right)+N\left(y_{j}+\frac{h}{2}f(x_{j},y_{j})\right.$
$\displaystyle\left.+\frac{h^{2}}{4}\left(K(x_{j},x_{0},y_{0})+K(x_{j},x_{j},y_{j})+K(x_{j+1},x_{0},y_{0})\right)\right.$
$\displaystyle\left.+\frac{h^{2}}{2}\left(\sum_{i=1}^{j-1}K(x_{j},x_{i},y_{i})+\sum_{i=1}^{j}K(x_{j+1},x_{i},y_{i})\right)\right).$

or

$\displaystyle y_{j+1}$ $\displaystyle=$ $\displaystyle
y_{j}+\frac{h}{2}f(x_{j},y_{j})+\frac{h^{2}}{4}\left(K(x_{j},x_{0},y_{0})+K(x_{j},x_{j},y_{j})+K(x_{j+1},x_{0},y_{0})\right)$
(3.5)
$\displaystyle+\frac{h^{2}}{2}\left(\sum_{i=1}^{j-1}K(x_{j},x_{i},y_{i})+\sum_{i=1}^{j}K(x_{j+1},x_{i},y_{i})\right)$
$\displaystyle+\frac{h}{2}f\left(x_{j+1},y_{j}+\frac{h}{2}f(x_{j},y_{j})+\frac{h^{2}}{4}\left(K(x_{j},x_{0},y_{0})+K(x_{j},x_{j},y_{j})+K(x_{j+1},x_{0},y_{0})\right)\right.$
$\displaystyle\left.+\frac{h^{2}}{2}\left(\sum_{i=1}^{j-1}K(x_{j},x_{i},y_{i})+\sum_{i=1}^{j}K(x_{j+1},x_{i},y_{i})\right)+\frac{h}{2}f\left(x_{j+1},y_{j}+\frac{h}{2}f(x_{j},y_{j})\right.\right.$
$\displaystyle\left.\left.+\frac{h^{2}}{4}\left(K(x_{j},x_{0},y_{0})+K(x_{j},x_{j},y_{j})+K(x_{j+1},x_{0},y_{0})\right)\right.\right.$
$\displaystyle\left.\left.+\frac{h^{2}}{2}\left(\sum_{i=1}^{j-1}K(x_{j},x_{i},y_{i})+\sum_{i=1}^{j}K(x_{j+1},x_{i},y_{i})\right)\right)\right.$
$\displaystyle\left.+\frac{h^{2}}{4}K\left(x_{j+1},x_{j+1},y_{j}+\frac{h}{2}f(x_{j},y_{j})+\frac{h^{2}}{4}\left(K(x_{j},x_{0},y_{0})+K(x_{j},x_{j},y_{j})+K(x_{j+1},x_{0},y_{0})\right)\right.\right.$
$\displaystyle\left.\left.+\frac{h^{2}}{2}\left(\sum_{i=1}^{j-1}K(x_{j},x_{i},y_{i})+\sum_{i=1}^{j}K(x_{j+1},x_{i},y_{i})\right)\right)\right)$
$\displaystyle+\frac{h^{2}}{4}K\left(x_{j+1},x_{j+1},y_{j}+\frac{h}{2}f(x_{j},y_{j})+\frac{h^{2}}{4}\left(K(x_{j},x_{0},y_{0})+K(x_{j},x_{j},y_{j})+K(x_{j+1},x_{0},y_{0})\right)\right.$
$\displaystyle\left.+\frac{h^{2}}{2}\left(\sum_{i=1}^{j-1}K(x_{j},x_{i},y_{i})+\sum_{i=1}^{j}K(x_{j+1},x_{i},y_{i})\right)+\frac{h}{2}f\left(x_{j+1},y_{j}+\frac{h}{2}f(x_{j},y_{j})\right.\right.$
$\displaystyle\left.\left.+\frac{h^{2}}{4}\left(K(x_{j},x_{0},y_{0})+K(x_{j},x_{j},y_{j})+K(x_{j+1},x_{0},y_{0})\right)\right.\right.$
$\displaystyle\left.\left.+\frac{h^{2}}{2}\left(\sum_{i=1}^{j-1}K(x_{j},x_{i},y_{i})+\sum_{i=1}^{j}K(x_{j+1},x_{i},y_{i})\right)\right)\right.$
$\displaystyle\left.+\frac{h^{2}}{4}K\left(x_{j+1},x_{j+1},y_{j}+\frac{h}{2}f(x_{j},y_{j})+\frac{h^{2}}{4}\left(K(x_{j},x_{0},y_{0})+K(x_{j},x_{j},y_{j})+K(x_{j+1},x_{0},y_{0})\right)\right.\right.$
$\displaystyle\left.\left.+\frac{h^{2}}{2}\left(\sum_{i=1}^{j-1}K(x_{j},x_{i},y_{i})+\sum_{i=1}^{j}K(x_{j+1},x_{i},y_{i})\right)\right)\right)$

If we set

$\displaystyle M_{1}$ $\displaystyle=$ $\displaystyle
y_{j}+\frac{h}{2}f(x_{j},y_{j})+\frac{h^{2}}{4}\left(K(x_{j},x_{0},y_{0})+K(x_{j},x_{j},y_{j})+K(x_{j+1},x_{0},y_{0})\right)$
(3.6)
$\displaystyle+\frac{h^{2}}{2}\left(\sum_{i=1}^{j-1}K(x_{j},x_{i},y_{i})+\sum_{i=1}^{j}K(x_{j+1},x_{i},y_{i})\right),$
$\displaystyle M_{2}$ $\displaystyle=$ $\displaystyle
M_{1}+\frac{h}{2}f(x_{j+1},M_{1})+\frac{h^{2}}{4}K(x_{j+1},x_{j+1},M_{1})$
(3.7)

then equation (3.5) becomes

$\displaystyle
y_{j+1}=M_{1}+\frac{h}{2}f(x_{j+1},M_{2})+\frac{h^{2}}{4}K(x_{j+1},x_{j+1},M_{2}).$
(3.8)

##  4 Analysis of Numerical Method

###  4.1 Error Analysis

###### Theorem 4.1

The numerical method (3.8) is of third order.

Proof: Suppose $y_{j+1}$ is an approximation to $y(x_{j+1})$. By using
Eq.(3.2) and (3.8), we obtain

$\displaystyle\mid y(x_{j+1})-y_{j+1}\mid$ $\displaystyle=$
$\displaystyle\mid\frac{h}{2}f(x_{j+1},y_{j+1})+\frac{h^{2}}{4}K(x_{j+1},x_{j+1},y_{j+1})+O(h^{3})$
$\displaystyle-\frac{h}{2}f(x_{j+1},M_{2})-\frac{h^{2}}{4}K(x_{j+1},x_{j+1},M_{2})\mid$
$\displaystyle\leq$ $\displaystyle\frac{h}{2}\mid
f(x_{j+1},y_{j+1})-f(x_{j+1},M_{2})\mid+\frac{h^{2}}{4}\mid
K(x_{j+1},x_{j+1},y_{j+1})$
$\displaystyle-K(x_{j+1},x_{j+1},M_{2})\mid+O(h^{3})$ $\displaystyle\leq$
$\displaystyle\left(\frac{h}{2}L_{1}+\frac{h^{2}}{4}L_{2}\right)\mid
y_{j+1}-M_{2}\mid+O(h^{3}).$

Using Eq.(3.7) and (3.8), we get

$\displaystyle\mid y(x_{j+1})-y_{j+1}\mid$ $\displaystyle\leq$
$\displaystyle\left(\frac{h}{2}L_{1}+\frac{h^{2}}{4}L_{2}\right)\mid
M_{1}+\frac{h}{2}f(x_{j+1},M_{2})+\frac{h^{2}}{4}K(x_{j+1},x_{j+1},M_{2})$
$\displaystyle-
M_{1}-\frac{h}{2}f(x_{j+1},M_{1})-\frac{h^{2}}{4}K(x_{j+1},x_{j+1},M_{1})\mid+O(h^{3})$
$\displaystyle\leq$
$\displaystyle\left(\frac{h}{2}L_{1}+\frac{h^{2}}{4}L_{2}\right)\left(\mid\frac{h}{2}f(x_{j+1},M_{2})-\frac{h}{2}f(x_{j+1},M_{1})\mid\right.$
$\displaystyle\left.+\mid\frac{h^{2}}{4}K(x_{j+1},x_{j+1},M_{2})-\frac{h^{2}}{4}K(x_{j+1},x_{j+1},M_{1})\mid\right)+O(h^{3}).$

Using Eq.(3.6) and (3.7) ,we get

$\displaystyle\mid y(x_{j+1})-y_{j+1}\mid$ $\displaystyle\leq$
$\displaystyle\left(\frac{h}{2}L_{1}+\frac{h^{2}}{4}L_{2}\right)^{2}\mid\frac{h}{2}f(x_{j+1},M_{1})+\frac{h^{2}}{4}K(x_{j+1},x_{j+1},M_{2})+O(h^{3})$
$\displaystyle\leq$ $\displaystyle
h^{3}\left(\left(\frac{1}{2}L_{1}+\frac{h}{4}L_{2}\right)^{2}\mid\frac{1}{2}f(x_{j+1},M_{1})+\frac{h}{4}K(x_{j+1},x_{j+1},M_{2})\mid\right)+O(h^{2}).$

$\Rightarrow$ The numerical method (3.8) is of third order.

###### Corollary 1

The numerical method (3.8) is convergent.

Proof: By Theorem 4.1 and definition (2.1) the numerical method (3.8) is
convergent.

###  4.2 Stability Analysis of Numerical Method

Consider the test equation

$y^{\prime}=\alpha y+\beta\int_{x_{0}}^{x}y(t)dt.$ (4.1)

Applying numerical method (3.8) to this equation, we get

$y_{j+1}=M_{1}+\frac{h\alpha}{2}M_{2}+\frac{h^{2}\beta}{4}M_{2},$ (4.2)

where

$\displaystyle M_{1}$ $\displaystyle=$
$\displaystyle\frac{h^{2}\beta}{2}y_{0}+h^{2}\beta\sum_{i=1}^{j-1}y_{i}+\left(1+\frac{h\alpha}{2}+\frac{3h^{2}\beta}{4}\right)y_{j},$
(4.3) $\displaystyle M_{2}$ $\displaystyle=$
$\displaystyle\left(1+\frac{h\alpha}{2}+\frac{h^{2}\beta}{4}\right)M_{1}.$
(4.4)

If we set $u=h\alpha$ and $v=h^{2}\beta$ then Eq. (4.2) can be written as

$\displaystyle y_{j+1}$ $\displaystyle=$
$\displaystyle\left(\frac{v}{2}+\frac{uv}{4}+\frac{u^{2}v}{8}+\frac{v^{2}}{8}+\frac{uv^{2}}{8}+\frac{v^{3}}{32}\right)y_{0}$
(4.5)
$\displaystyle+\left(v+\frac{uv}{2}+\frac{u^{2}v}{4}+\frac{v^{2}}{4}+\frac{uv^{2}}{4}+\frac{v^{3}}{16}\right)\sum_{i=1}^{j-1}y_{i}$
$\displaystyle+\left(1+u+\frac{u^{2}}{2}+\frac{u^{3}}{8}+v+\frac{3uv}{4}+\frac{5u^{2}v}{16}+\frac{v^{2}}{4}+\frac{7uv^{2}}{32}+\frac{3v^{3}}{64}\right)y_{j}.$

Simplifying Eq.(4.5), we get

$y_{j}=b_{1}y_{j-1}+b_{2}y_{j-2},$ (4.6)

where

$\displaystyle b_{1}$ $\displaystyle=$ $\displaystyle
2+u+\frac{u^{2}}{2}+\frac{u^{3}}{8}+v+\frac{3uv}{4}+\frac{5u^{2}v}{16}+\frac{v^{2}}{4}+\frac{7uv^{2}}{32}+\frac{3v^{3}}{64},$
(4.7) $\displaystyle\textrm{and}\quad b_{2}$ $\displaystyle=$
$\displaystyle-1-u-\frac{u^{2}}{2}-\frac{u^{3}}{8}-\frac{uv}{4}-\frac{u^{2}v}{16}+\frac{uv^{2}}{32}+\frac{v^{3}}{64}.$
(4.8)

The characteristic equation [40] of the difference Eq. (4.6) is

$r^{2}-b_{1}r-b_{2}=0.$ (4.9)

The characteristic roots are given by

$r_{1}=\frac{b_{1}+\sqrt{b_{1}^{2}+4b_{2}}}{2},\quad
r_{2}=\frac{b_{1}-\sqrt{b_{1}^{2}+4b_{2}}}{2}.$ (4.10)

The zero solution of system (4.6) is asymptotically stable if $\mid
r_{1}\mid<1$ and $\mid r_{2}\mid<1$. Thus the stability region of numerical
method (3.8) is given by $\mid r_{1}\mid<1$ and $\mid r_{2}\mid<1$ as shown in
Fig.1.

---
Fig.1 Stability region for numerical method (3.8)

###  4.3 Bifurcations Analysis

We discuss bifurcation analysis for the following test equation given in [41]

$y^{\prime}=-\int_{x_{0}}^{x}e^{-\alpha(x-t)}y(t)dt.$ (4.11)

As discussed in [41], Eq.(4.11) with convolution kernel and with exponentially
fading memory is practically important. The qualitative behavior of the
Eq.(4.11) is described as below [41]:
Table 1: Bifurcation in IDE (4.11).
Region Qualitative behavior $\alpha\geq 2$ Solutions are converging to zero
without oscillation $0<\alpha<2$ Solutions are damped oscillations converging
to zero $-2<\alpha<0$ Divergent solutions exhibiting unbounded oscillations
$\alpha\leq-2$ Solutions become unbounded without oscillations

Now we apply the numerical scheme (3.8) to Eq.(4.11) and discuss the
bifurcations.

$y_{j+1}=-a_{0}e^{-\alpha hj}y_{0}-a_{1}\sum_{i=1}^{j-1}e^{-\alpha
h(j-1)}y_{i}+a_{2}y_{j},$ (4.12)

where

$\displaystyle a_{0}$ $\displaystyle=$
$\displaystyle\frac{h^{2}}{4}\left(1-\frac{h^{2}}{4}+\frac{h^{4}}{16}\right)\left(1+e^{-\alpha
h}\right),$ $\displaystyle a_{1}$ $\displaystyle=$
$\displaystyle\frac{h^{2}}{2}\left(1-\frac{h^{2}}{4}+\frac{h^{4}}{16}\right)\left(1+e^{-\alpha
h}\right),$ $\displaystyle a_{2}$ $\displaystyle=$
$\displaystyle\left(1-\frac{h^{2}}{4}+\frac{h^{4}}{16}\right)\left(1-\frac{h^{2}}{2}+\frac{h^{2}}{4}e^{-\alpha
h}\right).$

Eq. (4.11) can be written as

$y_{j+2}-b_{1}y_{j+1}+b_{2}y_{j}=0,$ (4.13)

where

$\displaystyle b_{1}$ $\displaystyle=$ $\displaystyle e^{-\alpha h}+a_{2},$
$\displaystyle b_{2}$ $\displaystyle=$ $\displaystyle(a_{1}+a_{2})e^{-\alpha
h}.$

The characteristic equation of the difference Eq. (4.13) is

$r^{2}-b_{1}r+b_{2}=0.$ (4.14)

The roots of this equation are

$r=\frac{b_{1}\pm\sqrt{b_{1}^{2}-4b_{2}}}{2}.$ (4.15)

###### Theorem 4.2

[40] All solutions of the difference equation

$y_{j+2}+p_{1}y_{j+1}+p_{2}y_{j}=0,$ (4.16)

converge to zero (i.e.,the zero solution is asymptotically stable) if and only
if $max\\{|r_{1}|,|r_{2}|\\}<1$, where $r_{1},r_{2}$ are the roots of the
characteristic equation $r^{2}+p_{1}r+p_{2}=0$ of Eq. (4.16).

Define

$\displaystyle\alpha_{0}$ $\displaystyle=$
$\displaystyle\frac{1}{h}\text{Log}\left(1+\frac{h^{6}}{64}\right),$
$\displaystyle\alpha_{1}$ $\displaystyle=$
$\displaystyle\frac{1}{h}\text{Log}\left(\frac{2\left(128+160h^{2}-32h^{4}+8h^{6}-h^{8}+8h\sqrt{1024+8h^{6}-2h^{8}}\right)}{\left(-4+h^{2}\right)^{2}\left(16-4h^{2}+h^{4}\right)}\right)\quad\textrm{and}$
$\displaystyle\alpha_{2}$ $\displaystyle=$
$\displaystyle\frac{1}{h}\text{Log}\left(-\frac{2\left(-128-160h^{2}+32h^{4}-8h^{6}+h^{8}+8h\sqrt{1024+8h^{6}-2h^{8}}\right)}{\left(-4+h^{2}\right)^{2}\left(16-4h^{2}+h^{4}\right)}\right).$

Using Theorem 4.2 we have the following observations:
Case 1: The characteristics roots of the Eq. (4.15) are real and distinct

$\displaystyle r_{1}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(1+e^{-h\alpha}-\frac{h^{2}}{2}-\frac{1}{2}e^{-h\alpha}h^{2}+\frac{h^{4}}{8}+\frac{1}{8}e^{-h\alpha}h^{4}-\frac{h^{6}}{64}-\frac{1}{32}e^{-h\alpha}h^{6}+\right.$
$\displaystyle\left.\sqrt{-4e^{-h\alpha}\left(1+\frac{h^{6}}{64}\right)+\left(1+e^{-h\alpha}-\frac{h^{2}}{2}-\frac{1}{2}e^{-h\alpha}h^{2}+\frac{h^{4}}{8}+\frac{1}{8}e^{-h\alpha}h^{4}-\frac{h^{6}}{64}-\frac{1}{32}e^{-h\alpha}h^{6}\right)^{2}}\right)$
$\displaystyle r_{2}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(1+e^{-h\alpha}-\frac{h^{2}}{2}-\frac{1}{2}e^{-h\alpha}h^{2}+\frac{h^{4}}{8}+\frac{1}{8}e^{-h\alpha}h^{4}-\frac{h^{6}}{64}-\frac{1}{32}e^{-h\alpha}h^{6}-\right.$
$\displaystyle\left.\sqrt{-4e^{-h\alpha}\left(1+\frac{h^{6}}{64}\right)+\left(1+e^{-h\alpha}-\frac{h^{2}}{2}-\frac{1}{2}e^{-h\alpha}h^{2}+\frac{h^{4}}{8}+\frac{1}{8}e^{-h\alpha}h^{4}-\frac{h^{6}}{64}-\frac{1}{32}e^{-h\alpha}h^{6}\right)^{2}}\right),$

when
Case 1.1: $\alpha>\alpha_{1}$.
In this case the characteristic roots are of magnitude less than unity (cf.
Region I in Fig.2). Therefore, all the solutions of difference Eq. (4.13)
converge to zero without oscillations.
Case 1.2: $\alpha<\alpha_{2}$.
In this case the characteristic roots are of magnitude greater than unity.
This is shown in Fig.2 as Region IV. In this region, all the solutions will
tend to infinity without oscillations.
Case 2: The characteristics roots are real and equal

$r=r_{1}=r_{2}=\frac{1}{2}\left(1+e^{-h\alpha}-\frac{h^{2}}{2}-\frac{1}{2}e^{-h\alpha}h^{2}+\frac{h^{4}}{8}+\frac{1}{8}e^{-h\alpha}h^{4}-\frac{h^{6}}{64}-\frac{1}{32}e^{-h\alpha}h^{6}\right),$

when
Case 2.1: $\alpha=\alpha_{1}$.
In this case the characteristic root $r\in(0,1)$. Hence the solutions of Eq.
(4.13) tend to zero without oscillations as in case 1.1. This case is shown in
Fig.2 as boundary line of Region I and Region II.
Case 2.2: $\alpha=\alpha_{2}$.
In this case the characteristic root $r>1$ and the solutions are unbounded
(cf. boundary line of Region III and Region IV in Fig.2).
Case 3: The characteristics roots are complex

$\displaystyle r_{1}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(1+e^{-h\alpha}-\frac{h^{2}}{2}-\frac{1}{2}e^{-h\alpha}h^{2}+\frac{h^{4}}{8}+\frac{1}{8}e^{-h\alpha}h^{4}-\frac{h^{6}}{64}-\frac{1}{32}e^{-h\alpha}h^{6}\right.$
$\displaystyle\left.+i\sqrt{4e^{-h\alpha}\left(1+\frac{h^{6}}{64}\right)-\left(1+e^{-h\alpha}-\frac{h^{2}}{2}-\frac{1}{2}e^{-h\alpha}h^{2}+\frac{h^{4}}{8}+\frac{1}{8}e^{-h\alpha}h^{4}-\frac{h^{6}}{64}-\frac{1}{32}e^{-h\alpha}h^{6}\right)^{2}}\right)$
$\displaystyle r_{2}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(1+e^{-h\alpha}-\frac{h^{2}}{2}-\frac{1}{2}e^{-h\alpha}h^{2}+\frac{h^{4}}{8}+\frac{1}{8}e^{-h\alpha}h^{4}-\frac{h^{6}}{64}-\frac{1}{32}e^{-h\alpha}h^{6}\right.$
$\displaystyle\left.-i\sqrt{4e^{-h\alpha}\left(1+\frac{h^{6}}{64}\right)-\left(1+e^{-h\alpha}-\frac{h^{2}}{2}-\frac{1}{2}e^{-h\alpha}h^{2}+\frac{h^{4}}{8}+\frac{1}{8}e^{-h\alpha}h^{4}-\frac{h^{6}}{64}-\frac{1}{32}e^{-h\alpha}h^{6}\right)^{2}}\right),$

when $\alpha_{2}<\alpha<\alpha_{1}$.
Note that $\mid r_{1}\mid=\mid r_{2}\mid$.
Case 3.1: If $\alpha_{0}<\alpha<\alpha_{1}$ then $\mid r_{1}\mid<1$. The
solutions of Eq. (4.13) are oscillatory and converge to zero (Region II in
Fig.2).
Case 3.2: If $\alpha=\alpha_{0}$ then $\mid r_{1}\mid=1$ and bifurcation
occurs. This is shown by the boundary line of Region II and Region III in
Fig.2.
Case 3.3: If $\alpha_{2}<\alpha<\alpha_{0}$
then $\mid r_{1}\mid>1$. This implies that the solutions of Eq. (4.13) are
oscillatory and diverge to infinity (Region III in Fig.2).

---
Fig.2: Bifurcation digram for Eq. (4.13).

Note: It is clear from Table 1 and Fig.2 that the bifurcations in numerical
scheme (4.13) coincide with those in IDE (4.11) for $h=0$.

##  5 Illustrative examples

###### Example 5.1

[8, 10] Consider the Volterra integro-differential equation

$y^{\prime}(x)=1+2x-y(x)+\int_{x_{0}}^{x}x(1+2x)e^{t(x-t)}y(t)dt,\quad
y(0)=1.$ (5.1)

The exact solution of Eq.(5.1) is $y(x)=e^{x^{2}}$. We compare errors obtained
using different numerical methods given in [8, 10], viz. third order Runge-
Kutta method, predictor-corrector method, multi-step method and Day method
with our method for $h=0.1$ and $h=0.025$ in Table 2 and Table 3 respectively.
It is observed that the error in our method is very less as compared with
other numerical methods.
Note that the error values in Table 2 and Table 3 for all other methods are
taken from [8, 10].

Table 2: Comparison of numerical results of example 6.1 with h=0.1
x Our method Runge-Kutta Predictor- Multi-step Day method $O(h^{3})$
$O(h^{3})$ Corrector $O(h^{3})$ $O(h^{3})$ 0.1 $2.3\times 10^{-5}$ $5.5\times
10^{-4}$ $5.0\times 10^{-5}$ $5.0\times 10^{-5}$ $5.4\times 10^{-4}$ 0.2
$1.2\times 10^{-4}$ $1.2\times 10^{-3}$ $5.2\times 10^{-5}$ $8.1\times
10^{-4}$ $1.1\times 10^{-3}$ 0.3 $3.1\times 10^{-4}$ $1.9\times 10^{-3}$
$2.7\times 10^{-4}$ $7.5\times 10^{-4}$ $1.9\times 10^{-3}$ 0.4 $6.2\times
10^{-3}$ $2.8\times 10^{-3}$ $6.6\times 10^{-4}$ $2.1\times 10^{-3}$
$2.8\times 10^{-3}$ 0.5 $1.0\times 10^{-3}$ $3.9\times 10^{-3}$ $1.2\times
10^{-3}$ $2.5\times 10^{-3}$ $4.0\times 10^{-3}$ 0.6 $1.7\times 10^{-3}$
$5.8\times 10^{-3}$ $2.0\times 10^{-3}$ $4.7\times 10^{-3}$ $5.6\times
10^{-3}$ 0.7 $2.8\times 10^{-3}$ $7.3\times 10^{-3}$ $3.1\times 10^{-3}$
$6.0\times 10^{-3}$ $7.9\times 10^{-3}$ 0.8 $4.4\times 10^{-3}$ $1.0\times
10^{-2}$ $4.8\times 10^{-3}$ $1.0\times 10^{-2}$ $1.1\times 10^{-2}$ 0.9
$6.8\times 10^{-3}$ $1.4\times 10^{-2}$ $7.3\times 10^{-3}$ $1.4\times
10^{-2}$ $1.5\times 10^{-2}$ 1 $1.0\times 10^{-2}$ $2.0\times 10^{-2}$
$1.1\times 10^{-2}$ $2.1\times 10^{-2}$ $2.2\times 10^{-2}$

Table 3: Comparison of numerical results of example 6.1 with h=0.025
x Our method Runge-Kutta Predictor- Multi-step Day method $O(h^{3})$
$O(h^{3})$ Corrector $O(h^{3})$ $O(h^{3})$ 0.1 $2.1\times 10^{-6}$ $3.3\times
10^{-5}$ $2.2\times 10^{-6}$ $7.7\times 10^{-6}$ $3.3\times 10^{-5}$ 0.2
$9.0\times 10^{-6}$ $7.1\times 10^{-5}$ $1.1\times 10^{-5}$ $2.6\times
10^{-5}$ $7.1\times 10^{-5}$ 0.3 $2.1\times 10^{-5}$ $1.2\times 10^{-4}$
$2.7\times 10^{-5}$ $5.7\times 10^{-5}$ $1.1\times 10^{-4}$ 0.4 $4.1\times
10^{-5}$ $1.7\times 10^{-4}$ $5.0\times 10^{-5}$ $1.0\times 10^{-4}$
$1.7\times 10^{-4}$ 0.5 $7.0\times 10^{-5}$ $2.4\times 10^{-4}$ $8.4\times
10^{-5}$ $1.7\times 10^{-4}$ $2.5\times 10^{-4}$ 0.6 $1.7\times 10^{-4}$
$3.3\times 10^{-4}$ $1.3\times 10^{-4}$ $2.6\times 10^{-4}$ $3.5\times
10^{-4}$ 0.7 $1.8\times 10^{-4}$ $4.5\times 10^{-4}$ $2.0\times 10^{-4}$
$4.0\times 10^{-4}$ $5.0\times 10^{-4}$ 0.8 $2.8\times 10^{-4}$ $6.3\times
10^{-4}$ $3.1\times 10^{-4}$ $5.9\times 10^{-4}$ $7.0\times 10^{-4}$ 0.9
$4.3\times 10^{-4}$ $8.7\times 10^{-4}$ $4.6\times 10^{-4}$ $8.8\times
10^{-4}$ $1.0\times 10^{-3}$ 1 $6.6\times 10^{-4}$ $1.2\times 10^{-3}$
$7.0\times 10^{-4}$ $1.3\times 10^{-3}$ $1.4\times 10^{-3}$

###### Example 5.2

[11] Consider the nonlinear Volterra integro-differential equation

$y^{\prime}(x)=2x-\frac{1}{2}sin(x^{4})+\int_{0}^{x}x^{2}tcos(x^{2}y(t))dt,\quad
y(0)=0.$ (5.2)

The exact solution of Eq.(5.2) is $y(x)=x^{2}$. We compare the maximum
absolute errors obtained using meshless method given in [11] and CPU times
used for different values of n in the interval $[0,1]$. It is observed that
the error in our method is very less as compared with linear meshless method.
Further, our method is more time efficient than both linear and quadratic
meshless method.

Table 4: Maximum absolute errors and CPU times used for different values of n.
n Our method Time Linear Meshless Time Quadratic Meshless Time method method 5
$1.49\times 10^{-3}$ 0.00 $2.83\times 10^{-3}$ 0.42 $2.21\times 10^{-4}$ 0.58
9 $4.54\times 10^{-4}$ 0.00 $7.22\times 10^{-4}$ 0.68 $3.41\times 10^{-5}$
0.79 17 $1.26\times 10^{-4}$ 0.00 $2.25\times 10^{-4}$ 0.76 $9.98\times
10^{-6}$ 1.16 33 $3.35\times 10^{-5}$ 0.01 $6.93\times 10^{-5}$ 0.87
$3.26\times 10^{-6}$ 1.85 65 $8.26\times 10^{-6}$ 0.04 $2.19\times 10^{-5}$
0.95 $1.87\times 10^{-6}$ 2.15 129 $2.19\times 10^{-6}$ 0.15 $6.58\times
10^{-6}$ 1.40 $1.28\times 10^{-6}$ 2.91

###### Example 5.3

Consider the VIDE [11]

$y^{\prime}(x)=1-\frac{x}{2}+\frac{xe^{-x^{2}}}{2}+\int_{0}^{x}xte^{-y^{2}(t)}dt,\quad
y(0)=0.$ (5.3)

The exact solution of Eq.(5.3) is $y(x)=x$. We compare the maximum absolute
errors obtained using meshless method given in [11] and CPU time used for
different values of n in the interval $[0,1]$. It is observed that our method
is more accurate as compared with linear meshless and more time efficient than
both linear and quadratic meshless method.
Table 5: Maximum absolute errors and CPU times used for different values of n.
n Our method Time Linear Meshless Time Quadratic Meshless Time method method 5
$8.96\times 10^{-3}$ 0.00 $5.84\times 10^{-3}$ 0.45 $1.96\times 10^{-4}$ 0.46
9 $2.33\times 10^{-4}$ 0.00 $1.75\times 10^{-3}$ 0.57 $2.42\times 10^{-5}$
0.58 17 $6.13\times 10^{-5}$ 0.00 $4.88\times 10^{-4}$ 0.65 $3.48\times
10^{-6}$ 0.73 33 $1.59\times 10^{-5}$ 0.01 $1.30\times 10^{-4}$ 0.72
$6.55\times 10^{-7}$ 0.87 65 $4.10\times 10^{-6}$ 0.04 $3.21\times 10^{-5}$
1.68 $5.75\times 10^{-7}$ 1.66 129 $1.04\times 10^{-6}$ 0.17 $8.19\times
10^{-6}$ 4.83 $6.05\times 10^{-7}$ 5.61

We compare our solution with exact solution for $h=0.01$ in Fig.3. It is
observed that our solution coincides with exact solution.

---
Fig.3: Comparison of solutions of Eq.(5.3) for $h=0.01$.

###### Example 5.4

Consider the VIDE

$y^{\prime}(x)=1+\int_{0}^{x}e^{-t}y^{2}(t)dt,\quad y(0)=1.$ (5.4)

The exact solution of Eq.(5.4) is $y(x)=e^{x}$. We compare our solution with
exact solution for $h=0.01$ in Fig.4. It is observed that our solution is well
in agreement with exact solution.

---
Fig.4: Comparison of solutions of Eq.(5.4) for $h=0.01$.

##  6 Software Package

In this section we provide software package based on Mathematica-$10$ to solve
VIDEs using numerical method.

This software is used for solving VIDEs of the form (5.3). One has to provide
the value of $f(x,y)$ and $K(x,t,y)$ in first two windows, initial condition
$y_{0}$, step size $h$ and number of steps $n$ in third, fourth and fifth
windows respectively. The last window gives required solution curve. We solve
Eq.(5.2) using this software package as shown in Fig.5.

---
Fig.5: Software Package

##  7 Conclusions

A new method proposed in this paper is a combination of trapezium rule and
DJM. The DJM is used to find the approximation to the unknown term on the
right side of the difference formula obtained by using trapezium rule. This
combination is proved to be very efficient. The error in this third order
method is very less as compared to other similar methods. However, the time
taken for the computation was very less. The stability and bifurcation
analysis also show the power of this method. An ample number of examples show
the wide range of applicability of this method.

Acknowledgements:
S. Bhalekar acknowledges CSIR, New Delhi for funding through Research Project
[25(0245)/15/EMR-II].

## References

  * [1] V. Volterra, Theory of Functionals and of Integral and Integro-differential Equations, Dover, New York, 1959.
  * [2] V. Volterra, Lecons sur la theorie mathematique de la luttle pour la vie Gauthier - Villars, Paris, 1931.
  * [3] V. Lakshmikantham and M. R. M. Rao, Theory of Integro-Differential Equations, Gordon and Breach Science Publishers, Amsterdam, 1995.
  * [4] D. O’Regan and M. Meehan, Existence Theory for Nonlinear Integral and Integro-differential Equations, Kluwer Academic, Dordrecht, 1998.
  * [5] M. Kostic, Abstract Volterra Integro-Differential Equations, CRC Press, Boca Raton, 2015.
  * [6] J. Matthys, A-Stable Linear Multistep Methods for Volterra Integro-Differential Equations, Numer. Math. 27, (1976) 85–94.
  * [7] H. Brunner, High-order methods for the numerical solution of Volterra integro-differential equations, J. Comput. Appl. Math., 15 (1986) 301–309.
  * [8] J. T. Day, Note on the numerical solution of integro-differential equations. The Computer Journal, 9(4), (1967) 394–395.
  * [9] P. Linz, A method for solving nonlinear Volterra integral equations of the second kind, Mathematics of Computation, 23(107) (1969) 595–599.
  * [10] M. A. Wolfe, and G. M. Phillips, Some methods for the solution of non-singular Volterra integro-differential equations, The Computer Journal, 11(3), (1968) 334–336.
  * [11] M. Dehghan, R. Salehi, The numerical solution of the non-linear integro-differential equations based on the meshless method, J. Comput. Appl. Math., 2012, 236(9), 2367-2377.
  * [12] I. Caraus, F. Al Faqih Approximate solution of singular integro-differential equations in generalized Holder spaces, Numer. Algor., 45(1-4) (2007) 205–215.
  * [13] L. Saeedi, A. Tari and S. H. M. Masuleh, Numerical solution of a class of the nonlinear Volterra integro-differential equations, J. Appl. Math. & Informatics, 31 (2013) 65–77.
  * [14] H. Brunner and J. D. Lambert, Stability of numerical methods for Volterra integro-differential equations, Computing 12 (1974) 75–89.
  * [15] C. T. H. Baker, A. Makroglou and E. Short, Regions of stability in the numerical treatment of Volterra integro-differential equations, SIAM J. Numer. Anal., 16(6) (1979) 890–910.
  * [16] H. Brunner, A survey of recent advances in the numerical treatment of Volterra integral and integro-differential equations, J. Comput. Appl. Math., 8 (1982) 213–229.
  * [17] F. Bloom, Ill-posed Problems for Integro-differential Equations in Mechanics and Electromagnetic Theory, SIAM, Philadelphia, 1981.
  * [18] W. E. Kastenberg and P.L.Chambre, On the’ stability of nonlinear space dependent reactor kinetics, Nucl. Sci. Eng. 31(1968), 67-79.
  * [19] MacCamy, R.C. , A model for one dimensional nonlinear viscoelasticity, Q. Appl. Math. 35 (1977) 21–33.
  * [20] A. Narain and D. Joseph, Classification of linear viscoelastic solids based on a failure criterion, J. Elasticity, 14(1984), 19-26.
  * [21] W. Hrusa, Global existence and asymptotic stability for semilinear hyperbolic Volterra equation with large initial data, SI AM J. Math. Anal. 16(1985), 110-134.
  * [22] J. Pruss, Positivity and regularity of hyperbolic Volterra equations in Banach spaces, Math. Ann. 279(1987), 317-344.
  * [23] V. Capasso, Asymptotic stability for an integrodifferential reaction diffusion system, J. Math. Anal. Appl. 103(1984), 575-588.
  * [24] V. Daftardar-Gejji, H. Jafari, An iterative method for solving non linear functional equations, J. Math. Anal. Appl. 316 (2006) 753–763.
  * [25] A.A. Kilbas, H.M. Srivastava, J.J. Trujillo, Theory and Applications of Fractional Differential Equations, in: North-Holland Mathematics Studies, vol. 204, 2006.
  * [26] V. Daftardar-Gejji, S. Bhalekar, Solving fractional diffusion-wave equations using the New Iterative Method, Frac. Calc. Appl. Anal. 11 (2008) 193–202.
  * [27] S. Bhalekar, V. Daftardar-Gejji, New iterative method: Application to partial differential equations, Appl. Math. Comput. 203 (2008) 778–783.
  * [28] V. Daftardar-Gejji, S. Bhalekar, Solving fractional boundary value problems with Dirichlet boundary conditions, Comput. Math. Appl. 59 (2010) 1801–1809.
  * [29] S. T. Mohyud-Din, A. Yildirim, M. M. Hosseini, An iterative algorithm for fifth-order boundary value problems, World Appl. Sci. J. 8 (2010) 531–535.
  * [30] S. Bhalekar, V. Daftardar-Gejji, Solving evolution equations using a new iterative method, Numer. Methods Partial Differential Equations 26 (2010) 906–916.
  * [31] S. Bhalekar, V. Daftardar-Gejji, Solving a system of nonlinear functional equations using revised new iterative method, Int. J. Computa. Math. Sci. 6 (2012) 127–131.
  * [32] J. Patade, S. Bhalekar, Approximate analytical solutions of Newell-Whitehead-Segel equation using a new iterative method, World J. Modell. Simul. 11 (2) (2015) 94-103.
  * [33] S. Bhalekar, V. Daftardar-Gejji, Solving fractional order logistic equation using a new iterative method, Int. J. Differ. Equ. 2010 (2012) (Art. ID 975829).
  * [34] S. Bhalekar, V. Daftardar-Gejji, Numeric-analytic solutions of dynamical systems using a new iterative method, J. Appl. Nonlin. Dyn. 1 (2012) 141–158.
  * [35] V. Daftardar-Gejji, Y. Sukale, S Bhalekar, A new predictor-corrector method for fractional differential equations, Appl. Math. Comput. 244 (2014) 158–182.
  * [36] V. Daftardar-Gejji, Y. Sukale, S Bhalekar, Solving fractional delay differential equations: A new approach, Fract. Calc. Appl. Anal. 16 (2015) 400–418.
  * [37] J. Patade, S. Bhalekar, A new numerical method based on Daftardar-Gejji and Jafari technique for solving differential Equations, World J. Modell. Simul. 11 (2015) 256–271.
  * [38] S. Bhalekar , V. Daftardar-Gejji, Convergence of the new iterative method, Int. J. Differ. Equ. 2011 (2011).
  * [39] M. K. Jain, Numerical methods for scientific and engineering computation, New Age International, New Delhi 2012.
  * [40] S. Elaydi, An introduction to difference equations, Springer Science & Business Media, 2005.
  * [41] J. T. Edwards, N. J. Ford, J. A. Roberts, Bifurcations in numerical methods for Volterra integro-differential equations, Int. J. bifurcation and chaos. 13(11)(2003) 3255–71.


[affiliation=1]Viet AnhTrinh [affiliation=2]RosySouthwell
[affiliation=1]YiwenGuan [affiliation=1]XinluHe [affiliation=2]ZhiyongWang
[affiliation=1]JacobWhitehill

# Discrete Multimodal Transformers with a Pretrained Large Language Model for
Mixed-Supervision Speech Processing

###### Abstract

Recent work on discrete speech tokenization has paved the way for models that
can seamlessly perform multiple tasks across modalities, e.g., speech
recognition, text to speech, speech to speech translation. Moreover, large
language models (LLMs) pretrained from vast text corpora contain rich
linguistic information that can improve accuracy in a variety of tasks. In
this paper, we present a decoder-only Discrete Multimodal Language Model
(DMLM), which can be flexibly applied to multiple tasks (ASR, T2S, S2TT, etc.)
and modalities (text, speech, vision). We explore several critical aspects of
discrete multimodal models, including the loss function, weight
initialization, mixed training supervision, and codebook. Our results show
that DMLM benefits significantly, across multiple tasks and datasets, from a
combination of supervised and unsupervised training. Moreover, for ASR, it
benefits from initializing DMLM from a pretrained LLM, and from a codebook
derived from Whisper activations.

###### keywords:

discrete speech tokens, multimodal LLM

##  1 Introduction

The recent introduction of discrete speech tokens [1, 2] has facilitated the
integration of large language model (LLM) capabilities into speech models.
Consequently, these speech models not only possess knowledge from the labeled
speech-text pairs on which they were trained, but also inherit the extensive
linguistic knowledge from the LLM. Furthermore, discrete tokenization opens
the possibility for the seamless integration of speech, text, and vision
within a single framework, enabling the model to perform multiple translation
tasks from one modality to another (see Related Work section below). However,
within this new research direction, crucial questions arise regarding how to
manage training when dealing with multi-modal inputs encompassing speech,
text, and/or vision; whether mixed supervision (supervised and/or
unsupervised) on multiple tasks can benefit prediction accuracy; and how the
discrete codebook should be created.

In this paper, we contribute to the growing field of discrete speech models
and (1) propose a Discrete Multimodal Language Model (DMLM) (see Figure 1),
which is a unified model that can flexibly translate data from one modality to
another to perform multiple tasks for speech processing and more. DMLM is one
of the first discrete token-based decoder-only models that can translate to
and from text, speech, and images. As such, it can benefit from pretrained
large language models, such as OPT [3], that were optimized on vast text
corpora. We illustrate DMLM on the tasks of automatic speech recognition
(ASR), text-to-speech (T2S), speech-to-text translation (S2TT) from 21
different languages, and image captioning (I2T). In addition, we (2) propose a
length-normalized tri-modal loss, a modification of the standard cross-entropy
loss for large language models. This adaptation enhances the efficiency of
training multi-modal systems, particularly when diverse modalities exhibit
significant differences in sequence length. We (3) conduct extensive
experiments to explore how mixed-supervision training improves performance on
ASR, T2S, S2TT and I2T. Finally, (4) we introduce a discrete codebook based on
Whisper [4] and compare with current start-of-the-art codebook [2] for the ASR
task.

Figure 1: Discrete Multimodal Language Model (DMLM): The input (text, speech,
image) is tokenized by a modality-specific tokenizer/codec. The discrete token
sequence ${\bf x}_{1},{\bf x}_{2},\ldots$ is concatenated to a task
description and then passed to DMLM for processing. Training is conducted
based on next-token prediction. Inference can be performed for a variety of
tasks (ASR, S2TT, T2S, I2T, etc.).

##  2 Related Work

Discrete Speech Tokens: Motivated by advances in text generation via large
language models as well as the goal of crafting more robust and space-
efficient speech representations, researchers have recently explored ways of
discretizing input speech [1, 5, 2]. One encoder-decoder method [1] employs an
encoder to compress speech into embeddings, and then a decoder to reconstruct
the input speech signal. The embeddings undergo quantization and are expressed
as a combination of basis vectors in a codebook. While this method can yield
high-quality speech with a specific number of quantizers (e.g., 8 or 12
quantizers), it requires a substantial number of tokens (codebook IDs) to
capture just one second of speech – typically in the range of $8\times 65$ or
$8\times 75$ discrete tokens per second. An alternative discretization
strategy [2] involves clustering on a particular layer of a pre-trained speech
model, using the clustering ID as the discrete token index. In particular, [2]
cluster the activations of a speech model trained through self-supervised
learning with contrastive and diversity loss [6, 7].

Speech Processing Based on Pretrained LLMs: Advancements in large language
models (LLMs) have prompted efforts to leverage pretrained LLMs for speech-
related tasks [8, 2]. For instance, [8] utilize PaLM-2 [9], a non-open-source
LLM, by fine-tuning it on speech recognition and speech-to-speech translation
(S2ST) datasets. However, they employ loss masking on the Automatic Speech
Recognition (ASR) task. Loss masking is a special case of the modal-norm loss
function that we use (see Section 3); in particular, it corresponds to setting
one or more of the $\lambda$ values to 0. In our work, instead of hard-coding
these $\lambda$ to a particular value, we optimize them using hyperparameter
search to reach potentially more accurate models. VoxtLM [10] use OPT, which
is an open-source LLM [3], training it on unsupervised speech data with
balanced sampling, ensuring that the number of speech tokens is comparable to
the number of text tokens to address length differences. Our approach, on the
other hand, ensures that the loss is independent of token length for each
modality (speech, text), which may enhance training stability across batches.

Multimodal processing: Efforts in training multimodal models have
predominantly focused on the integration of text and vision. An effective
approach to combining text and vision is to leverage pretrained single modals
for each modality separately. Subsequently, one can either perform cross-
attention between these domains within the embedding space [11] or concatenate
the embeddings of the two modalities [12]. This methodology has been extended
to incorporate audio in [13].

While various approaches exist for combining latent spaces across different
modalities, some efforts have explored the integration of discrete tokens [14,
15, 16]. Specifically, discrete tokens have been employed in tasks related to
vision and text [15, 17, 18]. What sets our work apart is our emphasis on a
tri-modal approach, involving speech, vision, and text. In parallel with our
research, [16] has also introduced a tri-modal framework. However, our
approach distinguishes itself by harnessing the capabilities of LLM and
addressing multilingual aspects, particularly in tasks such as speech-to-
speech translation. Another advantageous aspect of our model is its ability to
facilitate input as a combination of three modalities: speech, vision, and
text. In contrast, in [16], the input is limited to either speech, vision, or
text individually, lacking the capability to accommodate a combination of
these modalities.

##  3 Discrete Multimodal Language Model

Our discrete multimodal language model (DMLM) is a Transformer decoder in the
manner of open-source models such as OPT [3]; as such, it can harness
knowledge from LLM pretraining on large text corpora. The model both inputs
and outputs discrete tokens across multiple modalities. It can also translate
to and from multiple languages. The model can be used for both speech- and
non-speech processing tasks.

Let $\mathcal{T}$, $\mathcal{S}$, and $\mathcal{I}$ represent the (mutually
disjoint) sets of discrete tokens for the text, speech, and image modalities,
respectively, and let $\mathcal{D}=\mathcal{T}\cup\mathcal{S}\cup\mathcal{I}$.
These tokens are computed from raw inputs by modality-specific
codecs/tokenizers, e.g., Seamless [2] for audio, OPT-125 for text, and DALL-E
[15] for images. Associated with each discrete input token ${\bf x}_{i}$ is a
$d$-dimensional embedding vector, which we can implement within the first
layer of the decoder by extending the existing embedding matrix (of shape
$d\times|\mathcal{T}|$) of a pretrained LLM by adding more columns, with a
final shape of $d\times|\mathcal{D}|$. Moreover, the space of possible output
tokens is also $\mathcal{D}$. To convert output tokens to speech, we can use a
vocoder such as [19].

We frame each task (ASR, T2S, S2TT, etc.) as a next-token prediction problem.
The data is organized as a sequence of tokens ${\bf x}_{1},{\bf
x}_{2},\ldots$, and depending on the specific task at hand, the number and
order of modalities may vary. During training the model to translate from a
modality A to modality B, the general format is: Task InputTokens EndA
OutputTokens EndB

where EndA and EndB represents the end delimiter tokens for the source and
target modality, respectively. Note that more advanced uses of the model are
also possible (e.g., input both audio and images to perform context-dependent
speech recognition), but we leave these to future work.

Loss function: Naively applying the cross-entropy loss summed over the
multimodal sequence of tokens ${\bf x}_{1},{\bf x}_{2},\ldots$ leads to poor
results due to large differences in the number of tokens between domains. For
ASR, for example, the number of audio tokens is typically much larger than the
corresponding number of text tokens. To account for these differences, we
devised a new loss function that normalizes the lengths of tokens within each
modality:

$\mathcal{L}=-\frac{\lambda_{\mathcal{S}}}{n_{\mathcal{S}}}\sum_{j\in\mathcal{S}}{p_{j}\log\hat{p}_{j}}-\frac{\lambda_{\mathcal{T}}}{n_{\mathcal{T}}}\sum_{j\in\mathcal{T}}{p_{j}\log\hat{p}_{j}}-\frac{\lambda_{\mathcal{I}}}{n_{\mathcal{I}}}\sum_{j\in\mathcal{I}}{p_{j}\log\hat{p}_{j}}$

where $\hat{p}$ is the vector of softmax probabilities over the entire output
token space ($\mathcal{D}$), and $\lambda_{\mathcal{T}}$,
$\lambda_{\mathcal{S}}$, and $\lambda_{I}$ are the weights of text, speech,
and image modalities, respectively. Meanwhile, $n_{\mathcal{T}}$,
$n_{\mathcal{S}}$, $n_{\mathcal{I}}$ represent the number of text, speech and
image tokens. This loss function resembles that in DALL-E, albeit with a
distinction: In DALL-E [15], the loss is weighted between text and vision
tokens, and the weights $\lambda$ are determined heuristically. In our
approach, we utilize hyperparameter search to determine the contribution of
each modality to the loss.

Example inputs & outputs: Figure 2 shows example inputs & outputs of a DMLM
model, trained on the datasets mentioned in Section 4, when applied to
different tasks. Note that for these examples, different model weights were
optimized for each task; however, in principle a single set of DMLM model
weights could be trained to perform all tasks simultaneously.

Figure 2: Example inputs/outputs of a trained DMLM.

Also, for the specific task of ASR, we performed supervised training
LibriSpeech (960 hours), TED-LIUM (777 hours), SPGI (5000 hours), along with
unsupervised on a 3% subset of Gigaspeech, and thereby produced an ASR model
that attains 10.6% word error rate (WER) on LibriSpeech-other. This is
slightly better than the [10] VOXTLM model, which was trained on 4x as much
supervised and 160x as much unsupervised data and attains a WER of 11.9%.

##  4 Experiments

Before presenting experiments on (1) the loss function weights $\lambda$, (2)
model pretraining, (3) unsupervised training, and (4) tokenization codebooks,
we first describe the experimental setup.

Datasets: For speech recognition (ASR), we used the following datasets: (1)
Librispeech [20], encompassing 960 hours (460 hours clean, 500 hours noisy) of
speech from public-domain audio books; (2) AMI Meeting [21], featuring 100
hours of speech from meetings; (3) TED-LIUM ASR corpus [22] containing a total
of 777 hours of recordings of TED talks; (4) SPGISpeech [23], encompassing
5000 hours of corporate earnings calls along with their associated
transcriptions; and (5) Gigaspeech [24]. For text-to-speech (T2S) tasks, we
used LibriTTS-R [25], containing 585 hours of speech. For the speech-to-text-
translation (S2TT) tasks, we used the CoVoST2 dataset [26]. CoVoST2 consists
of 861 hours of transcribed speech from 21 languages and the corresponding
English translation. Finally, for image captioning (I2T) we used COCO [27].

Network, Hyperparameter Settings, & Codecs: We employ the OPT model [3] as our
pretrained large language model (LLM) with 125 million parameters. We leverage
Hugging Face [28] as our training framework. Our experiments run with a batch
size of 4, a learning rate set at $5\times 10^{-5}$, weight decay of $0.0001$,
and a patience of 1. Model training is conducted on a single Nvidia $40$ GB or
$80$ GB A100 GPU.

To convert from raw waveform to discrete speech tokens, we use the pretrained
Seamless codec [2] except in Experiment 4, where we compare to Whisper [4]. To
convert from images to discrete tokens and vice-versa, we use DALL-E [15].

Evaluation and Metrics: For speech recognition, we utilize the standard WER
metric. For text to speech, we transcribe the synthesized speech and use
pretrained Whisper large-v2 to transcribe and evaluate the character error
rate (CER). For image captioning, we use CIDEr [29], BLEU-4 [30] and METEOR
[31]. For speech-to-text translation, we use BLEU-4 and sacreBLEU [32]

###  4.1 Experiment 1: Tri-Modal Loss Function

Table 1: Effect of $\lambda_{\mathcal{S}},\lambda_{\mathcal{T}}$ on ASR WER(%) of LibriSpeech-clean test set, after training on LibriSpeech (100 hours). Speech weight $\lambda_{\mathcal{S}}$ | Text weight $\lambda_{\mathcal{T}}$ | LS-c
---|---|---
0.25 | 0.93 | 12.9
1.00 | 1.00 | 21.9
0.00 | 1.00 | 15.7

Early during development of DMLM we found that the loss function played a
critical role in ensuring the stability of training. In particular, it was
necessary to select the $\lambda$ to account for the varying lengths of tokens
depending on modality. Hence, in our first experiment, we assess how
$\lambda_{\mathcal{S}},\lambda_{\mathcal{T}}$ affect the test accuracy of DMLM
when applied to speech recognition. In particular, we train DMLM for ASR on
100 hours of LibriSpeech, and we compare (1) heuristically setting
$\lambda_{\mathcal{S}}=\lambda_{\mathcal{T}}$ (similar to Votxtlm [10]); (2)
setting $\lambda_{\mathcal{S}}=0,\lambda_{\mathcal{T}}=1.00$ (as with Google
AudioPaLM [8]); and (3) automatic hyperparameter selection using Ax
(https://ax.dev). With Ax, we conducted 25 trials using early stopping on the
LibriSpeech-clean development set with WER as the metric to choose the best
model.

Results: Ax selected hyperparameters $\lambda_{\mathcal{S}}=0.25$ and
$\lambda_{\mathcal{T}}=0.93$. We thus used these values in this as well as all
other experiments in this paper. Furthermore, as a heuristic choice, we also
set $\lambda_{\mathcal{I}}=0.25$ for image captioning in Experiment 3.

Table 1 shows the Word Error Rate (WER) on the LibriSpeech-clean test set for
different values of $\lambda_{\mathcal{S}},\lambda_{\mathcal{T}}$. It is clear
that applying a loss to the predictions of both speech and text tokens yields
the best WER (12.9%), outperforming both the setting (used in AudioPalm [8])
where only text is modeled (15.7%) and the setting (used in VoxtLM [10]) with
balanced weights between speech and text (21.9%).

###  4.2 Experiment 2: Pretraining the LLM

A potential benefit of a decoder-only model such as DMLM is that it can
harness the knowledge from pre-trained LLMs such as OPT. Here, we compare a
DMLM with random initialization to a model initialized from OPT. In
particular, we train DMLM for speech recognition (ASR) on all 960 hours of
LibriSpeech, and we evaluated it on the LibriSpeech-clean, LibriSpeech-other,
and AMI test sets.

Results: The results in Table 2 illustrate a consistent accuracy improvement
due to LLM pre-training. Notably, the LLM pre-training makes a large
difference for the out-of-domain AMI dataset (WER of 94.8% vs. 658.7%).

Table 2: Effect of model pretraining on ASR WER(%), comparing random initial model weights to OPT-pretrained weights, followed by training/fine-tuning on LibriSpeech (960 hours). Initialization | LS-clean | LS-other | AMI
---|---|---|---
Random | 7.1 | 14.3 | 658.7
With LLM | 5.6 | 12.5 | 94.8

###  4.3 Experiment 3: Mixed-Supervision Training

In contrast to state-of-the-art speech recognition systems such as Whisper,
DMLM can be trained in an unsupervised manner for next-token prediction,
without needing labeled pairs that link across modalities. In this experiment
we explore the utility of mixed-supervision training for speech recognition
(ASR), speech synthesis (T2S), and image captioning (I2T) tasks.

For ASR, we compare training only supervised on LibriSpeech (either 100 hours
or 960 hours), to training both supervised on LibriSpeech as well as
unsupervised with either speech-only data from the AMI training set,
LibriSpeech-500 speech-only data, or Gigaspeech (we used just 3% of the whole
dataset due to computational constraints) text-only data. For T2S, we train
either just supervised on LibriTTS-R, or both supervised and unsupervised
using Gigaspeech (4% of the whole dataset). For I2T, we compare training just
supervised on COCO, to supervised on COCO combined with LibriSpeech 960 (text-
only) data.

Results:

Table 3: Effect of unsupervised training (plus supervised training on LibriSpeech-100) on ASR WER(%) on the LibriSpeech-other test set. Training set | Test set
---|---
Supervised | Unsupervised | LS-o
LS-100 | None | 28.6
LS-100 | 20% of LS-500 (Speech only) | 25.1
Table 4:  Effect of unsupervised training (plus supervised training on LibriSpeech-960) on ASR WER(%) on the AMI test set. Training set | Test set
---|---
Supervised | Unsupervised | AMI
LS-960 | None | 94.8
LS-960 | AMI (Speech only) | 92.6
LS-960 | 3% Gigaspeech (Text only) | 59.3

Table 3 shows that unsupervised training with speech-only data from
LibriSpeech-500, in conjunction with supervised training on LibriSpeech-100,
improves WER on LibriSpeech-other for ASR. Similarly, Table 4 shows that
leveraging speech-only data from AMI, in conjunction with supervised training
on LibriSpeech-960, improves WER on the AMI test set. Notably, unsupervised
training on a text-only dataset (3% of Gigaspeech) decreases the WER even
further to 59.3%.

Table 5: Effect of unsupervised training (plus supervised training on LibriTTS-R) on T2S CER(%) on the LibriTTS-R test set Training set | Test set
---|---
Supervised | Unsupervised | LibriTTS-R
LibriTTS-R | None | 22.9
LibriTTS-R | 4% Gigaspeech (Speech only) | 20.2
Table 6: Effect of mixed-supervision training on S2TT B@4 (BLEU-4) and sacreBLEU. BLEU and sacreBLEU scores are calculated based on the average performance across 21 languages within the CoVoST2 X$\rightarrow$en test set. Training set | Test set
---|---
Supervised | Unsupervised | B@4($\uparrow$) | sacreBLEU ($\uparrow$)
CoVoST2 | None | 5.71 | 5.75
CoVoST2 | LS-960 (speech) | 5.74 | 5.76
| CoVoST2 (S2TT)
---
\+ LS-960 (ASR)
None | 6.18 | 6.22
Table 7: Effect of unsupervised training (plus supervised training on COCO) on I2T metrics on the COCO test set. B@4: BLEU-4, M: METEOR. Training set | Test set
---|---
Supervised | Unsupervised |  CIDEr ($\uparrow$) |  B@4 ($\uparrow$) |  M ($\uparrow$)
COCO | None | 58.1 | 19.5 | 19.2
COCO | LS-960 (text) | 61.5 | 20.6 | 19.9
COCO | SBU (image) | 62.8 | 21.1 | 20.1

T2S: Table 5 shows that, by adding additional 4% of Gigaspeech for speech-only
unsupervised training, the CER is reduced from 22.9% to 20.2%.

S2TT: Table 6 shows a small improvement of speech-to-text translation with
unsupervised speech-only data from LibriSpeech-960 (line #2), as well as a
larger improvement (line #3) by mixing two forms of supervised training
(CoVoST2 for S2TT and LibriSpeech-960 for ASR).

I2C: Table 7 shows that, by adding the text-only data from LS-960, we improve
the CIDEr score from 58.1 to 61.5. Also, including image-only data from SBU
further elevates the CIDEr score to 62.8.

###  4.4 Experiment 4: Codebook

Finally, we investigate the use of an alternative codebook of speech tokens
derived from the state-of-the-art Whisper [4] model. We use the same approach
as [2], namely K-means clustering on the activations from a single layer of
the encoder – we used the final layer of Whisper large-v2. We use the mini-
batch k-means algorithm in scikit-learn [33] with k=10000 and 64 utterances
per minibatch. Tokens are then generated by passing audio to the encoder and
assigning a token based on the nearest cluster for each timestep. For
comparison, we train equivalent k-means models on XLS-R as used in [2]. We use
Librispeech-960h to train the k-means for both Whisper and XLS-R. Results:

Table 8: WER(%) using Whisper-derived tokens  Model | K-means data | Training set | LS-c | LS-o
---|---|---|---|---
Seamless | LS-960 | LS-960 | 31.8 | 61.2
Whisper | LS-960 | LS-960 | 14.9 | 23.1

We find that using Whisper features as the basis for discrete tokens yields a
lower WER on Librispeech test than when the units are trained on XLS-R (Table
8). This improvement may be attributed to the fact that Whisper was trained
with ASR as one of its tasks, thereby learning a more effective representation
for this task.

##  5 Discussion and Conclusions

We introduced DMLM, a decoder-only framework that can harness pre-trained LLMs
and can be used for multiple tasks across multiple modalities. DMLM uses
discrete speech, vision and text tokens as the inputs and treats all the tasks
as a next-token prediction task, thus enabling the framework to train with
mixed supervision in a unified manner. We demonstrate the effectiveness of
using unsupervised data on various of tasks: speech recognition, speech
synthesis, speech to text translation and image captioning. We further
investigate the choice of embedding space leading to discrete speech tokens,
in particular between Whisper and XLS-R for speech recognition task.

Our work leads to interesting future research directions, e.g., audiovisual
ASR through multimodal input streams.

## References

  * [1] N. Zeghidour, A. Luebs, A. Omran, J. Skoglund, and M. Tagliasacchi, ``Soundstream: An end-to-end neural audio codec,'' _IEEE/ACM Transactions on Audio, Speech, and Language Processing_ , vol. 30, pp. 495–507, 2021.
  * [2] L. Barrault, Y.-A. Chung, M. C. Meglioli, D. Dale, N. Dong, P.-A. Duquenne, H. Elsahar, H. Gong, K. Heffernan, J. Hoffman _et al._ , ``Seamlessm4t-massively multilingual & multimodal machine translation,'' _arXiv preprint arXiv:2308.11596_ , 2023.
  * [3] S. Zhang, S. Roller, N. Goyal, M. Artetxe, M. Chen, S. Chen, C. Dewan, M. Diab, X. Li, X. V. Lin _et al._ , ``Opt: Open pre-trained transformer language models,'' _arXiv preprint arXiv:2205.01068_ , 2022.
  * [4] A. Radford, J. W. Kim, T. Xu, G. Brockman, C. McLeavey, and I. Sutskever, ``Robust speech recognition via large-scale weak supervision,'' in _International Conference on Machine Learning_ , 2023.
  * [5] Z. Borsos, R. Marinier, D. Vincent, E. Kharitonov, O. Pietquin, M. Sharifi, D. Roblek, O. Teboul, D. Grangier, M. Tagliasacchi _et al._ , ``Audiolm: a language modeling approach to audio generation,'' _IEEE/ACM Transactions on Audio, Speech, and Language Processing_ , 2023.
  * [6] A. Babu, C. Wang, A. Tjandra, K. Lakhotia, Q. Xu, N. Goyal, K. Singh, P. von Platen, Y. Saraf, J. M. Pino, A. Baevski, A. Conneau, and M. Auli, ``Xls-r: Self-supervised cross-lingual speech representation learning at scale,'' in _Interspeech_ , 2021.
  * [7] A. Baevski, Y. Zhou, A. Mohamed, and M. Auli, ``wav2vec 2.0: A framework for self-supervised learning of speech representations,'' _Advances in neural information processing systems_ , vol. 33, pp. 12 449–12 460, 2020.
  * [8] P. K. Rubenstein, C. Asawaroengchai, D. D. Nguyen, A. Bapna, Z. Borsos, F. d. C. Quitry, P. Chen, D. E. Badawy, W. Han, E. Kharitonov _et al._ , ``Audiopalm: A large language model that can speak and listen,'' _arXiv preprint arXiv:2306.12925_ , 2023.
  * [9] R. Anil, A. M. Dai, O. Firat, M. Johnson, D. Lepikhin, A. Passos, S. Shakeri, E. Taropa, P. Bailey, Z. Chen _et al._ , ``Palm 2 technical report,'' _arXiv preprint arXiv:2305.10403_ , 2023.
  * [10] S. Maiti, Y. Peng, S. Choi, J.-w. Jung, X. Chang, and S. Watanabe, ``Voxtlm: unified decoder-only models for consolidating speech recognition/synthesis and speech/text continuation tasks,'' _arXiv preprint arXiv:2309.07937_ , 2023.
  * [11] J.-B. Alayrac, J. Donahue, P. Luc, A. Miech, I. Barr, Y. Hasson, K. Lenc, A. Mensch, K. Millican, M. Reynolds _et al._ , ``Flamingo: a visual language model for few-shot learning,'' _Advances in Neural Information Processing Systems_ , vol. 35, pp. 23 716–23 736, 2022.
  * [12] H. Liu, C. Li, Q. Wu, and Y. J. Lee, ``Visual instruction tuning,'' _Advances in neural information processing systems_ , vol. 36, 2024.
  * [13] S. Wu, H. Fei, L. Qu, W. Ji, and T.-S. Chua, ``Next-gpt: Any-to-any multimodal llm,'' _arXiv preprint arXiv:2309.05519_ , 2023.
  * [14] G. Team, R. Anil, S. Borgeaud, Y. Wu, J.-B. Alayrac, J. Yu, R. Soricut, J. Schalkwyk, A. M. Dai, A. Hauth _et al._ , ``Gemini: a family of highly capable multimodal models,'' _arXiv preprint arXiv:2312.11805_ , 2023.
  * [15] A. Ramesh, M. Pavlov, G. Goh, S. Gray, C. Voss, A. Radford, M. Chen, and I. Sutskever, ``Zero-shot text-to-image generation,'' in _International Conference on Machine Learning_. PMLR, 2021, pp. 8821–8831.
  * [16] M. Kim, J.-w. Jung, H. Rha, S. Maiti, S. Arora, X. Chang, S. Watanabe, and Y. M. Ro, ``Tmt: Tri-modal translation between speech, image, and text by processing different modalities as different languages,'' _arXiv preprint arXiv:2402.16021_ , 2024.
  * [17] Y. Ge, S. Zhao, Z. Zeng, Y. Ge, C. Li, X. Wang, and Y. Shan, ``Making llama see and draw with seed tokenizer,'' _arXiv preprint arXiv:2310.01218_ , 2023.
  * [18] T. Du, Y. Wang, and Y. Wang, ``On the role of discrete tokenization in visual representation learning,'' in _The Twelfth International Conference on Learning Representations_ , 2023.
  * [19] J. Kong, J. Kim, and J. Bae, ``Hifi-gan: Generative adversarial networks for efficient and high fidelity speech synthesis,'' _Advances in Neural Information Processing Systems_ , vol. 33, pp. 17 022–17 033, 2020.
  * [20] V. Panayotov, G. Chen, D. Povey, and S. Khudanpur, ``Librispeech: an asr corpus based on public domain audio books,'' in _2015 IEEE international conference on acoustics, speech and signal processing (ICASSP)_. IEEE, 2015, pp. 5206–5210.
  * [21] J. Carletta, S. Ashby, S. Bourban, M. Flynn, M. Guillemot, T. Hain, J. Kadlec, V. Karaiskos, W. Kraaij, M. Kronenthal _et al._ , ``The ami meeting corpus: A pre-announcement,'' in _International workshop on machine learning for multimodal interaction_. Springer, 2005, pp. 28–39.
  * [22] F. Hernandez, V. Nguyen, S. Ghannay, N. Tomashenko, and Y. Esteve, ``Ted-lium 3: Twice as much data and corpus repartition for experiments on speaker adaptation,'' in _Speech and Computer: 20th International Conference, SPECOM 2018, Leipzig, Germany, September 18–22, 2018, Proceedings 20_. Springer, 2018, pp. 198–208.
  * [23] P. K. O'Neill, V. Lavrukhin, S. Majumdar, V. Noroozi, Y. Zhang, O. Kuchaiev, J. Balam, Y. Dovzhenko, K. Freyberg, M. D. Shulman _et al._ , ``Spgispeech: 5,000 hours of transcribed financial audio for fully formatted end-to-end speech recognition,'' _arXiv preprint arXiv:2104.02014_ , 2021.
  * [24] G. Chen, S. Chai, G. Wang, J. Du, W.-Q. Zhang, C. Weng, D. Su, D. Povey, J. Trmal, J. Zhang _et al._ , ``Gigaspeech: An evolving, multi-domain asr corpus with 10,000 hours of transcribed audio,'' _arXiv preprint arXiv:2106.06909_ , 2021.
  * [25] Y. Koizumi, H. Zen, S. Karita, Y. Ding, K. Yatabe, N. Morioka, M. A. U. Bacchiani, Y. Zhang, W. Han, and A. Bapna, ``Libritts-r: Restoration of a large-scale multi-speaker tts corpus,'' 2023.
  * [26] C. Wang, A. Wu, and J. Pino, ``Covost 2 and massively multilingual speech-to-text translation,'' _arXiv preprint arXiv:2007.10310_ , 2020.
  * [27] T.-Y. Lin, M. Maire, S. Belongie, J. Hays, P. Perona, D. Ramanan, P. Dollár, and C. L. Zitnick, ``Microsoft coco: Common objects in context,'' in _Computer Vision–ECCV 2014: 13th European Conference, Zurich, Switzerland, September 6-12, 2014, Proceedings, Part V 13_. Springer, 2014, pp. 740–755.
  * [28] T. Wolf, L. Debut, V. Sanh, J. Chaumond, C. Delangue, A. Moi, P. Cistac, T. Rault, R. Louf, M. Funtowicz _et al._ , ``Huggingface's transformers: State-of-the-art natural language processing,'' _arXiv preprint arXiv:1910.03771_ , 2019.
  * [29] R. Vedantam, C. Lawrence Zitnick, and D. Parikh, ``Cider: Consensus-based image description evaluation,'' in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2015, pp. 4566–4575.
  * [30] K. Papineni, S. Roukos, T. Ward, and W.-J. Zhu, ``Bleu: a method for automatic evaluation of machine translation,'' in _Proceedings of the 40th annual meeting of the Association for Computational Linguistics_ , 2002, pp. 311–318.
  * [31] S. Banerjee and A. Lavie, ``Meteor: An automatic metric for mt evaluation with improved correlation with human judgments,'' in _Proceedings of the acl workshop on intrinsic and extrinsic evaluation measures for machine translation and/or summarization_ , 2005, pp. 65–72.
  * [32] M. Post, ``A call for clarity in reporting BLEU scores,'' in _Proceedings of the Third Conference on Machine Translation: Research Papers_. Association for Computational Linguistics, Oct. 2018, pp. 186–191. [Online]. Available: https://aclanthology.org/W18-6319
  * [33] D. Sculley, ``Web-scale k-means clustering,'' in _Proceedings of the 19th International Conference on World Wide Web_ , ser. WWW '10. New York, NY, USA: Association for Computing Machinery, 2010, p. 1177–1178. [Online]. Available: https://doi.org/10.1145/1772690.1772862