\def\firstpage{11}
\def\lastpage{30}%number of pages
\def\mmacorresp{olga.stikoniene@mif.vu.lt}
\def\mmamid{23.18029}
\def\mmayear{24}
\def\firstfootnote{0}
% mma style needs information bellow ======================================
\def\mmaauthor{O. \v{S}tikonien\.{e} \and M. Sapagovas}%Example:% F. Author\sep S. Author \and Th. N. Author
\def\mmatitle{Alternating direction implicit method  for \\ Poisson equation with integral conditions}%Example:% Title of Article
\def\mmarunauthor{O. \v{S}tikonien\.{e} \and M. Sapagovas}%Example:% F. Author\sep S. Author \and Th. N. Author
\def\mmaruntitle{ADI Method for Poisson Equation with Integral Conditions}%Example:% Runtitle of Article
\def\mmaabstract{In this paper, we investigate the convergence of the Peaceman-Rachford Alternating Direction Implicit  method for the system of difference equations,
 approximating the two-dimensional elliptic equations in rectangular domain with nonlocal integral conditions.  The main goal of
the paper is the analysis of spectrum structure of difference eigenvalue problem with nonlocal conditions.
The convergence of iterative method is proved in the case when the system of eigenvectors is complete.
The main results are generalized for the system of difference equations, approximating the differential problem with truncation error $\mathcal{O}(h^4)$.
}
\def\mmakeywords{elliptic equation; integral boundary conditions; finite-difference method; iterative method; eigenvalue problem}%Example:% discrete problem, function
\def\mmaams{35J25; 65N06; 65N22; 35P15}%Example:%34B24; 34B10
\def\mmareceived{November 14, 2022}
\def\mmarevised{Month day, year}
\def\mmaacepted{Month day, year}
\def\mmapublished{Month day, year}

\begin{filecontents*}{mma.def}
\endinput
\end{filecontents*}
% ************* Author's  references  ******************%
\begin{filecontents*}{x.bib}
@preamble{ "\newcommand{\nosort}[1]{}" }

@BOOK{29,
  author =       {K.E. Atkinson},
  title =        {An Introduction to Numerical Analysis},
  publisher =    {Wiley},
  year =         {1988},
  }

@ARTICLE{1964_Beals,
   author = "R.W. Beals",
   title =  "Nonlocal elliptic boundary value problems",
   journal = "Bulletin of the American Mathematical Society",
   year = 1964,
   volume = 70,
   number = "",
   pages = "693--696",
   doi = {10.1090/S0002-9904-1964-11166-6}
}

@ARTICLE{1969_Bitsadze_Samarskii,
   author = "A.V. Bitsadze and A.A. Samarskii" ,
   title = "Some elementary generalizations of linear elliptic boundary value problems",
   journal = "Soviet Math. Dokl.",
   year = 1969,
   volume =10,
   pages = "398--400",
   note ={(in Russian)},
  }

@BOOK{1981_Gordeziani,
  author =       {N. Gordeziani},
  title =        {Solution Methods for a Class of Nonlocal Boundary Value Problems},
  publisher =    {Tbilisi},
  year =         {1981},
  note ={(in Russian)},
  }

@ARTICLE{1987_Sapagovas_Ciegis_1,
   author = "M.P. Sapagovas and R.Yu. \v{C}iegis",
   title = " The numerical solution of some nonlocal problems",
   journal = "Litovsk. Matem. sborn.",
   year = "1987",
   volume = 27,
   number = 2,
   pages = "348--356",
   note ={(in Russian)},
}

@ARTICLE{2002_Sapagovas,
   author = "M.P. Sapagovas" ,
   title = "The eigenvalues of some problem with a nonlocal condition",
   journal = "Differential Equations",
   year = 2002,
   volume = 38,
   number = 7,
   pages = "1020--1026",
   doi={10.1023/A:1021115915575},
   }

@ARTICLE{1990_Il'in_Moiseev,
   author = "V.A. Il'in and E.I. Moiseev",
   title = "Two-dimensional nonlocal boundary value problems for {P}oisson's operator in differential and difference variants",
   journal = "Mat. Mod.",
   year = 1990,
   volume = 2,
   number = 8,
   pages = "139--156",
    note ={(in Russian)},
 }

@ARTICLE{2001_Berikelashvili,
   author = "G. Berikelashvili",
   title =  "Finite difference schemes for some mixed boundary problems",
   journal = "Proceedings of A.Razmadze Mathematical Institute",
   year = 2001,
   volume = 127,
   number = "",
   pages = "77--87",
  }

@ARTICLE{2002_Wang,
   author = "Y. Wang",
   title = "Solutions to nonlinear elliptic equations with a nonlocal boundary condition",
   journal = "Electron. J. Differential Equations",
   year = 2002,
   volume = 2002,
   number = 05,
   pages = "1--16",
   url={http://ejde.math.txstate.edu/Volumes/2002/05/wang.pdf}
}

@ARTICLE{2013_Volkov_Dosiyev_Buranay,
   author = "E.A. Volkov and A.A. Dosiyev and S.C. Buranay",
   title = "On the solution of a nonlocal problem",
   journal = "Computers \& Mathematics with Applications",
volume = {66},
number = {3},
   year = 2013,
   pages = "330--338",
doi = {10.1016/j.camwa.2013.05.010},
  }



@ARTICLE{2014_Ashyralyev_Ozturk,
   author = "A. Ashyralyev and E. Ozturk",
   title = "On a difference scheme of second order of accuracy for the {B}itsadze-{S}amarskii type nonlocal boundary-value problem",
   journal = "Boundary Value Problems",
   year = 2014,
   number = 14,
   volume  = 2014,
   pages = "1--19",
   doi = "10.1186/1687-2770-2014-14"
}

@ARTICLE{2018_Zhou_Yu,
   author = "L. Zhou and H. Yu",
   title = "Error estimate of a high accuracy difference scheme for {P}oisson equation with two integral boundary conditions",
   journal = "Advances in Difference Equations",
   year = 2018,
   volume = "2018",
   number = 225,
   pages = "1--11",
   doi = "10.1186/s13662-018-1682-z"
}

@ARTICLE{2009LMJ49n3o,
   author = "M. Sapagovas and O. \v{S}tikonien\.e",
  title = "A fourth-order alternating-direction method for difference schemes with nonlocal condition",
  journal = {Lithuanian Mathematical Journal},
  year = "2009",
  volume = 49,
  number = 3,
  pages = "309--317",
   doi = {10.1007/s10986-009-9057-5},
  }
@ARTICLE{2011DE47n8,
   author = "M. Sapagovas and A. \v{S}tikonas and  O. \v{S}tikonien\.e",
   title = "Alternating Direction Method for the {P}oisson Equation with Variable Weight Coefficients in an Integral Condition",
  journal = {Differential Equations},
   year = 2011,
   volume = 47,
   number = 8,
   pages = "1176--1187",
   doi = {10.1134/S0012266111080118},
}

@ARTICLE{2014NA19n3o,
   author = "O. {\v{S}}tikonien\.{e} and M. Sapagovas and R. \v{C}iupaila",
   title =  " On iterative methods for some elliptic equations with nonlocal conditions",
   journal = "Nonlinear Analysis: Modelling and Control",
   year = 2014,
   volume = 19,
   number  = 3,
   pages = "517--535",
   doi = {10.15388/NA.2014.3.14},
}


@ARTICLE{2008_Sapagovas_1,
   author = "M.P. Sapagovas",
   title = {Difference method of increased order of accuracy for the {P}oisson equation with nonlocal conditions},
   journal = {Differential Equations},
   year = 2008,
   volume = 44,
   number = 7,
   pages = "1018--1028",
   doi = {10.1134/S0012266108070148}
 }

@ARTICLE{2013_Ashyralyev_Ozturk,
   author = "A. Ashyralyev and E. Ozturk",
   title = "On a difference scheme of fourth-order of accuracy for the {B}itsadze-{S}amarskii type nonlocal boundary-value problem",
   journal = "Mathematical Methods in the Applied Sciences",
   year = 2013,
   volume = {36},
   number = {8},
   pages = "935--955",
   doi = {10.1002/mma.2650},
  }


@ARTICLE{2019_Dosiyev,
   author = "A.A. Dosiyev",
   title = "Difference method of fourth order accuracy for the {L}aplace equation with multilevel nonlocal conditions",
   journal = "Journal of Computational and Applied Mathematics",
   year = 2019,
  volume = 354,
  number = "",
   pages = "587--596",
   doi = {10.1016/j.cam.2018.04.046},
  }

@ARTICLE{2019_Dosiyev_Reis,
   author = "A.A. Dosiyev and R. Reis",
   title = "A fourth-order accurate difference {D}irichlet problem for the approximate solution of {L}aplace’s equation with integral boundary condition",
   journal = "Advances in Difference Equations",
   year = 2019,
   volume = 2019,
   number = 340,
   pages = "587--596",
   doi = {10.1186/s13662-019-2282-2},
  }

@ARTICLE{2017_Sapagovas_Griskoniene_Stikoniene,
   author = "M. Sapagovas and V. Gri\v{s}konien\.e and O. \v{S}tikonien\.e",
   title = "Application of {M}-matrices theory to numerical investigation of a nonlinear elliptic equation with an integral condition",
   journal = "Nonlinear Analysis: Modelling and Control",
   year = 2017,
   volume=22,
   number = 4,
   pages = "489--504",
   doi = {10.15388/NA.2017.4.5},
  }

  @ARTICLE{2017_Sapagovas_Stikoniene_KJ_RC,
   author = "M. Sapagovas  and O. \v{S}tikonien\.e and K. Jakub\.elien\.e and R. \v{C}iupaila",
   title = "Finite difference method for boundary value problem for nonlinear elliptic equation with nonlocal conditions",
   journal = "Boundary Value Problems",
   year = 2019,
   volume = "2019",
   number = 94,
   pages = "1--16",
   doi = {10.1186/s13661-019-1202-4},
  }

@ARTICLE{2021_Ciupaila_Pupalaige_Sapagovas,
   author = "R. {\v{C}}iupaila and K. Pupalaig\.e and M. Sapagovas",
   title = "On the numerical solution for nonlinear elliptic equations with variable weight coefficients in an integral boundary conditions",
   journal = "Nonlinear Analysis: Modelling and Control",
   volume=26,
   year = "2021",
   number = 4,
   pages = "738--758",
   doi = {10.15388/namc.2021.26.23929},
  }

  @article{2011_Avalishvili_Gordeziani,
  author={G. Avalishvili and M. Avalishvili and D. Gordeziani},
  title="On a nonlocal problem with integral boundary conditions for a multidimensional elliptic equation",
  journal="Applied Mathematics Letters",
  volume=24,
  number=4,
  pages="566--571",
  year=2011,
  doi = {10.1016/j.aml.2010.11.014},
}


  @article{2011_Avalishvili_Gordeziani_1,
  author={G. Avalishvili and M. Avalishvili and D. Gordeziani},
  title="On Integral Nonlocal Boundary Value Problems for some Partial Differential Equations",
  journal="Bulletin of the Georgian National Academy of Sciences",
  volume=5,
  number=1,
  pages="31--37",
  year=2011,
  doi = {10.1016/j.aml.2010.11.014},
}

@ARTICLE{2013_Berikelashvili_Khomeriki,
   author = "G. Berikelashvili and N. Khomeriki",
  title = "On a numerical solution of one nonlocal boundary-value problem with mixed {D}irichlet-{N}eumann conditions",
  journal = {Lithuanian Mathematical Journal},
  year = "2013",
  volume = 53,
  number = 4,
  pages = "367--380",
  doi = {10.1007/s10986-013-9214-8},
}


@ARTICLE{2014_Berikelashvili_Khomeriki,
   author = "G. Berikelashvili and  N. Khomeriki",
   title =  "On the convergence rate of a difference solution of the {P}oisson equation with fully nonlocal constraints",
   journal = "Nonlinear Analysis: Modelling and Control",
   year = 2014,
   volume = 19,
   number  = 3,
   pages = "367--381",
   doi = {10.15388/NA.2014.3.4},
   }
   @ARTICLE{2005MMA10n4,
   author = "S. Pe\v{c}iulyt\.e and O. \v{S}tikonien\.e and A. \v{S}tikonas",
   title =  "{S}turm-{L}iouville Problem for Stationary Di{ff}erential Operator with Nonlocal Integral Boundary Condition",
   journal = "Math. Model. Anal.",
   year = 2005,
   volume = 10,
   number  = 4,
   pages = "377--392",
     doi={10.1080/13926292.2005.9637295},
 }

 @ARTICLE{2009_Jachimaviciene_et_al,
  author = "J. Jachimavi\v{c}ien\.e and  \v{Z}. Jesevi\v{c}i\={u}t\.e  and M. Sapagovas",
  title = "The stability of finite-difference schemes for a pseudoparabolic equation with nonlocal conditions",
  journal = {Numerical Functional Analysis and Optimization},
  year = "2009",
  volume = 30,
  number = "9",
  pages = "988--1001",
  doi= {10.1080/01630560903405412},
}



@ARTICLE{2019NA24n4,
   author = " K. Bingel\.e and A. Bankauskien\.e and A. \v{S}tikonas" ,
   title = "Spectrum curves for a discrete {S}turm–{L}iouville problem with one integral boundary condition",
   journal = "Nonlinear Analysis: Modelling and Control",
   year = 2019,
   volume = 24,
   number = 5,
   pages = "755--774",
   doi ={10.15388/NA.2019.5.5},
}

@ARTICLE{2008NA13n4,
   author = "S. Pe\v{c}iulyt\.e and O. \v{S}tikonien\.e and A. \v{S}tikonas",
   title =  "Investigation of negative critical points of the characteristic function for problems with nonlocal boundary conditions",
   journal = "Nonlinear Analysis: Modelling and Control",
   year = 2008,
   volume = 13,
   number  = 4,
   pages = "467--490",
   doi ={10.15388/NA.2008.13.4.14552},
}
@ARTICLE{2007MMA12n2,
   author = "S. Pe\v{c}iulyt\.e and A. \v{S}tikonas",
   title =  "On Positive Eigenfunctions of {S}turm-{L}iouville Problem with Nonlocal Two-Point Boundary Condition",
   journal = "Mathematical Modelling and Analysis",
   year = 2007,
   volume = 12,
   number  = 2,
   pages = "215--226",
   doi   = {10.3846/1392-6292.2007.12.215-226},
}

@ARTICLE{2022LMJ62n4,
   author = "A. {\v{S}}tikonas and E. \c{S}en",
   title =  "Asymptotic analysis of {S}turm--{L}iouville problem with {N}eumann and nonlocal two-point boundary conditions",
   journal = "Lithuanian Mathematical Journal",
   year = 2022,
   volume = 62,
   number = 4,
   pages = "519--541",
   doi ={10.1007/s10986-022-09577-6},
}

@BOOK{45, 
author="A.A. Samarskii", 
title="The Theory of Difference Schemes",
publisher="Moscow, Nauka", 
year="1989",
note = {(in Russian)},
 }


\end{filecontents*}
% ************* end Author's references ******************%
\documentclass{mma}%Copyediting,Article in press
%\input{mma.copyed.tex}
%\input{mma.proof.tex}
%\input{mma.press.tex}
%\input{mma.def.tex}
\MMAJCLS
 \makeatletter
 \spnewtheorem*{west}{Test the West}{\bfseries}{\itshape}


%  \usepackage{showkeys}



% \mathtoolsset{showonlyrefs}

%\usepackage{cleveref}
%\usepackage{refcheck}

 \@ifundefined{BibTeX}{%
 \def\BibTeX{\mbox{{\sc Bib}\TeX}}}{}
 

 
 \makeatother
 


% ************* Author's  definitions  ******************%************ Author's  definitions  ******************%
%\usepackage{mandi}
%%\usepackage[OMLmathrm,OMLmathbf]{isomath} % This gives me correct \mathbf{\pi}, but italic \mathrm
%\usepackage{isomath}
%\usepackage{refcheck}

\makeatletter
\newif\if@borderstar
\def\bordermatrix{\@ifnextchar*{%
\@borderstartrue\@bordermatrix@i}{\@borderstarfalse\@bordermatrix@i*}%
}
\def\@bordermatrix@i*{\@ifnextchar[{\@bordermatrix@ii}{\@bordermatrix@ii[()]}}
\def\@bordermatrix@ii[#1]#2{%
\begingroup
\m@th\@tempdima8.75\p@\setbox\z@\vbox{%
\def\cr{\crcr\noalign{\kern 2\p@\global\let\cr\endline }}%
\ialign {$##$\hfil\kern 2\p@\kern\@tempdima & \thinspace %
\hfil $##$\hfil && \quad\hfil $##$\hfil\crcr\omit\strut %
\hfil\crcr\noalign{\kern -\baselineskip}#2\crcr\omit %
\strut\cr}}%
\setbox\tw@\vbox{\unvcopy\z@\global\setbox\@ne\lastbox}%
\setbox\tw@\hbox{\unhbox\@ne\unskip\global\setbox\@ne\lastbox}%
\setbox\tw@\hbox{%
$\kern\wd\@ne\kern -\@tempdima\left\@firstoftwo#1%
\if@borderstar\kern2pt\else\kern -\wd\@ne\fi%
\global\setbox\@ne\vbox{\box\@ne\if@borderstar\else\kern 2\p@\fi}%
\vcenter{\if@borderstar\else\kern -\ht\@ne\fi%
\unvbox\z@\kern-\if@borderstar2\fi\baselineskip}%
\if@borderstar\kern-2\@tempdima\kern2\p@\else\,\fi\right\@secondoftwo#1 $%
}\null \;\vbox{\kern\ht\@ne\box\tw@}%
\endgroup
}


\makeatother
% ************* end Author's definitions ******************%
\begin{document}
\bibliographystyle{plainmma}
\begin{topmatter}
%
%\thankstext{t1}{Title thanks}
%add runauthor at beginning of Your file
%add runtitle at beginning of Your file
%add title at beginning of Your file
%add add all authors at beginning of Your file
\title{\mmatitletext}%Example with thanks:% \mmatitletext\thanks{...}\thanks{article thanks}
\author{Olga \v{S}tikonien\.{e}${}^{a}\COR{\mmacorresp}\ID{0000-0002-0302-3449}$ \and Mifodijus Sapagovas${}^{b} \ID{0000-0002-7139-3468}$}%Example:% First Author${}^{a}$\sep Second Author${}^{a,b}$ \and Third New Author 
%\thanks{author thanks}
\vskip 0.2cm
{\setstretch{0.75}
\institution{${}^{a}$Institute of Applied Mathematics, Vilnius University, Lithuania}
\institution{${}^{b}$Institute of Data Science and Digital Technologies, Vilnius University, Lithuania}
}
\Keywords%\Keywords and AMS add keywords at beginning of Your file
\Email{\mmacorresp}
%\email{${}^{a}$\mmacorresp;${}^{b}$mifodijus.sapagovas@mif.vu.lt}
\end{topmatter}
\MMAJarticle



%
% ************* Text entry area ******************%

\section{Introduction}\label{s:1}


 Boundary value problems for differential equations with various types of nonlocal conditions are currently being studied quite intensively in the theory of differential equations and numerical analysis.

The study of numerical methods for elliptic equations with nonlocal conditions is strongly influenced by two causes.
Firstly, over the past few decades, new mathematical models with nonlocal conditions have been developed for applications in physics, thermoelasticity, ecology, biotechnology,  etc. %biochemistry, engineering,
Secondly, investigating the problems of pure mathematics, several scientific articles have been published on the generalization of classical boundary conditions for elliptic equations \cite{1964_Beals, 1969_Bitsadze_Samarskii}.


The first results on the solution of a two-dimensional elliptic equation with a nonlocal condition were obtained  in \cite{1981_Gordeziani, 1990_Il'in_Moiseev, 1987_Sapagovas_Ciegis_1, 2002_Sapagovas}. This condition was later named the Bitsadze--Samarskii nonlocal condition.
These papers began to investigate  iterative methods for system of difference equations with nonlocal conditions. We note one characteristic feature of such systems.
Due to the nonlocal condition, the matrix of the system of difference equations is not symmetric. But quite often it has some nice properties, for example, all the eigenvalues of the matrix are positive.

   

Many articles  are devoted to the estimation of the error of the finite difference method and the convergence of elliptic equations with various types of
nonlocal conditions  \cite{1990_Il'in_Moiseev, 2001_Berikelashvili, 2002_Wang, 2013_Volkov_Dosiyev_Buranay, 2014_Ashyralyev_Ozturk, 2018_Zhou_Yu}. The
alternating  direction method for system of difference equations with nonlocal conditions is examined in the papers  \cite{2009LMJ49n3o, 2011DE47n8,
2014NA19n3o}.
 In many cases, the matrix of system of difference equations has properties which are typical to M-matrices. Therefore, the theory of M-matrices can be
 applied to the study and solution of problems with nonlocal conditions   \cite{2014NA19n3o,  2017_Sapagovas_Griskoniene_Stikoniene,
 2017_Sapagovas_Stikoniene_KJ_RC, 2021_Ciupaila_Pupalaige_Sapagovas}.
The works     \cite{2008_Sapagovas_1, 2013_Ashyralyev_Ozturk, 2019_Dosiyev, 2019_Dosiyev_Reis}   are   devoted to high-precision finite difference methods for the simplest elliptic equations with nonlocal conditions.

 We consider the nonlocal boundary value problem for two-dimensional Poisson equation in rectangular domain
\begin{equation}\label{eq:1.1}
 \dfrac{\partial^2 u}{\partial x^2}+ \dfrac{\partial^2 u}{\partial y^2}=-f(x, y), \quad (x, y)\in \Omega=\{0<x<1,\quad 0<y<1\},
\end{equation}
 with the following integral and Dirichlet boundary conditions
\begin{align}
 %  \dfrac{\partial^2 u}{\partial x^2}+ \dfrac{\partial^2 u}{\partial y^2}=&f(x, y), \quad (x, y)\in D=\{0<x,y<1\}, \nonumber\\
	&   \int_0^{\xi_1} u(x, y)dx =\nu_a(y),  \quad \int_{\xi_2}^1 u(x, y)dx = \nu_b(y), \label{eq:1.2}\\
	& u(x, 0) = \nu_c(x), \quad u(x, 1) = \nu_d(x), \label{eq:1.3}
   \end{align}
where $0<\xi_1<1$, \ $0<\xi_2<1$.  Let us consider, all three cases that are possible: $\xi_1<\xi_2$, $\xi_1=\xi_2$ and $\xi_1>\xi_2$.

The purpose of this article is to study the Peaceman-Rachford Alternating Direction Implicit (ADI) method for a system of difference equations approximating a differential problem \eqref{eq:1.1}--\eqref{eq:1.3}. To the best of the authors' knowledge,
 %(As far as the authors know),
 iterative methods for the system of difference equations in the case of nonlocal conditions \eqref{eq:1.2} have not been studied.

To examine the convergence conditions for the ADI method, we analyze in sufficient detail the structure of the spectrum of the corresponding  difference
problem. The structure of the spectrum for a differential problem with other types of nonlocal conditions is considered in many papers (see, for example,
\cite{2002_Sapagovas,  2009LMJ49n3o,  2011DE47n8,2008NA13n4,2007MMA12n2,2022LMJ62n4}. 
As in previous our papers \cite{2009LMJ49n3o, 2011DE47n8, 2014NA19n3o}, 
      we proved the convergence of ADI method in the case of the system of eigenvectors is complete. But in the present paper we take some comments and examples about the convergence of ADI method without this condition.

The further structure of this paper is as follows.
The difference problem corresponding to the differential problem is formulated in Section \ref{s:2}, where the ADI method is also introduced.
The structure of the spectrum of difference problems is discussed in Section~\ref{s:3}. The convergence of the ADI method is demonstrated in Section~\ref{s:4}. In Section~\ref{s:5}, a higher order  finite difference method is considered.
In Section~\ref{s:6}, numerical results are provided to verify the accuracy and efficiency of the proposed algorithms.
 The last Section~\ref{s:7} presents comments and conclusions.


\section{Difference problem}\label{s:2}


Consider a uniform mesh  in $x$ and $y$ with   step size  $h=1/{N}$ ($1<N\in \mathbb{N}$):
\begin{equation*}
 \overline{\omega}^{h}_{x}:=\{x_i=ih, i=\overline{0,N}\}; \quad  \overline{\omega}^{h}_{y}:=\{y_j=jh, j=\overline{0,N}\},\quad \overline{\omega}^h:=\overline{\omega}^{h}_{x}\times\overline{\omega}^{h}_{y} .
\end{equation*}
We use the following notation
\begin{equation*}
\delta^2_{x}U_{i\,j}:=\frac{U_{i-1,j}-2U_{ij}+U_{i+1,j}}{h^2},\quad
\delta^2_{y}U_{i\,j}:=\frac{U_{i,j-1}-2U_{ij}+U_{i,j+1}}{h^2}.
\end{equation*}
Let us replace the differential problem \eqref{eq:1.1}--\eqref{eq:1.3} with the following difference problem  on the mesh $\overline{\omega}^h$
\begin{align}
 &-\delta^2_{x}U_{ij}-\delta^2_{y}U_{ij}={F}_{ij},\quad \quad &i,j=\overline{1, N-1}, \label{eq:2.1}\\
   \noalign{\smallskip}
  &h\Big( \frac{U_{0j}+U_{rj}}{2}+\sum_{i=1}^{r-1}U_{ij}\Big)=\mu^a_{j},\quad &j=\overline{1, N-1}, \label{eq:2.2}\\
   &h\Big( \frac{U_{sj}+U_{Nj}}{2}+\sum_{i=s+1}^{N-1}U_{ij}\Big)=\mu^b_{j},\quad &j=\overline{1, N-1}, \label{eq:2.3}\\
    &U_{i0}=\mu^c_{i}, \quad U_{iN}=\mu^d_{i}, \quad &i=\overline{1, N-1}, \label{eq:2.4}
  \end{align}
The integral conditions \eqref{eq:1.2} are approximated by the trapezoidal rule.
For simplicity, we assume that the values $\xi_1$ and $\xi_2$ are such that $\xi_1=rh$ and $\xi_2=sh$, \quad   $r, s \in \mathbb{N}$,   $ 0< r<N$, $0<s<N$. Assume that $N$, $r$ and $s$ are even numbers. Note, that it is not a strong restriction, as we can always halve the step size $h$.



The existence and uniqueness of the solution of the differential problem \eqref{eq:1.1}--\eqref{eq:1.3}  are investigated in \cite{2011_Avalishvili_Gordeziani,
2011_Avalishvili_Gordeziani_1}. The error  estimate and convergence of the solution of the finite difference method are presented in \cite{2013_Berikelashvili_Khomeriki, 2014_Berikelashvili_Khomeriki}.

The corresponding
difference scheme for this problem  under the condition that the desired solution belongs to the Sobolev space $W_2^s (1<s\le3)$ has been investigated in \cite{2014_Berikelashvili_Khomeriki}.

The system \eqref{eq:2.1}--\eqref{eq:2.4} has $N(N-1)$ equations \eqref{eq:2.1}--\eqref{eq:2.3}  and the same number of unknowns $U_{ij}$, $i=\overline{0, N-1}$, $j=\overline{1, N-1}$.
First of all, for  theoretical study we will write  the system \eqref{eq:2.1}--\eqref{eq:2.4} in a more compact  matrix form
\begin{equation}\label{eq:2.5}
\mathbf{A}_1\mathbf{U}+\mathbf{A}_2\mathbf{U}=\mathbf{\Phi},
\end{equation}
 where $\mathbf{A}_{1}, \mathbf{A}_{2}$ are $(N-1)^2$ order matrices and  \begin{align*}
  \mathbf{U}:=& (U_{11},\dots,U_{N-1,1},U_{12},\dots,U_{N-1,N-1})^T,\\
  \mathbf{\Phi}:=& (\Phi_{11},\dots,\Phi_{N-1,1},\Phi_{12},\dots,\Phi_{N-1,N-1})^T \
  \end{align*}
are $(N-1)^2$-vectors.  
 For this purpose, we will express for each $j=\overline{1, N-1}$ the unknowns $U_{0j}$ and $U_{Nj}$ from nonlocal conditions \eqref{eq:2.2}--\eqref{eq:2.3} through remaining unknowns $U_{ij}$, $i=\overline{1, N-1}$:
 \begin{align}
  U_{0j}=&-U_{rj}-2\sum_{i=1}^{r-1}U_{ij}+2h^{-1}\mu^a_{j},\quad &j=\overline{1, N-1}, \label{eq:2.6}\\
  U_{Nj}=&-U_{sj}-2\sum_{i=s+1}^{N-1}U_{ij}+2h^{-1}\mu^b_{j},\quad &j=\overline{1, N-1}. \label{eq:2.7}
   \end{align}
   Putting expressions \eqref{eq:2.6}--\eqref{eq:2.7} into Equation \eqref{eq:2.1} as $i=1$ or $i=N-1$,
 we get new system of equations, the order of which  and the number of the unknowns $U_{ij}$ , $i, j = \overline{1, N-1}$ are  equal to $(N-1)^2$:
\begin{align}
  &h^{-2}(2\sum_{i=1}^{r-1} U_{ij}+U_{rj}+2U_{1j}-U_{2j})-\delta^2_{y}U_{1j}=\Phi_{1j}:={F}_{1j} +2h^{-3}\mu^a_{j},
  \label{eq:2.8}\\
  &-\delta^2_{x}U_{ij}-\delta^2_{y}U_{ij}=\Phi_{ij}:={F}_{ij},\quad \quad i=\overline{2, N-2}, \nonumber
  \\
  &h^{-2}(-U_{N-2,j}+2U_{N-1,j}+U_{sj}+2\sum_{i=s+1}^{N-1}U_{ij})-\delta^2_{y }U_{N-1,j}\nonumber\\
  &=\Phi_{N-1,j}:={F}_{N-1,j} +2h^{-3}\mu^b_{j},\label{eq:2.10}
\end{align}
where  $j=\overline{1, N-1}$.
The system \eqref{eq:2.8}--\eqref{eq:2.10} together with \eqref{eq:2.6}--\eqref{eq:2.7} and boundary conditions \eqref{eq:2.4} is equivalent to the system \eqref{eq:2.1}--\eqref{eq:2.4}.

 Next, we will consider the system \eqref{eq:2.8}--\eqref{eq:2.10}, \eqref{eq:2.4}. We will write  the system \eqref{eq:2.1}--\eqref{eq:2.4} in the  matrix form \eqref{eq:2.5}. For this purpose let us define matrices of order $(N-1)$:
\begin{align}\label{eq:2.11}
\mathbf{\Lambda}=&\frac{1}{h^{2}}\left(
                                           \begin{array}{ccccc}
                                             2 & -1 &  &  &  \\
                                             -1 & 2 & -1 & \cdots & 0 \\
                                              & \ddots & \ddots & \ddots &  \\
                                              &  & -1 & 2 & -1 \\
                                              &  &  & -1 & 2 \\
                                           \end{array}
                                         \right), \\
 \mathbf{C}=&\frac{1}{h^{2}}
    \bordermatrix[{()}]{%
     & \mbox{\scriptsize{1}} & \mbox{\scriptsize{ 2}} & \cdots   &\mbox{\scriptsize{$ r-1$}} & \mbox{\scriptsize{$r$}} & \mbox{\scriptsize{$r+1$}} & \cdots   & \mbox{\scriptsize{$s-1$}} & \mbox{\scriptsize{$s$}} & \mbox{\scriptsize{$s+1$}} & \cdots   & \mbox{\scriptsize{$ N-1$ }} \cr
      & 2 & 2 & \cdots & 2 & 1 & 0 & \cdots & 0 & 0 & 0 & \cdots & 0 \cr
  & 0 & 0 & \cdots& 0 & 0 & 0 & \cdots & 0 & 0 &0 & \cdots & 0 \cr
   &   \vdots  &  \vdots  & \cdots &  \vdots  &  \vdots  &  \vdots  &  \cdots  &  \vdots  &  \vdots  &  \vdots  &  \cdots  &  \vdots \cr
     & 0 & 0 & \cdots& 0 & 0 & 0 & \cdots & 0 & 0 &0 & \cdots & 0 \cr
  & 0 & 0 & \cdots & 0 & 0 & 0 & \cdots & 0 & 1 & 2 & \cdots & 2 \cr
    }  .
  \nonumber
  \end{align}
Note that only the first and the last row of matrix $ \mathbf{C}$ are non-zero.
Here we indicate the column numbers on top of the matrix. Let us denote
\[
\mathbf{\Lambda}_x:=\mathbf{\Lambda}+\mathbf{C}, \quad \quad \mathbf{\Lambda}_y:=\mathbf{\Lambda}.
\]
 The system \eqref{eq:2.1}--\eqref{eq:2.4} can be written in   matrix form \eqref{eq:2.5} using matrices
\begin{align}
\mathbf{A}_1&:=\mathbf{I}_y\otimes\mathbf{\Lambda}_x=\!\left( \!\!
                \begin{array}{cccc}
                  \mathbf{\Lambda}_{x} &  &  &  \\
                   & \mathbf\Lambda_{x} &  &  \\
                   &  & \ddots &  \\
                   &  &  & \mathbf\Lambda_{x} \\
                \end{array}
              \!\!
              \right)\!,\label{eq:2.12}  \\
\noalign{\smallskip}
\mathbf{A}_2&:=\mathbf{\Lambda}_y\otimes \mathbf{I}_x=\left(
                                       \begin{array}{ccccc}
                                         2\mathbf{I}_y  & -\mathbf{I}_y  &  &  &  \\
                                         -\mathbf{I}_y  & 2\mathbf{I}_y  & -\mathbf{I}_y  &  &  \\
                                          & \ddots & \ddots & \ddots &  \\
                                          &  & -\mathbf{I}_y  & 2\mathbf{I}_y  & -\mathbf{I}_y  \\
                                          &  &  & -\mathbf{I}_y & 2\mathbf{I}_y  \\
                                       \end{array}
                                     \right),\nonumber%\label{eq:2.13}
\end{align}
where  $\mathbf{I}_x$ and $\mathbf{I}_y$ are the identity matrices of order $(N-1)$ (in our case of square domain $\mathbf{I}_x=\mathbf{I}_y$).

Our main goal is to study the ADI method for solving a system of difference equations.
We write the ADI method for system \eqref{eq:2.5}:
\begin{align}\label{eq:2.14}
  \frac{\mathbf{U}^{n+1/2}-\mathbf{U}^{n}}{\tau^1_{n}}+ {\mathbf{A}}_1 \mathbf{U}^{n+1/2}+
   {\mathbf{A}}_2 \mathbf{U}^{n}=&{\mathbf{\Phi}},\nonumber\\
  \noalign{\smallskip}
  \frac{\mathbf{U}^{n+1}-\mathbf{U}^{n+1/2}}{\tau^2_{n}}
  +{\mathbf{A}}_1\mathbf{U}^{n+1/2}+{\mathbf{A}}_2\mathbf{U}^{n+1}=&{\mathbf{\Phi}}.
\end{align}
where $\tau^1_{n}$, $\tau^2_{n}$,   $n=1,\dots $ are iteration parameters.  We give a explicit formula for determining iteration parameters in  Section~\ref{s:4}.



\section{The structure of the spectrum of the difference problem}\label{s:3}


Proof of the convergence of the method \eqref{eq:2.14} is based on the structure of the spectrum of the one-dimensional eigenvalue problems.
We consider two difference eigenvalue problems. First of these is problem  with nonlocal boundary conditions:
\begin{align}
-\frac{v_{i-1}-2v_{i}+v_{i+1}}{h^2}=&\eta v_{i}, \quad i=\overline{1, N-1}, \label{eq:3.1}\\
\frac {v_0+v_r}2+\sum \limits _{i=1}^{r-1}v_i=&0, \label{eq:3.2}\\
\frac {v_s+v_N}2+\sum \limits _{i=s+1}^{N-1}v_i=&0. \label{eq:3.3}
 \end{align}
The second problem is classical problem:
 \begin{align}
&-\frac{w_{j-1}-2w_{j}+w_{j+1}}{h^2}=\mu w_{j}, \quad j=\overline{1, N-1}, \label{eq:3.4}\\
&w_0=0, \quad w_N=0. \label{eq:3.5}
 \end{align}

In problem with nonlocal  conditions we express  $v_{0}$ and $v_{N}$ from  \eqref{eq:3.2}--\eqref{eq:3.3}  and substitute into \eqref{eq:3.1}. So we get
\begin{equation}\label{eq:3.6}
\mathbf\Lambda_x \mathbf{v}=\eta \mathbf{v},\quad \mathbf{v}=(v_1,\dots,v_{N-1})^T.
\end{equation}
We can write the problem \eqref{eq:3.4}--\eqref{eq:3.5} in matrix form
\begin{equation}\label{eq:3.7}
{\mathbf\Lambda}_y \mathbf{w}=\mu \mathbf{w},\quad \mathbf{w}=(w_1,\dots,w_{N-1})^T.
\end{equation}
The problem  \eqref{eq:3.7} has known solution
\begin{align}
\mu_l=&\frac{4}{h^2}\sin^2\frac{\pi  l h}{2}, &l=\overline{1, N-1}, \nonumber\\%\label{eq:3.8}\\
\mathbf{w}^l=&(w^l_1,\dots,w^l_{N-1})^T, \quad  w^l_j=\sin(\pi ljh), &  j,l=\overline{1, N-1}. \label{eq:3.9}
 \end{align}
%
% \[
% \mu_{kl}=\eta_k+\theta_l, \quad k,l=\overline{1, N-1}.
% \]
We will find eigenvalues and eigenvectors of the problem \eqref{eq:3.6}. As far as we know, this problem has not been investigated.
\begin{lem}\label{lem:1}
The eigenvalues of the problem \eqref{eq:3.6}, which satisfy conditions
\begin{equation*}
0<\eta_k<\frac{4}{h^2},
\end{equation*}
 can be expressed in the form
\begin{equation}\label{eq:3.10}
\eta_k=\frac{4}{h^2}\sin^2\frac{\alpha_k h}{2},\ k=\overline{1, N-3},
\end{equation}
where $\alpha_k$ is root of any  of equations
%\begin{align*}
\[
\sin  \frac{\alpha(1-\xi_2)}{2}=0,\quad
\sin  \frac{\alpha \xi_1}{2}=0,\quad
\sin  \frac{\alpha(1-\xi_1+\xi_2)}{2}=0.
\]
% \end{align*}
Eigenvectors $\bm{v}^k$ can be expressed as
 \begin{equation}\label{eq:3.11}
v^k_i=c_1\cos (\alpha_k ih) + c_2 \sin (\alpha_k ih),
\end{equation}
where $(c_1, c_2)$ is nontrivial solution of system
\begin{align}
&c_1h\frac{\sin(\alpha \xi_1)}{2\tan \frac{\alpha h}{2}}+c_2h\frac{1-\cos(\alpha \xi_1)}{2\tan \frac{\alpha h}{2}}=0, \label{eq:3.12}\\
&c_1h\frac{\sin\alpha -\sin( \alpha \xi_2)}{2\tan \frac{\alpha h}{2}}+c_2h\frac{\cos(\alpha \xi_2)-\cos\alpha}{2\tan \frac{\alpha h}{2}}=0. \nonumber
 \end{align}
 in the case $\alpha=\alpha_k$.
%$$ c_1^2+c_2^2\neq 0.$$  %\ l=\overline{1, N-1}.$$
\end{lem}
\begin{pf}
The inequality
\[
\abs{ 1-\frac{\eta h^2}{2} }<1
\]
is satisfied for $0<\eta<4/h^2$ in Equation \eqref{eq:3.1}.
So let us denote
\[
\cos (\alpha h)=1-\frac{\eta h^2}{2}, \quad 0<\alpha<\frac{\pi}{h},
\]
where $\alpha$ is  a new parameter instead of $\eta$.
Then general solution of Equation~\eqref{eq:3.1}  can be written in form of  \eqref{eq:3.9}.
This solution must obey conditions \eqref{eq:3.2} and \eqref{eq:3.3}. Putting the expression \eqref{eq:3.11} into \eqref{eq:3.2}  and \eqref{eq:3.3} after some elementary rearragements   we get \eqref{eq:3.12}.

Eigenvectors are defined by \eqref{eq:3.11}. So, $v_i\not\equiv 0$, namely $c_1^2+c_2^2\neq 0$. It means that  determinant $D$ of the system \eqref{eq:3.12} must be equal to zero.
After some calculations we have
\[
D=4\sin \frac{\alpha(1-\xi_1+\xi_2)}{2}\sin\frac{\alpha \xi_1}{2}\sin \frac{\alpha(1-\xi_1)}{2}=0.
\]
So, the statement of Lemma follows from here.


For even numbers $r$, $s$ and $N$ ($rh=\xi_1$, $sh=\xi_2$, $Nh=1$)
\begin{description}
  \item[-] equation $\sin \frac{\alpha \xi_1}{2}=0$ has $\frac{r}{2}-1$ different roots,
  \item[-] equation $\sin \frac{\alpha (1- \xi_2)}{2}=0$ has $\frac{N}{2}-\frac{s}{2}-1$ different roots,
  \item[-] equation $\sin \frac{\alpha (1- \xi_1+\xi_2)}{2}=0$ has $\frac{N}{2}-\frac{r}{2}+\frac{s}{2}-1$ different roots.
\end{description}
There are $N-3$ roots.
\begin{align}
  & \alpha_k=\frac{2\pi k}{\xi_1},  &k=1,2,\dots, \frac{r}{2}-1,  \label{eq:3.13}\\
  & \alpha_k=\frac{2\pi k}{1-\xi_2},  &k=1,2,\dots, \frac{N}{2}-\frac{s}{2}-1,  \label{eq:3.14}\\
   & \alpha_k=\frac{2\pi k}{1-\xi_1+\xi_2},  &k=1,2,\dots, \frac{N}{2}-\frac{r}{2}+\frac{s}{2}-1.  \label{eq:3.15}
\end{align}
 \end{pf}

%\end{corollary}
\begin{rem}
In each of the formulas \eqref{eq:3.13}--\eqref{eq:3.15}, the roots are different.
 However, the eigenvalue \eqref{eq:3.10}
\begin{equation*}%\label{eq17}
\eta_k=\frac{4}{h^2}\sin^2\frac{\alpha_k h}{2},\ k=\overline{1, N-3},
\end{equation*}
 can be multiple if some of the roots of different formulas coincide.
 For example, in the case $\xi_1=1-\xi_2$ the equations \eqref{eq:3.13} and \eqref{eq:3.14} have the same roots. In this case, the eigenvalues are multiples, and the system of eigenvectors may not be complete. 
\end{rem}

\begin{lem}\label{lem:2}
In the case of even  numbers $r$, $s$ and $N$
\begin{equation}\label{eq:3.16}
\eta=\frac{4}{h^2}
\end{equation}
is the eigenvalue of the problem \eqref{eq:3.6} of multiplicity  2  %is twice multiple eigenvalue
with  two  linearly independent corresponding eigenvectors $\mathbf{v}^1$ and $\mathbf{v}^2$:
\begin{equation}
v_i^1=(-1)^i ih,\quad  v_i^2=(-1)^i,\quad i=\overline{1, N-1}. \label{eq:3.17}
 \end{equation}
\end{lem}
\begin{pf}
The general solution of Equation \eqref{eq:3.1} in the case $\eta=4/h^2$ is

\begin{equation}\label{eq:3.18}
 v_i=(-1)^i(c_1+c_2ih),\quad i=\overline{1, N-1}.
\end{equation}
Putting this expression into conditions \eqref{eq:3.2} and  \eqref{eq:3.3} we have that these conditions are satisfied with all values of $c_1$ and $c_2$, if
$N$, $r$ and $s$ are even. So, it is possible to choose constants  $c_1$ and $c_2$ such that two linearly independent eigenvectors are defined by formula
\eqref{eq:3.18}, corresponding to eigenvalue  \eqref{eq:3.16}. In particularly, those are vectors  \eqref{eq:3.17}.  
\end{pf}
%
\begin{cor}\label{cor:1}
  It follows from  Lemmas \ref{lem:1} and \ref{lem:2}  that all $N-1$  eigenvalues of the difference eigenvalue problem \eqref{eq:3.6}  are positive. Depending on $\xi_1$ and $\xi_2$  the system of eigenvectors may be complete or not.
\end{cor}
%

\section{The convergence of the ADI method}\label{s:4}


Let us write the iterative method \eqref{eq:2.14}  as a matrix equation
\begin{equation}\label{eq:4.01}
  \mathbf{U}^{n+1}=\mathbf{S}_{n}\mathbf{U}^{n}+{\mathbf{\Psi}},
  \end{equation}
where
%
\begin{equation}\label{eq:4.02}
  \mathbf{S}_n=(\mathbf{I}+\tau^2_{n}{\mathbf{A}}_2)^{-1}(\mathbf{I}-\tau^2_{n}{\mathbf{A}}_1)(\mathbf{I}+ \tau^1_{n}{\mathbf{A}}_1)^{-1}(\mathbf{I}- \tau^1_{n}{\mathbf{A}}_2),
\end{equation}
where  $\mathbf{I}=\mathbf{I}_y\otimes \mathbf{I}_x$ is the identity matrix of order $(N-1)^2$.
We will prove that spectral radius $\varrho(\mathbf{S}_n)<1$.



\begin{lem}\label{lem:3}
$\mathbf{A}_1$ and  $\mathbf{A}_2$ are commuting matrices
\begin{equation*}
 \mathbf{A}_1 \mathbf{A}_2=\mathbf{A}_2 \mathbf{A}_1.
\end{equation*}
\end{lem}

\begin{pf}
It is easy to check that
\begin{equation*}
\mathbf{A}_1 \mathbf{A}_2=\mathbf{A}_2 \mathbf{A}_1={\mathbf\Lambda}_y\otimes {\mathbf\Lambda}_x=\frac{1}{h^2}\left (\begin{array}{ccccc}
    2{\mathbf\Lambda}_{x} & -{\mathbf\Lambda}_{x} &  & \mathrm{O}& \\
   -{\mathbf\Lambda}_{x} & 2{\mathbf\Lambda}_{x} &-{\mathbf\Lambda}_{x} &  & \\
     &  \ddots& \ddots & \ddots & \\
     &  &-{\mathbf\Lambda}_{x} &2{\mathbf\Lambda}_{x} & -{\mathbf\Lambda}_{x} \\
     & \mathrm{O} & &  -{\mathbf\Lambda}_{x} &2{\mathbf\Lambda}_{x}
 \end{array}\right ).
\end{equation*}
Using eigenvectors $\mathbf{v}^k$, $\mathbf{w}^l,\  k, l,=\overline{1, N-1}$ of the eigenvalue problems \eqref{eq:3.1}--\eqref{eq:3.3} and \eqref{eq:3.4}--\eqref{eq:3.5}
we define vectors
\begin{equation}\label{eq:4.1}
\mathbf{U}^{kl}=\mathbf{w}^l\otimes \mathbf{v}^k, \quad\text{where}\ U^{kl}_{ij}=v^k_iw^l_j,\quad  i, j=\overline{1, N-1}.   \qquad 
\end{equation}
\end{pf}
\begin{lem}\label{lem:4}
Matrices $\mathbf{A}_1$,   $\mathbf{A}_2$, $\mathbf{A}_1+ \mathbf{A}_2$, $\mathbf{A}_1 \mathbf{A}_2$ and $\mathbf{A}_2 \mathbf{A}_1$  have the same system of eigenvectors (not necessarily full) $\mathbf{U}^{kl}=\mathbf{w}^l\otimes \mathbf{v}^k$.
\end{lem}

\emph{Proof.}
It follows from definition of $\mathbf{A}_1$ and  $\mathbf{U}^{kl}$ by formulas \eqref{eq:2.12} and \eqref{eq:4.1} that
\begin{equation}\label{eq:4.2}
\mathbf{A}_1 \mathbf{U}^{kl}=\eta_k \mathbf{U}^{kl}.
\end{equation}
Furthermore, $\eta$ and $\mathbf{v}$ are eigenvalue and eigenvector of matrix ${\mathbf\Lambda}_x$; $\mu$ and $\mathbf{w}$ are eigenvalue and eigenvector of matrix ${\mathbf\Lambda}_y$. So, from properties of tensor product we obtain that

\begin{equation}\label{eq:4.3}
(\mathbf{A}_1+\mathbf{A}_2) \mathbf{U}^{kl}=(\mathbf{I}_y\otimes{\mathbf\Lambda}_x+{\mathbf\Lambda}_y\otimes \mathbf{I}_x)(\mathbf{w}^l\otimes \mathbf{v}^k)=(\mu_k+\eta_l) \mathbf{U}^{kl}.
\end{equation}
According to Lemma \ref{lem:1}
\begin{equation*}
(\mathbf{A}_1\mathbf{A}_2) \mathbf{U}^{kl}=(\mathbf{A}_2\mathbf{A}_1) \mathbf{U}^{kl}=({\mathbf\Lambda}_y\otimes{\mathbf\Lambda}_x) \mathbf{U}^{kl}
=\mu_k\eta_l \mathbf{U}^{kl}.
 \end{equation*}
 It follows from \eqref{eq:4.2} and \eqref{eq:4.3} that
 \begin{equation*}
\mathbf{A}_2 \mathbf{U}^{kl}=\eta_l \mathbf{U}^{kl}. \quad \qed
\end{equation*}

\begin{rem}\label{rem:1}
We note that the system of eigenvectors $\mathbf{w}^l, l=1, \dots, N-1$ is  always complete. So, if the system of eigenvectors $\mathbf{v}^k, k=1, \dots, N-1$ is complete, then the system of eigenvectors  $\mathbf{U}^{kl}=\mathbf{w}^l\otimes \mathbf{v}^k$ is complete, too.     
    It is follows from properties of tensor product.
\end{rem}

Now, we can prove the statement on convergence of the iterative method.
\begin{thm}\label{thm:1}
If $\tau^1_n=\tau^2_n=:\tau_n$, where $\tau_n\in [\beta_1, \beta_2]$, $\beta_1>0$, $\beta_2<\infty$, and the system of eigenvectors of the problem \eqref{eq:3.6} is complete,  then   ADI method converges.
\end{thm}
\emph{Proof.}
Let be $\mathbf{U}^*$ the solution of the difference problem \eqref{eq:2.1}--\eqref{eq:2.4}. We denote the error of iterative method as
\[
\mathbf{Z}^n:=\mathbf{U}^*-\mathbf{U}^n.
\]  
From \eqref{eq:4.01} and Remark \ref{rem:1} follows that
\[
\mathbf{Z}^{n+1}=\mathbf{S}_n\mathbf{Z}^n=\prod_{j=0}^n\mathbf{S}_{n-j}
\mathbf{Z}^0=\sum_{k, l=1}^{N-1}c_{kl}\cdot \prod_{j=0}^n\mathbf{S}_{n-j}\mathbf{U}^{kl}.
\]
Hence, for any vector norm we have
\begin{equation}\label{4.4}
\lVert \mathbf{Z}^{n+1}\rVert \le\varrho\Big(\prod_{j=0}^n \mathbf{S}_{n-j}\Big) \sum_{k, l=1}^{N-1}|c_{kl}| \lVert\mathbf{U}^{kl}\rVert.
\end{equation}
Let us estimate the factor $ \varrho(\prod_{j=0}^n \mathbf{S}_{n-j})$ and prove that it tends to 0 when $n\rightarrow\infty $.  
It means that the ADI method converges.

Since $\mathbf{A}_1$ and $\mathbf{A}_2$ commute and have the same system of eigenvectors, then
\[
\lambda\Big(\prod_{j=0}^n \mathbf{S}_{n-j}\Big)=\prod_{j=0}^n \lambda(\mathbf{S}_{n-j}), \qquad \varrho\Big(\prod_{j=0}^n \mathbf{S}_{n-j}\Big)=\prod_{j=0}^n \varrho(\mathbf{S}_{n-j}).
\]

From Lemma \ref{lem:4}  and \eqref{eq:4.02}, it follows that an  eigenvalue of matrix $\mathbf{S}_n$
\begin{equation}\label{4.5}
\lambda(\mathbf{S}_n)=\frac{\big(1-\tau_n\lambda(\mathbf{A}_1)\big)\big(1-\tau_n\lambda(\mathbf{A}_2)\big)}
{\big(1+\tau_n\lambda(\mathbf{A}_1)\big)\big(1+\tau_n\lambda(\mathbf{A}_2)\big)}.
\end{equation}
If taking concrete values $\tau_n$ and $\lambda(\mathbf{A}_1)$,  the inequality $1-\tau_n\lambda(\mathbf{A}_1)\ge0$ is true, then
\[
\abs{\frac{1-\tau_n\lambda(\mathbf{A}_1)}{1+\tau_n\lambda(\mathbf{A}_1)}}=\frac{1-\tau_n\lambda(\mathbf{A}_1)}{1+\tau_n\lambda(\mathbf{A}_1)} \le \frac{1-\beta_1\lambda(\mathbf{A}_1)}{1+\beta_1\lambda(\mathbf{A}_1)}
=\varrho_1<1, 
\]
here $\varrho_1$ depends only on $\beta_1$ and $\lambda(\mathbf{A}_1)$ but does not depend on $n$. 
Analogously, if $1-\tau_n\lambda(\mathbf{A}_1)\le0$ then
\[
\abs{\frac{1-\tau_n\lambda(\mathbf{A}_1)}{1+\tau_n\lambda(\mathbf{A}_1)} }=\frac{\tau_n\lambda(\mathbf{A}_1)-1}{1+\tau_n\lambda(\mathbf{A}_1)} \le \frac{\beta_2\lambda(\mathbf{A}_1)-1}{1+\beta_2 \lambda(\mathbf{A}_1)}
=\varrho_2<1, 
\]
where $\varrho_2$ depends only on $\beta_2$ and $\lambda(\mathbf{A}_2)$ but does not depend on $n$.  

The second factor in the formula \eqref{4.5} is estimated similarly
\[
\abs{\frac{1-\tau_n\lambda(\mathbf{A}_2)}{1+\tau_n\lambda(\mathbf{A}_2)} } \le \varrho_0<1,
\]
where $\varrho_0$ depends  on $\beta_1$,  $\beta_2$ and $\lambda(\mathbf{A}_2)$ but does not depend on $n$.  
Finally we get from \eqref{4.5}
\[
\varrho(\mathbf{S}_{n})=\max |\lambda(\mathbf{S}_{n})|\le \tilde{\varrho}<1,
\]
where $\tilde{\varrho}$ depends  on $\beta_1$,  $\beta_2$, $\lambda(\mathbf{A}_1)$ and $\lambda(\mathbf{A}_2)$ but does not depend on $n$. Because $\tilde{\varrho}<1$ then 
\[
\prod_{j=0}^n\varrho (\mathbf{S}_{n-j})\le \Tilde{\varrho}^{n+1} \rightarrow 0, \quad  n \rightarrow \infty.   \qquad {\qed}
\]
\begin{rem}\label{rem:2}
If parameters of ADI method $\tau_1$ and $\tau_2$ does not depend on $n$, then Theorem \ref{thm:1} is  correct even without the assumption of completeness of the system of eigenvectors. We have that when $\tau_n^1=\tau_n^2\equiv \tau>0$,  then the iterative method is stationary
\begin{equation}\label{4.5a}
 \mathbf{U}^{n +1}=\mathbf{S}\mathbf{U}^n+{\mathbf{\Psi}},
\end{equation}
where
%
\begin{equation*}
  \mathbf{S}=(\mathbf{I}+\tau{\mathbf{A}}_2)^{-1}(\mathbf{I}-
  \tau{\mathbf{A}}_1)(\mathbf{I}+\tau{\mathbf{A}}_1)^{-1}
  (\mathbf{I}- \tau{\mathbf{A}}_2).
\end{equation*}
A necessary and sufficient condition for the convergence of a stationary process \eqref{4.5a} with any initial data $\mathbf{U}^{0}$ is  $\varrho(\mathbf{S})\le 1$. 
\end{rem}

It is known that when applying the ADI method, optimal parameters $\tau^1_n$ and $\tau^2_n$ are usually used.


For obtaining optimal set of iteration parameters in the case of commuting
operators we suppose that  eigenvalues of ${\mathbf{A}}_1$ and ${\mathbf{A}}_2$ satisfy inequalities
\begin{equation*}
 \delta_1 \leq  \lambda({\mathbf{A}}_1) \leq \Delta_1, \quad
 \delta_2 \leq  \lambda({\mathbf{A}}_2) \leq \Delta_2,
\end{equation*}

 Following \cite[Ch.X, \S4]{45}, the set of ADI parameters by Jordan gives
 \begin{equation*}
\tau_1^n=\dfrac{qw_n+\varkappa}{1+w_np}, \quad  \tau_2^n=\dfrac{qw_n-\varkappa}{1-w_np}, \quad  n=\overline{1, m},
\end{equation*}
 where $m$ is number of iterations
 \begin{equation*}
m=m(\varepsilon)\approx \dfrac{1}{\pi^2} {\ln} \dfrac{4}{\varepsilon} \ln\dfrac{4}{\eta},
\end{equation*}
and we use notation
\begin{align*}
&t=\sqrt{\dfrac{(\Delta_1-\delta_1)(\Delta_2-\delta_2)} {(\Delta_1+\delta_2)(\Delta_2+\delta_1)}}, \quad
\kappa=\dfrac{(\Delta_1-\delta_1)\Delta_2}{(\Delta_2+\delta_1)\Delta_1}, \quad
\varkappa=\dfrac{\Delta_1-\Delta_2+(\Delta_1+\Delta_2)p}{2\Delta_1\Delta_2}, \\
&\eta=\dfrac{1-t}{1+t}, \quad p=\dfrac{\kappa-t}{\kappa+t}, \quad q=\varkappa+\dfrac{1-p}{\Delta_1}
\end{align*}
and
\[
w_n=\dfrac{(1+2\theta)(1+\theta^{\sigma})}  {2\theta^{\sigma/2}(1+\theta^{1-\sigma}+\theta^{1+\sigma})}, \quad n=\overline{1, m},
\]
where 
\begin{equation*}
\theta:=\dfrac{1}{16}\eta^2(1+\dfrac{1}{2}\eta^2), \quad
\sigma:=\sigma_n=\dfrac{2n-1}{2m}, \quad n=\overline{1, m}.
\end{equation*}


Such algorithm for the optimal choice of parameters $\tau^1_n$ and $\tau^2_n$ is proposed for the case of symmetric matrices ${\mathbf{A}}_1$ and ${\mathbf{A}}_2$, that is, when the system of eigenvectors for the finite difference scheme is complete. 

Let's briefly discuss what such a choice of parameters means in the case when the system of eigenvectors is incomplete (for problem \eqref{eq:2.1}--\eqref{eq:2.4}.  

From the definition of the ADI method \eqref{eq:2.14}, it follows that 
\begin{equation}\label{4.6}
\mathbf{Z}^{n}= \mathbf{S}_{n-1}\mathbf{Z}^{n-1}= \prod_{j=0}^{n-1} \mathbf{S}_{n-1-j}\mathbf{Z}^{0},
\end{equation}
where $\mathbf{Z}^{n}=\mathbf{U}^*-\mathbf{U}^n$ and $\mathbf{U}^*$ is the exact solution of the finite difference problem.

Without details, we note that the choice of optimal parameters $\tau^1_n$ and $\tau^2_n$  is based on the solution of the minimax problem, namely the finding of the minimum value of the spectral radius of the matrix $\prod_{j=0}^{n-1} \mathbf{S}_{n-1-j}$. 
Since for a symmetric matrix the spectral radius can be considered as the norm of the matrix, then in the case of symmetric matrices ${\mathbf{A}}_1$ and ${\mathbf{A}}_2$, it follows from {\color{blue}\eqref{4.6}}
\begin{equation}\label{4.7}
\norm{\mathbf{Z}^{n}}\le  \varrho\Big(\prod_{j=0}^{n-1} \mathbf{S}_{n-1-j}\Big)\norm{\mathbf{Z}^{0}}\rightarrow 0, \quad  n \rightarrow \infty. 
\end{equation}
However, if the system of eigenvectors of the finite difference problem \eqref{eq:2.1}--\eqref{eq:2.4} is incomplete, then \eqref{4.7} does not follow from \eqref{4.6}.
In this case, we can use the following statement from linear algebra \cite{29}.

Let $\mathbf{A}$ be an arbitrary square matrix.
If $\varepsilon>0$ is given, then there is matrix norm $\norm{ \mathbf{A}}_*$ such that 
\[
\norm{ \mathbf{A}}_*\le \varrho(\mathbf{A})+\varepsilon.
\]
It means, particularly, that from inequality $\varrho(\mathbf{A})+\varepsilon <1$ follows $\norm{ \mathbf{A}}_*<1$. 



These arguments, in our opinion, are not proof of the convergence of the ADI method with optimal parameters. But, at least, this is a strong motivation for the corresponding numerical experiment with optimal parameters in the case when the system of eigenvectors is incomplete.

\section{The higher-order method }\label{s:5}
Now, we will consider the difference problem approximating the differential problem \eqref{eq:1.1}--\eqref{eq:1.3} with truncation error $\mathcal{O}(h^4)$.
So, let us replace the differential problem with the following difference problem
\begin{align}
 &-(\delta^2_{x}U_{ij}+\delta^2_{y}U_{ij}+\frac{h^2}{6}\delta^2_{x}\delta^2_{y}U_{ij})=\tilde{F}_{ij},\nonumber\\
 &\tilde{F}_{ij}=F_{ij}+\frac{h^2}{12}\big(\delta^2_{x}+\delta^2_{y}\big)F_{ij},  \quad i,j=\overline{1, N-1},   \label{eq:5.1}\\
   \noalign{\smallskip}
  &\frac{h}{3}\Big( {U_{0j}+U_{rj}}+4\sum_{i=1}^{{r}/{2}}U_{2i-1, j}+2\sum_{i=1}^{{r}/{2}-1}u_{2i, j}\Big)=\mu^a_{j},\quad j=\overline{1, N-1}, \label{eq:5.2}\\
   &\frac{h}{3}\Big( {U_{sj}+U_{Nj}}+4\sum_{i=s/2+1}^{{N}/{2}}U_{2i-1, j}+2\sum_{i=s/2+1}^{{N}/{2}-1}U_{2i, j}\Big)=\mu^b_{j}, j=\overline{1, N-1}, \label{eq:5.3}\\
    &U_{i0}=\mu^c_{i}, \quad U_{iN}=\mu^d_{i}, \quad i=\overline{1, N-1}, \label{eq:5.4}
  \end{align}
In this case one dimensional eigenvalue problem with nonlocal boundary conditions is
\begin{align}
&-\frac{v_{i-1}-2v_{i}+v_{i+1}}{h^2}=\eta v_{i}, \quad i=\overline{1, N-1}, \label{eq:5.5}\\
&\frac{h}{3}\Big( {v_{0}+v_{r}}+4\sum_{i=1}^{\frac{r}{2}}v_{2i-1}+2\sum_{i=1}^{\frac{r}{2}-1}v_{2i}\Big)=0, \label{eq:5.6}\\
&\frac{h}{3}\Big( {v_{s}+v_{N}}+4\sum_{i=s+1}^{\frac{N}{2}}v_{2i-1}+2\sum_{i=s+1}^{\frac{N}{2}-1}v_{2i}\Big)=0. \label{eq:5.7}
 \end{align}
We rewrite the eigenvalue problem  \eqref{eq:5.5}--\eqref{eq:5.7} in equivalent matrix form. For this we express the values $v_0$ ir $v_N$ from the condition
 \eqref{eq:5.6}--\eqref{eq:5.7}. Putting these expressions into Equation  \eqref{eq:5.5} we can rewrite problem  \eqref{eq:5.5}--\eqref{eq:5.7} as follows
 \begin{equation}\label{eq:5.8}
{\mathbf\Lambda}_x \mathbf{v}=\eta \mathbf{v}.
\end{equation}
 Similarly, as in Section~\ref{s:2}, we have
 %\[
 ${\mathbf\Lambda}_x={\mathbf\Lambda}+\mathbf{C}/h^2$,
 %\]
 where ${\mathbf\Lambda}$ is  defined in \eqref{eq:2.11} and $\mathbf{C}$ is equal to
 \begin{equation*}
  \!\!\!\!\bordermatrix[{()}]{%
     & \mbox{\scriptsize{1}} & \mbox{\scriptsize{ 2}} & \mbox{\scriptsize{ 3}} & \cdots   &\mbox{\scriptsize{$ r-1$}} & \mbox{\scriptsize{$r$}} &
     \mbox{\scriptsize{$r+1$}} & \cdots   & \mbox{\scriptsize{$s-1$}} & \mbox{\scriptsize{$s$}} & \mbox{\scriptsize{$s+1$}} & \cdots   & \mbox{\scriptsize{$
     N-2$ }} & \mbox{\scriptsize{$ N-1$ }} \cr
   & 4 & 2 & 4 & \cdots & 4 & 1 & 0 & \cdots & 0 & 0 & 0 & \cdots & 0  &0\cr
  & \vdots & \vdots & \vdots & \cdots & \vdots & \vdots & \vdots & \cdots & \vdots & \vdots &\vdots &\cdots & \vdots & \vdots  \cr
   & 0 & 0 & 0& \cdots & 0 & 0 & 0 & \cdots & 0 & 1 & 4 & \cdots & 2& 4 \cr
    } .
 \end{equation*}
  We  note, that two eigenvalue problems \eqref{eq:5.5}--\eqref{eq:5.7} and  \eqref{eq:5.8} are equivalent.
\begin{lem}\label{lem:5}
The eigenvalues of the problem \eqref{eq:5.5}--\eqref{eq:5.7}, which satisfy inequality
%\[
$0<\eta_k<{4}/{h^2}$
%\]
can be expressed as
\begin{equation*}
\eta_k=\frac{4}{h^2}\sin^2\frac{\alpha_k h}{2},\ k=\overline{1, N-3},
\end{equation*}
where $\alpha_k$ is root of any  of equations
\begin{align*}
&\sin  \frac{\alpha(1-\xi_2)}{2}=0, \quad
\sin  \frac{\alpha \xi_1}{2}=0, \quad
\sin  \frac{\alpha(1-\xi_1+\xi_2)}{2}=0.
 \end{align*}
\end{lem}
%

\begin{pf}
The proof is analogous to the proof of Lemma \ref{lem:1}.
Putting general solution of Equation \eqref{eq:5.5}
\begin{equation*}\label{eq:5}
v_i=c_1\cos( \alpha ih)+c_2\sin (\alpha ih),\quad i=\overline{1,N-1},\quad 0<\alpha<\frac{\pi}{h},
\end{equation*}
into the  Equations \eqref{eq:5.6}--\eqref{eq:5.7} we get a system of equations
\begin{align}
&\frac{h}{3}c_1 \sin (\alpha \xi_1) \frac{2+\cos( \alpha h) }{\sin( \alpha h)} +\frac{h}{3} c_2 (1- \cos \alpha \xi_1) \frac{2+\cos (\alpha h) }{\sin( \alpha h)} =0, \label{eq:5.9}\\
&\frac{h}{3}c_1 \big(\sin \alpha -\sin( \alpha \xi_2)\big) \frac{2+\cos (\alpha h) }{\sin ( \alpha h)} +\frac{h}{3} c_2 \big(1- \cos (\alpha \xi_1)\big) \frac{2+\cos( \alpha h) }{\sin( \alpha h)} =0. \nonumber
 \end{align}
 The system \eqref{eq:5.9} is analogous to \eqref{eq:3.12}, only instead of a multiplier $\frac{1}{2}\tan \frac{\alpha h}{2}\neq 0 $  there is another multiplier $(2+\cos \big(\alpha h)\big)/\sin( \alpha h) \neq 0$.
The condition  $D=0$ in both cases means that
\[
\sin \frac{\alpha(1-\xi_1+\xi_2)}{2}\sin\frac{\alpha \xi_1}{2}\sin \frac{\alpha(1-\xi_1)}{2}=0.
\]
So, the statement of lemma follows from here.  
\end{pf}

\begin{cor}\label{c:2}
The eigenvalues of the problem  \eqref{eq:5.5}--\eqref{eq:5.7} satisfying the  inequalities $0<\eta_k<4/{h^2}$  are the same as for the problem \eqref{eq:3.1}--\eqref{eq:3.2} when the trapezoidal formula is used to approximate the integral condition.
\end{cor}
\begin{lem}\label{lem:6}
In the case of even  numbers $n$, $s$ and $N$
\begin{equation*}
\eta=\frac{4}{h^2}\cosh^2\frac{\alpha h}{2},\
\end{equation*}
where $\alpha$ is root of equation
\begin{equation*}
\cosh (\alpha h)=2,
\end{equation*}
is the eigenvalue of the problem \eqref{eq:5.5}--\eqref{eq:5.7} of multiplicity  2  %is twice multiple eigenvalue
with  two  linearly independent corresponding eigenvectors $\mathbf{v}^1$ and $\mathbf{v}^2$:
\begin{align*}
&v_i^1=(-1)^i  \cosh(\alpha ih), \quad  v_i^2=(-1)^i  \sinh(\alpha ih),\quad  i=\overline{1,N-1},
 \end{align*}
 where $\alpha$ is connected with $\eta$ by formula
 \begin{equation*}
\cosh (\alpha h)=\frac{\eta h^2}{2}-1.
\end{equation*}
\end{lem}
\begin{pf}
In the case $\eta>0$, the general solution of the Equation \eqref{eq:5.5} is
\begin{equation*}
v_i=(-1)^i\big(c_1 \cosh(\alpha ih) +c_2 \sinh(\alpha ih)\big),\quad i=\overline{1,N-1}.
\end{equation*}
After substitution this expression into conditions \eqref{eq:5.6} and \eqref{eq:5.7}, we get the system of equations
\begin{align*}
&\frac{h}{3}c_1 \sinh( \alpha \xi_1) \frac{-2+\cosh( \alpha h) }{\sinh( \alpha h)} +\frac{h}{3} c_2 \big(-1+ \cosh( \alpha \xi_1)\big) \frac{-2+\cosh( \alpha h) }{\sinh (\alpha h)} =0, \\
&\frac{h}{3}c_1 \big(\sinh \alpha -\sinh( \alpha \xi_2)\big) \frac{-2+\cosh (\alpha h) }{\sinh (\alpha h)} \nonumber\\
& \quad +\frac{h}{3} c_2 \big(\cosh \alpha-\cosh(\alpha \xi_2)\big) \frac{-2+\cosh( \alpha h) }{\sin( \alpha h)} =0. \nonumber
 \end{align*}
It follows that
\begin{align*}
D&=\frac{h^2}{9}\left(\frac{-2+\cosh(\alpha h) }{\sinh( \alpha h)}\right)^2\left|                                                                        \begin{array}{ll}
\sinh( \alpha \xi_1)  & -1+ \cosh (\alpha \xi_1) \\
\sinh \alpha -\sinh( \alpha \xi_2) & \cosh \alpha-\cosh( \alpha \xi_2) \\
\end{array}                                                                      \right|   \\
&=\frac{h^2}{9}\left(\frac{-2+\cosh( \alpha h) }{\sinh( \alpha h)}\right)^2 4 \sinh  \frac{\alpha(1-\xi_1+\xi_2)}{2}  \sinh  \frac{\alpha(\xi_2-1)}{2} \sinh \frac{\alpha \xi_1}{2}=0.
 \end{align*}
 So, $D=0$ if
 $\cosh (\alpha h)=2$
 or
% \begin{equation*}
$\alpha=\frac{\mathrm{arccosh } \ 2}{h}$.
%\end{equation*}
Then it follows that
\begin{equation}\label{eq:5.10}
\eta=\frac{4}{h^2}\sinh ^2 \frac{\alpha h}{2}.
\end{equation}
Now, we can conclude that there are two linearly 
	independent vectors $\mathbf{v}^1$, $\mathbf{v}^2$:
\begin{align*}
&v_i^1=(-1)^i  \cosh(\alpha ih),\quad  v_i^2=(-1)^i  \sinh(\alpha ih),\quad  i=\overline{1,N-1},
 \end{align*}
corresponding to eigenvalue $\eta$, defined by \eqref{eq:5.10}.
\end{pf}
\begin{cor}\label{c:3}
All $N-1$ eigenvalues of the difference eigenvalue problem \eqref{eq:5.5}--\eqref{eq:5.7}
\[\eta_k=\frac{4}{h^2}\sin^2\frac{\alpha_k h}{2},\ k=\overline{1, N-3}
\]
from Lemma \ref{lem:5} and
\[
\eta_{N-2}=\eta_{N-1}=\frac{4}{h^2}\cosh^2\frac{\alpha h}{2}
\]
from Lemma \ref{lem:6} are positive.
\end{cor}

Now, we can formulate the alternating direction method for the system \eqref{eq:5.1}--\eqref{eq:5.4}. We rewrite this system in the matrix form
\begin{equation}\label{eq:5.11}
(\mathbf{A}_1+\mathbf{A}_2-\frac{h^2}{6}\mathbf{A}_1\mathbf{A}_2)\mathbf{U}=\mathbf{\Phi},
\end{equation}
where $\mathbf{A}_1=\mathbf{I}\otimes \mathbf\Lambda_x$, $\mathbf{A}_2=\mathbf\Lambda_y\otimes \mathbf{I}$ and
$\mathbf{\Phi}=\mathbf{\Phi}({\mathbf{F}},\mu^a,\mu^b,\mu^c,\mu^d,h)$.
  The newly defined matrices $\mathbf{A}_1$ and $\mathbf{A}_2$ have all the same properties as the matrices  defined by system \eqref{eq:2.1}--\eqref{eq:2.3}. In other words, Lemmas \ref{lem:3} and \ref{lem:4} are true for the matrices of the system of Equations \eqref{eq:5.1}--\eqref{eq:5.4}.

Since the matrices $\mathbf{A}_1$ and $\mathbf{A}_2$ commute, the system \eqref{eq:5.11} can be written in another form
\begin{equation}\label{eq:5.12}
(\mathbf{I}-\kappa \mathbf{A}_2)\mathbf{A}_1 \mathbf{U}+(\mathbf{I}-\kappa \mathbf{A}_1)\mathbf{A}_2 \mathbf{U}=\mathbf{\Phi},
\end{equation}
where $\kappa=h^2/12$. Further, we rewrite system \eqref{eq:5.12} in equivalent form
\begin{equation}\label{eq:5.13}
 \Bar{\mathbf{A}}_1 \mathbf{U}+ \Bar{\mathbf{A}}_2\mathbf{U}=\Bar{\mathbf{\Phi}},
\end{equation}
where
\begin{align*}
& \Bar{\mathbf{A}}_k=\mathbf{A}_k(\mathbf{I}-\kappa \mathbf{A}_k)^{-1}, \quad k=1, 2,  \\
&\Bar{\mathbf{\Phi}}=(\mathbf{I}-\kappa \mathbf{A}_1)^{-1}(\mathbf{I}-\kappa \mathbf{A}_2)^{-1}\mathbf{\Phi},
 \end{align*}
(see \cite{2009LMJ49n3o} for details).

Now, we write Peaceman-Rachford alternating direction method for the system \eqref{eq:5.13}:
\begin{align}
  &\frac{\mathbf{U}^{n+1/2}-\mathbf{U}^{n}}{\tau^1_{n}}+ \Bar{\mathbf{A}}_1 \mathbf{U}^{n+1/2}+
   \Bar{\mathbf{A}}_2 \mathbf{U}^{n}=\Bar{\mathbf{\Phi}},\nonumber\\
    &\frac{\mathbf{U}^{n+1}-\mathbf{U}^{n+1/2}}{\tau^2_{n}}
  +\Bar{\mathbf{A}}_1\mathbf{U}^{n+1/2}+\Bar{\mathbf{A}}_2\mathbf{U}^{n+1}=\Bar{\mathbf{\Phi}}.  \label{eq:5.14}
\end{align}

We will notice that if all eigenvalues of the matrices $\mathbf{A}_1$ and $\mathbf{A}_2$  are positive, then all eigenvalues of the matrices
$\Bar{\mathbf{A}}_1$ and $\Bar{\mathbf{A}}_2$ are also positive for sufficiently small values of $h$. Thus, it follows from Lemmas \ref{lem:5} and
\ref{lem:6}  that for  method \eqref{eq:5.8} the statements proven in Theorem \ref{thm:1} and Remark \ref{rem:2} are valid.





\section{Numerical experiments}\label{s:6}

In the section, some numerical examples will be computed to verify the numerical accuracy and
efficiency of difference schemes that we presented in the work. Truncation error analysis provides a widely applicable framework for analyzing the accuracy
of finite difference schemes. We consider three illustrative examples. The first example deals with the second order difference scheme. The second and third examples
demonstrate  empirical verification of the truncation error for higher order difference scheme   with or without a known exact solution, respectively. 

We consider a model problem
\begin{align*}
	& \dfrac{\partial^2 u}{\partial x^2}+ \dfrac{\partial^2 u}{\partial y^2}=e^{\pi x}\sin(\pi y)+\frac{x^2y^2}{2}, \quad (x, y)\in \Omega=\{0<x<1,\quad 0<y<1\},\\
	& \int_0^{\xi_1} u(x, y)dx =\frac{1}{6\pi}(\xi_1^3y^2\pi+6\sin(\pi y)e^{\pi \xi_1}-6\sin(\pi y)),\\
	& \int_{\xi_2}^1 u(x, y)dx =- \frac{1}{6\pi}(\xi_2^3y^2\pi-y^2\pi+6\sin(\pi y)e^{\pi \xi_2}-6\sin(\pi y)e^{\pi}), \\
	& u(x, 0) = 0, \quad u(x, 1) = {x^2}/{2},
\end{align*}
where $0<\xi_1<1$, \ $0<\xi_2<1$.
The exact solution is
\[
  u=e^{\pi x}\sin(\pi y)+{x^2y^2}/{2}.
\]
 We consider the uniform grids with the different mesh sizes $h = 1/2^{k}$, $\quad k=3, \dots, 7$ and analyze the convergence and the accuracy of the computed solution of the second and the fourth order difference schemes.
 Test problems were solved with the different values of the parameters
$\xi_1$, $\xi_2$.

We compute  the maximum norm of the error of the numerical solution with  respect to the exact
solution, which is defined as
%
\begin{equation*}
  \varepsilon_{h}=\max_{j=1,\cdots,n}\max_{i=1,\cdots,n}\abs{u(x_i, y_j)-u_{ij}}.
\end{equation*}
We define a number $p$ as
$p={\varepsilon_{2h}}/{\varepsilon_{h}}$,
which theoretically must be approximately $p\thickapprox 4$ for the second order difference scheme and $p\approx 16$ for the fourth order difference scheme.
\begin{exmp}\label{ex:1}
 The computational results for ADI method \eqref{eq:2.14}  are reported in Table \ref{tab:1}. We can clearly observe a second-order convergence in the maximum norm for all presented choices of $\xi_1$ and $\xi_2$. 
 From the last two columns of Table \ref{tab:1} it is clear that the number of iterations of the ADI method is quite accurately proportional to the value of $\log (1/h))$. When we were planning a numerical experiment, this is what we wanted to check (see the last paragraph of Section \ref{s:4}). Please note that $\xi_1$ and  $\xi_2$ are chosen  so that $\xi_1 = 1-\xi_2$, that is, the difference problem has multiple eigenvalues (see Corollary  \ref{cor:1}).
 \end{exmp}
\begin{exmp}\label{ex:2}
The outputs for the different values of the pairs of  parameters
$(\xi_1, \xi_2)$, respectively, together with the experimental convergence order for higher-order difference scheme \eqref{eq:5.1}--\eqref{eq:5.4}  are shown in
Table~\ref{tab:2}. We demonstrate a fourth-order convergence in the maximum norm for all choices of $\xi_1$ and $\xi_2$.

In the next example, we check that the empirically
observed convergence rates in experiments coincide with the theoretical value of it found
from truncation error analysis.
\end{exmp}
\begin{table}[H]
\beforecaption={\tabcolsep=0pt}
\caption{ Accuracy of the solution and the number of the iterations for the ADI method \eqref{eq:2.14}. }\label{tab:1} \centering
\raggedright\tabcolsep=5pt
\begin{tabular}{lllllcl}
 \hline
$\xi_1$ & $\xi_2$ & $h$ &  $\varepsilon_{h}$ &$p$ &  number of iter.& $\log(1/h)$  \\
\hline
  0.25    & 0.75  & $1/2^3$ & $ 0.3809$&     &  17& 2.079\\
   &  & $1/2^4$ & $ 0.09801$&  3.8863   & 22& 2.773 \\
  Case  &  &$1/2^5$ & $ 0.02468$& 3.9712  & 27 & 3.466\\
     $0<\xi_1<\xi_2<1$ &  & $1/2^6$ & $ 0.006182$& 3.9922  & 33 & 4.159\\
        &  & $1/2^7$ & $ 0.001546$&  3.9987 & 38 & 4.852\\
\hline
  0.5    & 0.5  & $1/2^3$     & $ 0.4385$&     &  15& 2.079\\
  &  & $1/2^4$  & $0.1128$&  3.8874   & 20& 2.773 \\
 Case   &  & $1/2^5$ & $ 0.02841$& 3.9704  & 25 & 3.466\\
  $0<\xi_1=\xi_2<1$    &  & $1/2^6$ & $ 0.007116$& 3.9924  & 31 & 4.159\\
        &  & $1/2^7$ & $ 0.001780$& 3.9978 & 36 & 4.852\\
 \hline
  0.75    & 0.25  & $1/2^3$ & $ 0.4484$&     &  13& 2.079\\
   &  & $1/2^4$ & $ 0.1154$&  3.8856   & 18& 2.773 \\
  Case  &  &$1/2^5$ & $ 0.02906 $& 3.9711  & 23 & 3.466\\
  $0<\xi_2<\xi_1<1$    &  & $1/2^6$ & $ 0.007278$& 3.9929  & 29 & 4.159\\
        &  & $1/2^7$ & $ 0.001820$& 3.9989 & 34 & 4.852\\
  \hline
\end{tabular}
\end{table}



\begin{table}[H]
\beforecaption={\tabcolsep=0pt}
\caption{ Accuracy of the solution and the number of the iterations for the  ADI method  \eqref{eq:5.14}. }\label{tab:2} \raggedright\tabcolsep=5pt
\begin{tabular}{lllllcl}
\hline
$\xi_1$ & $\xi_2$ & $h$ &  $\varepsilon_{h}$ &$p$ &  number of iter.& $\log(1/h)$ \\
\hline
 0.25    & 0.75  & 0.125     & $ 3.006\cdot 10^{-3}$&     &  24& 2.079\\
   &  & $1/2^4$ & $ 1.903\cdot 10^{-4}$&  15.7961    & 31& 2.773 \\
 Case   &  &$1/2^5$ & $ 1.193 \cdot 10^{-5}$& 15.9514  & 37 & 3.466 \\
    $0<\xi_1<\xi_2<1$  &  & $1/2^6$ & $ 7.463\cdot 10^{-7}$& 15.9855  & 44 & 4.159\\
        &  & $1/2^7$ & $ 4.638\cdot 10^{-8}$&  16.0910 & 51 & 4.852\\
  \hline

  0.5    & 0.5  & $1/2^3$ & $ 3.009\cdot 10^{-3}$&     &  21& 2.079\\
   &  & $1/2^4$ & $ 1.903\cdot 10^{-5}$&  15.8119    & 28& 2.773 \\
  Case  &  &$1/2^5$ & $ 1.193 \cdot 10^{-5}$& 15.9514  & 35 &3.466 \\
    $0<\xi_1=\xi_2<1$  &  & $1/2^6$ & $ 7.461\cdot 10^{-7}$& 15.9898  & 41 & 4.159\\
        &  & $1/2^7$ & $ 4.738\cdot 10^{-8}$& 15.7472 & 48 & 4.852\\
 \hline
 0.75    & 0.25  & $1/2^3$     & $ 3.010\cdot 10^{-3}$&     &  19& 2.079\\
   &  & $1/2^4$   & $ 1.903\cdot 10^{-5}$&  15.8171    & 26& 2.773 \\
 Case   &  &$1/2^5$   & $ 1.193 \cdot 10^{-5}$& 15.9514  & 32& 3.466\\
 $0<\xi_2<\xi_1<1$    &  & $1/2^6$ & $ 7.462\cdot 10^{-7}$& 15.9877  & 39 & 4.159\\
        &  & $1/2^7$ & $ 4.674\cdot 10^{-8}$& 15.9649 & 46 & 4.852\\
 \hline

\end{tabular}
\end{table}
  

\begin{exmp}\label{ex:3}
We solve this problem as test problem without a known solution. We use Runge’s rule for a practical error estimate for  higher-order method.
 We performed numerical experiments with the scheme and compared results with the ones demonstrated in Table \ref{tab:2}.
 As expected, there is  a fourth-order convergence in the maximum norm for all choices of $\xi_1$ and $\xi_2$. The results are recorded in Table \ref{tab:3}.

\begin{table}[H]
\beforecaption={\tabcolsep=0pt}
\caption{ Accuracy of the solution and the number of the iterations for the ADI method \eqref{eq:5.14} using Runge rule for error estimate.}\label{tab:3}
\raggedright\tabcolsep=5pt
\begin{tabular}{lllllcl}
\hline
$\xi_1$ & $\xi_2$ & $h$ &  $\varepsilon_{h}$ &$p$ &  number of iter.& $\log(1/h)$ \\
\hline
  0.25 & 0.75 & $1/2^4$ & $  1.87737299\cdot 10^{-4}$&    & 31& 2.773 \\
    &  &$1/2^5$ & $ 1.18900520\cdot 10^{-5}$& 15.7894  & 37 & 3.466 \\
   Case   &  & $1/2^6$ & $ 7.45593983\cdot 10^{-7}$& 15.9471  & 44 & 4.159\\
    $0<\xi_1<\xi_2<1$    &  & $1/2^7$ & $ 4.66591776\cdot 10^{-8}$&  15.9796 & 51 & 4.852\\
  \hline
 0.5  & 0.5 & $1/2^4$ & $ 1.87922737\cdot 10^{-4}$&      & 28& 2.773 \\
    &  &$1/2^5$ & $ 1.18929493 \cdot 10^{-5}$& 15.8012  & 35 &3.466 \\
 Case     &  & $1/2^6$ & $  7.45647480\cdot 10^{-7}$& 15.9498  & 41 & 4.159\\
  $0<\xi_1=\xi_2<1$      &  & $1/2^7$ & $ 4.66055705\cdot 10^{-8}$& 15.9991 & 48 & 4.852\\
 \hline
 0.75  & 0.25 & $1/2^4$   & $ 1.87955101\cdot 10^{-4}$&      & 26& 2.773 \\
    &  &$1/2^5$   & $ 1.18934557 \cdot 10^{-5}$&  15.8032 & 32& 3.466\\
  Case    &  & $1/2^6$ & $ 7.45651475\cdot 10^{-7}$& 15.9504  & 39 & 4.159\\
   $0<\xi_2<\xi_1<1$     &  & $1/2^7$ & $  4.66640055\cdot 10^{-8}$&  15.9792 & 46 & 4.852\\
 \hline

\end{tabular}
\end{table}
\end{exmp}

 
\section{Conclusions}\label{s:7}

We have proven that  for the second order and the fourth-order difference schemes, the spectrum of difference problem consists only from positive eigenvalues.  The convergence of the ADI method is proven with an additional assumption that the system of eigenvectors of difference problem is complete. A numerical experiment  showed that the ADI method with optimal parameters practically converges in all studied cases. In addition, it converges as quickly as in the case with Dirichlet conditions.

\subsection*{Acknowledgements}

The authors are grateful to Professor Raimondas Čiegis for the  valuable discussion  that improved the manuscript. 

%spell_to this point %********** End of text entry *****************
\bibliography{x}
%\emptypage
\end{document}
