\documentclass[reqno]{amsart} \usepackage{hyperref} \usepackage{graphicx} \AtBeginDocument{{\noindent\small \emph{Electronic Journal of Differential Equations}, Vol. 2016 (2016), No. 118, pp. 1--14.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu} \thanks{\copyright 2016 Texas State University.} \vspace{8mm}} \begin{document} \title[\hfilneg EJDE-2016/118\hfil Convergence of iterative methods for elliptic equations] {Convergence of iterative methods for elliptic equations with integral boundary conditions} \author[M. Sapagovas, O. \v{S}tikonien\.e, R. \v{C}iupaila, \v{Z}. Jok\v sien\.e \hfil EJDE-2016/118\hfilneg] {Mifodijus Sapagovas, Olga \v{S}tikonien\.e,\\ Regimantas \v{C}iupaila, \v{Z}ivil\.e Jok\v{s}ien\.e} \address{Mifodijus Sapagovas \newline Institute of Mathematics and Informatics, Vilnius University, Akademijos 4, LT-08663 Vilnius, Lithuania} \email{mifodijus.sapagovas@mii.vu.lt} \address{Olga \v{S}tikonien\.e \newline Faculty of Mathematics and Informatics, Vilnius University, Naugarduko 24, LT-03225 Vilnius, Lithuania} \email{olga.stikoniene@mif.vu.lt} \address{Regimantas \v{C}iupaila \newline Vilnius Gediminas Technical University, Sauletekio ave. 11, LT-10223 Vilnius, Lithuania} \email{regimantas.ciupaila@vgtu.lt} \address{\v{Z}ivil\.e Jok\v{s}ien\.e \newline Vytautas Magnus University, Vileikos str. 8, LT-44404 Kaunas, Lithuania.\newline Lithuanian University of Health Sciences, Eiveni\c{u} str. 4, LT-50009 Kaunas, Lithuania} \email{zivile.joksiene@fc.lsmuni.lt} \thanks{Submitted March 14, 2016. Published May 11, 2016.} \subjclass[2010]{35J25, 65N06, 65N22, 35P15} \keywords{Elliptic equation; finite-difference method; iterative method; \hfill\break\indent integral boundary conditions; eigenvalue problem; region of convergence} \begin{abstract} In this article, we consider the convergence of iterative method for the system of difference equations, approximating the elliptic two-dimensional equations with variable coefficients and integral boundary conditions. We investigate how convergence of iterative method depends on the structure of spectrum for difference operator with nonlocal conditions. The main goal of the paper is to analyze the influence of the monotonicity of the coefficient in the differential equation to extension (or reduction) of the region of convergence. \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{theorem}{Theorem}[section] \newtheorem{lemma}[theorem]{Lemma} \newtheorem{corollary}[theorem]{Corollary} \allowdisplaybreaks \section{Introduction}\label{s:1} Various phenomena of modern natural science can be described most conveniently in terms of differential equations with nonlocal boundary conditions. The theoretical investigation of nonlocal boundary value problems as well as numerical methods has been an important research area in various branches of mathematics. Some examples and details of application of such models can be found in monographs \cite{re41, re43, re45} and in many papers (see \cite{re43a, re27, re44, re3, re25, re42} and references therein). The first papers on the numerical methods for the elliptic equations with nonlocal conditions were published a few decades ago \cite{re1,re4,re3,re2}. Convergence of the finite difference method was one of the main issues considered. The newest results on various aspects of finite difference method could be found in \cite{re8, re6, re10, re9, re5, re7} (see also references given therein). Other numerical methods, different from finite difference methods, are presented in \cite{re13, re12, re14, re11}. Iterative methods for the system of difference equations approximating the elliptic equations with nonlocal conditions are insufficiently investigated. The first results having no posterior elaboration were obtained in papers \cite{re1, re15}. In \cite{re16, re19, re17, re18} the dependence of the convergence of the iterative methods on the structure of spectrum of the difference operators with nonlocal conditions was investigated. Presently, the eigenvalue problem for differential and difference operators with nonlocal conditions is one of the actively researched problems in the field of differential equations and numerical analysis. The structure of spectrum of such operators is relatively complicated even for simple eigenvalue problem \cite{re20}. Furthermore, this structure strongly depends on the change of parameters or functions in the nonlocal conditions \cite{re24, re21, re23, re22} (see also references given therein). More exhaustive list of references may be found in the review paper \cite{re25}. The eigenvalue problem for the differential and difference operator with the nonlocal conditions is related not only to the convergence of iterative methods, but also to the stability of the difference schemes for parabolic and hyperbolic equations \cite{re26, re28, re30, re29, re27}. In \cite{re33, re31, re32} the eigenvalue problem is investigated in connection with the existence, uniqueness and multiplicity of the solutions of the problems with nonlocal conditions. In the case of nonlocal conditions, the eigenvalue problems usually are considered for the differential problems with the constant coefficients. When the coefficients in the differential equation are variable, it becomes more difficult to investigate the structure of spectrum of the operator. Some specific results in the case of variable coefficients are obtained in \cite{re34, re36, re37, re11}. Some other methods of theoretical investigation or numerical analysis, different from spectral analysis, for various differential equations with variable coefficients and nonlocal conditions, are presented in \cite{re8, re6, re47, re13, re9, re12, re46, re14} etc. In the case of nonlocal boundary conditions the matrix of the system of difference equations is usually non-symmetric (one of the exceptions -- periodical boundary conditions). However, often it is possible to investigate when the difference operator has only positive eigenvalues. This allows to determine the region of convergence for iterative methods. Iterative methods for some linear and nonlinear elliptic equations with nonlocal conditions were investigated in \cite{re40, re39, re11}. Our objectives in this paper are to investigate some characteristic properties of the difference operator with nonlocal conditions (when the coefficients in the differential equation are variable) and to determine the convergence of iterative methods using these properties. Our goal is to investigate, how the monotonicity of the coefficient of the differential equation influences the expansion or reduction of the region of convergence. The paper is organized as follows. In Section \ref{s:2}, the differential and difference problems are formulated. The characteristic properties of the difference operator with nonlocal conditions are considered in Section \ref{s:3}. In Section \ref{s:4}, the main results are proposed. We determine, how the monotonicity of the coefficient of the differential equation shifts spectrum of the operator and influences the convergence of the iterative method. Numerical results are presented in Section \ref{s:5}. \section{Statement of the problem}\label{s:2} We consider two-dimensional elliptic equation with nonlocal integral conditions \begin{gather} \big(p(x)u_x\big)_x + \big(q(y)u_y\big)_y=f(x, y), \quad (x, y)\in D=\{0 0$, $q(y)> 0$. The problem is solved by using finite difference method. In the domain of the differential problem, we introduce the grids \begin{gather*} \overline{\omega}^{h}_{x}=\{x_0=0, \; x_1, \dots, x_N=1\}, \quad h=x_i-x_{i-1}=1/N, \\ \overline{\omega}^{h}_{y}=\{y_0=0, \; y_1, \dots, y_N=1\}, \quad h=y_j-y_{j-1}=1/N, \\ \omega^{h}_{x}=\{x_1,\dots,x_{N-1}\}, \quad \omega^{h}_{y}=\{y_1,\dots,y_{N-1}\}\\ \omega^{h}_{1/2, x}=\{x_{1/2}, x_{3/2},\dots,x_{N-1/2}\}, \quad \omega^{h}_{1/2, y}=\{y_{1/2}, y_{3/2},\dots,y_{N-1/2}\}. \end{gather*} In the closed domain $\overline{{D}}$ we use the grids $\overline{\omega}^h=\overline{\omega}^{h}_{x}\times\overline{\omega}^{h}_{y}$, ${\omega}^h={\omega}^{h}_{x}\times{\omega}^{h}_{y}$ If $\omega$ is one of these grids, then the space $\mathcal{F}(\omega)$ of grid functions can be defined on it. Let $U_{ij}$ approximate $u(x_i, y_j)$, where $x_i=ih$ and $y_j=jh$. We introduce the following grid operators: \begin{gather*} \partial_{\bar{x}}, \partial_{x}, \partial_{\bar{y}}, \partial_{y}: \mathcal{F}(\overline{\omega}^h)\to\mathcal{F}(\omega^h),\\ \partial_{\bar{x}} U_{ij}:=\frac{U_{ij}-U_{i-1,j}}{h_x}, \quad \partial_{x} U_{ij}:=\frac{U_{i+1,j}-U_{ij}}{h_x},\\ \partial_{\bar{y}} U_{ij}:=\frac{U_{ij}-U_{i,j-1}}{h_y}, \quad \partial_{y} U_{ij}:=\frac{U_{i,j+1}-U_{ij}}{h_y}. \end{gather*} The functions $p$, $q$, $f$, $\varphi_1$, $\varphi_2$, $\varphi_3$ and $\varphi_4$ of the differential problem are approximated by grid functions $p_i$ (on the grid $\omega^{h}_{1/2, x}$), $q_{j}$ (on the grid $\omega^{h}_{1/2, y}$), $f_{ij}$ ( on the grid ${\omega}^h$), $(\varphi_1)_{j}$, $(\varphi_2)_{j}$ (on the grid ${\omega}^h_y$), $(\varphi_3)_i$ and $(\varphi_4)_i$ (on the grid $\overline{\omega}^h_x$). The system of finite difference equations corresponding to the differential problem is as follows: \begin{gather} \partial_{x}( p_{i-1/2}\partial_{\bar{x}} U_{ij}) +\partial_{y}( q_{j-1/2}\partial_{\bar{y}} U_{ij}) =f_{ij}, \quad(x_i,y_j)\in\omega^h, \label{2.5}\\ U_{0j}=\gamma_1 h\Big( \frac{U_{0j}+ U_{Nj}}{2} +\sum_{i=1}^{N-1} U_{ij}\Big)+(\varphi_1)_j, \quad y_j\in {\omega}^h_y, \label{2.6}\\ U_{Nj}=\gamma_2 h\Big( \frac{U_{0j}+ U_{Nj}}{2} +\sum_{i=1}^{N-1} U_{ij}\Big)+(\varphi_2)_j, \quad y_j\in {\omega}^h_y, \label{2.7}\\ U_{i0}= (\varphi_3)_i, \quad U_{iN}= (\varphi_4)_i, \quad x_i\in \overline{\omega}^h_x. \label{2.8} \end{gather} This system of difference equations has $(N+1)(N-1)$ equations and unknowns $U_{ij}$, \ $i=\overline{0, N}, \ j=\overline{1, N-1}.$ The system \eqref{2.5}--\eqref{2.8} can be rewritten in the matrix form. We express the values $U_{0j}$ and $U_{Nj}$ from \eqref{2.6} and \eqref{2.7} via the remaining unknowns \begin{gather} U_{0j}=\frac{\gamma_1 h}{D}\sum_{i=1}^{N-1}U_{ij} +(\Tilde{\varphi}_1)_j, \quad j=\overline{1, N-1}, \label{2.10}\\ U_{Nj}=\frac{\gamma_2 h}{D} \sum_{i=1}^{N-1}U_{ij} +(\Tilde{\varphi}_2)_j, \quad j=\overline{1, N-1}, \label{2.11} \end{gather} where \begin{gather*} (\Tilde{\varphi}_1)_j=\frac{1}{D}\big((\varphi_1)_j+\frac{(\gamma_1(\varphi_2)_j -\gamma_2(\varphi_1)_j)h}{2}\big), \\ (\Tilde{\varphi}_2)_j=\frac{1}{D}\big((\varphi_2)_j+\frac{(\gamma_2(\varphi_1)_j -\gamma_1(\varphi_2)_j)h}{2}\big), \\ D=1-\frac{(\gamma_1+\gamma_2) h}{2}. \end{gather*} The condition $D\neq 0$ is necessary and sufficient for representing $U_{0j}$ and $U_{Nj}$ in the form \eqref{2.10}, \eqref{2.11}. In the case $D<0$ the structure of spectrum of difference operator with nonlocal condition can be qualitatively different from the structure of spectrum of differential one \cite{re48}. So we only consider $D>0$. If $\gamma_1+\gamma_2> 0$, the step of the grid should be small enough, \begin{equation*}\label{2.9} h<\frac{2}{\gamma_1+\gamma_2}. \end{equation*} There is no limitation on the step $h$ when $\gamma_1+\gamma_2\le 0$. Now we substitute the expressions of $U_{0j}$ and $U_{Nj}$ into difference equations \eqref{2.5} for the cases $i=1$ and $i=N-1$. So, restructuring \eqref{2.5}--\eqref{2.7} in such a way we obtain new system of equations, the order of which is $(N-1)^2$ and the number of the unknowns $U_{ij}$, $i,j=\overline{1, N-1}$ is also equal to $(N-1)^2$. We write this system in a matrix form \begin{equation}\label{2.12} \mathbf{A}\mathbf{U}=\mathbf{F}, \end{equation} where the order of square matrix $\mathbf{A}$ and vector $\mathbf{U}$ is $(N-1)^2$. The unknowns $U_{0j}$ and $U_{Nj}$ can be obtained by the formulas \eqref{2.10}, \eqref{2.11} after solving the system \eqref{2.12} . Reduction of \eqref{2.5}--\eqref{2.8} to the system \eqref{2.12} with lower order has corresponding eigenvalue problem \begin{gather} \partial_{x}( p_{i-1/2}\partial_{\bar{x}} U_{ij}) +\partial_{y}( q_{j-1/2}\partial_{\bar{y}} U_{ij}) +\lambda U_{ij} =0, \quad i,j=\overline{1, N-1}, \label{2.13} \\ U_{0j}=\gamma_1 h\Big( \frac{U_{0j}+ U_{Nj}}{2}+\sum_{i=1}^{N-1} U_{ij}\Big), \quad j=\overline{1, N-1}, \label{2.14}\\ U_{Nj}=\gamma_2 h\Big( \frac{U_{0j}+ U_{Nj}}{2}+\sum_{i=1}^{N-1} U_{ij}\Big), \quad j=\overline{1, N-1}, \label{2.15} \\ U_{i0}=0, \quad U_{iN}=0, \quad i=\overline{1, N-1}, \label{2.16} \end{gather} This eigenvalue problem with condition \eqref{2.9} is equivalent one to the algebraic eigenvalue problem with the same matrix $\mathbf{A}$ of order $(N-1)^2$ \begin{equation}\label{2.17} \mathbf{A}\mathbf{U}=\lambda \mathbf{U} \end{equation} (for details, see \cite{re17, re19, re28}). Exactly, if from \eqref{2.14} and \eqref{2.15} we would express $U_{0j}$ and $U_{Nj}$ (formulas \eqref{2.10} and \eqref{2.11} with $(\Tilde{\varphi}_1)_j=(\Tilde{\varphi}_2)_j=0$) and put these expressions to \eqref{2.13} we would get \eqref{2.17} with the same matrix as $\mathbf{A}$ in \eqref{2.12}. If we would consider \eqref{2.5}--\eqref{2.8} as the system with the matrix of order $(N+1)(N-1)$, then we would receive that the eigenvalues of this problem are not the eigenvalues of \eqref{2.13}--\eqref{2.16}. We note that when solving \eqref{2.5}--\eqref{2.8} by the iterative method there is no necessity to reduce it to the form as in \eqref{2.12} with the matrix of the lower order. The iterative methods could be presented as for \eqref{2.5}--\eqref{2.8}, as well for the system \eqref{2.12}. At the same time, when investigating the convergence of the iterative methods, we exploit spectrum of $\mathbf{A}$, i.e. the structure of spectrum of \eqref{2.13}--\eqref{2.16}. \section{Spectrum of matrix $\mathbf{A}$}\label{s:3} We shall investigate some important properties of spectrum of matrix $\mathbf{A}$, necessary for the convergence of iterative methods. With this aim we reduce the two-dimensional eigenvalue problem \eqref{2.13}--\eqref{2.16} to two separate one-dimensional eigenvalue problems. To apply the method of separation of variables we seek nontrivial separated solutions of \eqref{2.13} that satisfy the boundary conditions \eqref{2.14}--\eqref{2.16} and have the form \begin{equation}\label{3.1} U_{ij}=V_iW_{j}, \quad i,j=\overline{0, N}. \end{equation} Substituting such a product solution into \eqref{2.13}--\eqref{2.16} and separating the variables we obtain two one-dimensional eigenvalues problems: \begin{gather} \partial_{x}( p_{i-1/2}\partial_{\bar{x}} V_i)+\eta V_i =0, \quad i=\overline{1, N-1}, \label{3.2}\\ V_0=\gamma_1 \langle \mathbf{1},\mathbf{V}\rangle, \quad V_N=\gamma_2 \langle \mathbf{1},\mathbf{V}\rangle, \label{3.3} \end{gather} where we denote \begin{equation*} \langle \mathbf{1},\mathbf{V}\rangle:=h\Big( \frac{V_0+ V_N}{2}+\sum_{i=1}^{N-1} V_i\Big) \end{equation*} and \begin{gather} \partial_{y}( q_{j-1/2}\partial_{\bar{y}} W_{j})+\mu W_{j} =0, \quad j=\overline{1, N-1}, \label{3.4}\\ W_0=W_N=0. \label{3.5} \end{gather} The equality for eigenvalues is \begin{equation}\label{3.6} \lambda_{kl}=\eta_k+\mu_l, \quad k,l=\overline{1, N-1}. \end{equation} All the eigenvalues for \eqref{3.4}, \eqref{3.5} with the Dirichlet conditions are positive and for the smallest eigenvalue $\mu_1$ it is true that \begin{equation}\label{3.7} \min_l \mu_l=\mu_1\ge \min q(y)\cdot \frac{4}{h^2}\sin^2 \frac{\pi h}{2}. \end{equation} Let us consider some properties of the spectrum of difference eigenvalue problem \eqref{3.2}, \eqref{3.3} with nonlocal conditions. Finite difference scheme becomes unstable, if there exists at least one negative eigenvalue. Therefore, it is important to investigate the conditions for the appearance of negative eigenvalue. First of all we find when $\eta=0$ is the eigenvalue of this problem. We write the general solution of \eqref{3.2} in the form \begin{equation}\label{3.8} V_i=c_1{V}^{(1)}_i+c_2{V}^{(2)}_i, \end{equation} where ${V}^{(1)}_i$ and ${V}^{(2)}_i$ are two linear-independent solutions of \eqref{3.2} with $\eta=0$: \begin{align*} &\partial_{x}( p_{i-1/2}\partial_{\bar{x}} V^{(1)}_i) =0, \quad i=\overline{1, N-1}, \\ %\noalign{\smallskip} &V^{(1)}_0=0, \quad V^{(1)}_N=1 \end{align*} and \begin{gather*} \partial_{x}( p_{i-1/2}\partial_{\bar{x}} V^{(2)}_i) =0, \quad i=\overline{1, N-1}, \\ V^{(2)}_0=1, \quad V^{(2)}_N=0. \end{gather*} Let us denote \begin{equation}\label{3.9} F_i:=h \sum_{l=1}^{i} \frac{1}{p_{l-0.5}}, \quad i=\overline{1, N}, \quad F_0=0. \end{equation} Then we can write \begin{equation}\label{3.10} V^{(1)}_i= \frac{ F_i}{F_N}, \quad V^{(2)}_i= \frac{ F_N-F_i}{F_N}, \quad i=\overline{0, N}. \end{equation} \begin{lemma}\label{lem:1} The number $\eta=0$ is the eigenvalue of \eqref{3.2}--\eqref{3.3} if and only if \begin{equation}\label{3.11} \alpha \gamma_1+\beta \gamma_2=1, \end{equation} where \begin{equation}\label{3.12} \beta=\frac{h}{F_N}\Big(\frac{F_N}{2} +\sum_{i=1}^{N-1}F_i \Big), \quad \alpha=1-\beta. \end{equation} \end{lemma} \begin{proof} We require that \eqref{3.8} satisfies nonlocal conditions \eqref{3.3} and $V_i\not\equiv 0$. Substituting \eqref{3.8} into conditions \eqref{3.3} we obtain \begin{gather*} -\gamma_1 \langle \mathbf{1},\mathbf{V}^{(1)}\rangle c_1 -(1-\gamma_1\langle \mathbf{1},\mathbf{V}^{(2)}\rangle )c_2=0 \\ (1-\gamma_2 \langle \mathbf{1},\mathbf{V}^{(1)}\rangle) c_1 -\gamma_2\langle \mathbf{1},\mathbf{V}^{(2)}\rangle c_2=0. \end{gather*} For the condition $V_i\not\equiv 0$ to hold, it is necessary and sufficient for the system determinant to be equal zero: \[ \begin{vmatrix} -\gamma_1 \langle \mathbf{1},\mathbf{V}^{(1)}\rangle & 1-\gamma_1\langle \mathbf{1},\mathbf{V}^{(2)}\rangle \\ 1-\gamma_2 \langle \mathbf{1},\mathbf{V}^{(1)}\rangle &-\gamma_2\langle \mathbf{1},\mathbf{V}^{(2)}\rangle \end{vmatrix}= 0. \] This implies \begin{equation}\label{29a} \gamma_1\langle \mathbf{1},\mathbf{V}^{(2)}\rangle+\gamma_2 \langle \mathbf{1},\mathbf{V}^{(1)}\rangle -1=0. \end{equation} The desired equality \eqref{3.11} now follows immediately from the definition of $\langle \mathbf{1},\mathbf{V}^{(1)}\rangle$ and $\langle \mathbf{1},\mathbf{V}^{(2)}\rangle$. \end{proof} \begin{theorem}\label{thm:1} If $p(x)$ is increasing on the interval $(0, 1)$, then $0<\alpha<1/2$. If $p(x)$ is decreasing, then $1/2<\alpha<1$. \end{theorem} \begin{proof} Let us take $p'>0$, it means that \begin{equation}\label{3.13} \frac{ 1}{p_{1/2}}>\frac{ 1}{p_{3/2}}>\dots >\frac{ 1}{p_{i+1/2}}>\dots >\frac{ 1}{p_{N-1/2}}. \end{equation} According to these inequalities, the definition of $F_i$ \eqref{3.9} and the properties of the arithmetic mean, we obtain \[ \frac{F_i}{i}>\frac{F_N}{N}, \quad i=\overline{1, N-1}. \] or $F_i>h i F_N$. So, \[ \beta=\frac{h}{F_N}\big(\frac{F_N}{2} +\sum_{i=1}^{N-1}F_i \big)>h\big(\frac{1}{2}+ \sum_{i=1}^{N-1}h i \big) =\frac{1}{2}. \] Since $F_i< F_N$, \[ \beta< h\big(\frac{1}{2}+(N-1)\big)<1. \] So, $0<\alpha<1/2$. The second part of the theorem, when $p'<0$, is proved analogously, using the properties \begin{gather*} \frac{ 1}{p_{i-1/2}}<\frac{ 1}{p_{i+1/2}}, \\ \frac{F_i}{i}<\frac{F_N}{N}, \quad i=\overline{1, N-1}. \end{gather*} \end{proof} We call the line \eqref{3.11} a \emph{characteristic line} of problem \eqref{3.2}--\eqref{3.3}. It is remarkable that the difference operator has a zero eigenvalue on this characteristic line. In the system of coordinates $(\gamma_1, \gamma_2)$ \eqref{3.11} is an equation of line, which crosses the point $\gamma_1=1, \gamma_2=1$ (see line (1) in Figure \ref{fig1}). In the case of constant $p(x)=1$, $\alpha=\beta=1/2$, \eqref{3.11} simplifies to \begin{equation}\label{31} \gamma_1+\gamma_2=2, \end{equation} see dashed line (0) in Figure \ref{fig1}. \begin{figure}[htb] \begin{minipage}[t]{0.47\textwidth} \begin{center} \includegraphics[width=0.65\textwidth]{fig1} \end{center} \scriptsize{(a) increasing case: $p'>0$, $p(x)=1/(1-0.5x)$, $\gamma_1^0=2.25>2$, $\gamma_2^0=1.80<2$, $\gamma_1^*=3.48>\gamma_1^0$, $\gamma_2^*=2.73>\gamma_2^0$;} \end{minipage}\qquad \begin{minipage}[t]{0.47\textwidth} \begin{center} \includegraphics[width=0.65\textwidth]{fig2} \end{center} \scriptsize{(b) decreasing case: $p'<0$, $p(x)=1-0.5x$, $\gamma_1^0=1.80<2$, $\gamma_2^0=2.23>2$, $\gamma_1^*=3.46>\gamma_1^0$, $\gamma_2^*=4.33>\gamma_2^0$.} \end{minipage} \caption{Plots of the characteristic lines for different $p(x)$: the line (1) -- $\eta=0$ (1D case), the line (2) -- $\lambda=0$ (2D case), dashed line - $\eta=0$ (1D case, $p(x)=1$). }\label{fig1} \end{figure} From Theorem \ref{thm:1} we have the following result. \begin{corollary}\label{col:1} If $p(x)$ is increasing function, then \eqref{3.11} is placed between the inclined line \eqref{31} and horizontal line $\gamma_2=1$. If $p(x)$ is decreasing, then \eqref{3.11} is placed between the inclined line \eqref{31} and vertical line $\gamma_1=1$. \end{corollary} If the point $(\gamma_1, \gamma_2)$ is on \eqref{3.11}, then \eqref{3.2}, \eqref{3.3} with these values has the eigenvalue $\eta=0$. When $\gamma_1= \gamma_2=0$, all the eigenvalues are real and positive. So, all the real eigenvalues of \eqref{3.2}, \eqref{3.3} are positive, when in the plane $(\gamma_1, \gamma_2)$ the point in interest is placed below the line \eqref{3.11}. One negative eigenvalue exists, when the point is placed above \eqref{3.11}. \section{Convergence of iterative methods}\label{s:4} We will use the results of Section \ref{s:3} for the investigation of convergence of iterative methods for the system of difference equations \eqref{2.5}--\eqref{2.8}, provided in matrix form \eqref{2.12}. For this particular system we use Chebyshev iterative method \begin{equation}\label{32} \mathbf{U}^{k+1}=\mathbf{U}^{k}-\tau_{k+1}(\mathbf{A}\mathbf{U}^{k}-\mathbf{F}), \end{equation} where \[ \tau_{k}=\frac{2}{1+t_k}, \quad t_k=\cos\big(\frac{(2k-1)\pi}{2m}\big),\quad k=\overline{1,m}. \] When all eigenvalues of the matrix $\mathbf{A}$ are positive, the convergence of this method is investigated in \cite{re16, re38}. To simplify the problem, we take $q(y)=1$ in \eqref{2.1}. Then the least eigenvalue of \eqref{3.4}, \eqref{3.5} is $\mu_1=\frac{4}{h^2}\sin^2 \frac{\pi h}{2}$. For sufficiently small $h$, the value of $\mu_1$ can be written approximately $\mu_1\approx \pi^2$. Note that $\mu_1 < \pi^2$ and $\mu_1\approx \pi^2 + O(h^2)$. Particularly, $|\mu_1 -\pi^2|<0.01$ for $h<0.03$. According to \eqref{3.6}, $\lambda_{kl}=0$ if and only if $\eta_{k}=-\mu_1\approx -\pi^2$. Let us consider, when the eigenvalue problem \eqref{3.2}, \eqref{3.3} may have negative eigenvalue $\eta_{k}=\beta$, depending on the values of parameters $\gamma_1, \gamma_2$. Here $\beta<0$ is a fixed number. We use the general solution of the equation \eqref{3.2} with \begin{equation*} V_i(\beta)=c_1{V}^{(1)}_i(\beta)+c_2{V}^{(2)}_i(\beta), \end{equation*} where ${V}^{(1)}_i(\beta)$ and ${V}^{(2)}_i(\beta)$ are two linear-independent solutions of the equation \eqref{3.2} with $\eta=\beta<0$. Then similarly to \eqref{29a}, we obtain that $\beta<0$ is the eigenvalue of \eqref{3.2}, \eqref{3.3} if and only if the point $(\gamma_1, \gamma_2)$ is on the line \begin{equation}\label{33} \gamma_1\langle \mathbf{1},\mathbf{V}^{(2)}(\beta)\rangle +\gamma_2 \langle \mathbf{1},\mathbf{V}^{(1)}(\beta)\rangle=1. \end{equation} So, locus, where \eqref{3.2}, \eqref{3.3} has the eigenvalue $\eta=-\pi^2$ is the line. Intersection points $(\gamma_1^*, 0)$ and $(0, \gamma_2^*)$ (see Figure \ref{fig1}(b)) where this line crosses the coordinate axes are important. We also denote the points, where \eqref{3.11} crosses coordinate axes: $(\gamma_1^0, 0)$ and $(0, \gamma_2^0)$. Note that \begin{equation}\label{34} \gamma_1^*>\gamma_1^0, \quad \gamma_2^*>\gamma_2^0. \end{equation} \begin{corollary}\label{col:2} If one-dimensional eigenvalue problem \eqref{3.2}, \eqref{3.3} has no complex eigenvalues, then the region of the convergence of the iterative method \eqref{32} according to the parameters $\gamma_1, \gamma_2$ is determined by the following condition: the point $(\gamma_1, \gamma_2)$ on the coordinate plane must be placed below the line \eqref{33}, crossing the points $(\gamma_1^*, 0)$ and $(0, \gamma_2^*)$, where \begin{itemize} \item[(i)] $\gamma_1^*>2$, $\gamma_2^*>1$, if $p(x)$ is increasing function, \item[(ii)] $\gamma_1^*>1$, $\gamma_2^*>2$, if $p(x)$ is decreasing function. \end{itemize} \end{corollary} If the condition prescribed in Corollary \ref{col:2} is fulfilled, all eigenvalues of the matrix $\mathbf{A}$ are positive. Therefore iterative method converges. The equation of the line, crossing the points $(\gamma_1^*, 0)$ and $(0, \gamma_2^*)$, can be written as follows \begin{equation}\label{36} \frac{\gamma_1}{\gamma_1^*}+\frac{\gamma_2}{\gamma_2^*}=1. \end{equation} So the phrase "the point $(\gamma_1, \gamma_2)$ on the coordinate plane must be placed below the line \eqref{33}, crossing the points $(\gamma_1^*, 0)$ and $(0, \gamma_2^*)$" can be replaced by condition: the inequality \begin{equation}\label{37} \frac{\gamma_1}{\gamma_1^*}+\frac{\gamma_2}{\gamma_2^*}<1 \end{equation} is true. As it could be seen from the numerical results provided below, the values of the parameters $\gamma_1, \gamma_2$, depending on the concrete expression of the function $p(x)$, might be higher than it was indicated in Corollary \ref{col:2} as lowest limits. Note that as in the case $p(x)\equiv 1$, there is present some compensative mechanism -- if one of the parameters ($\gamma_1$ or $\gamma_2$) is ``wrong one'', the convergence might be ensured by another parameter, i.e. convergence depends not on the each of parameters $\gamma_1, \gamma_2$ separately, but on the generalized parameter $\alpha \gamma_1+\beta \gamma_2$. We emphasize the role of monotonicity of the function $p(x)$. So, generally, $\gamma_1^0>\gamma_2^0$ and $\gamma_1^*>\gamma_2^*$ if $p'>0$ and $\gamma_1^0<\gamma_2^0$, $\gamma_1^*<\gamma_2^*$ if $p'<0$, what was all the time observed in the numerical experiment with various functions $p(x)$. It is important to clarify the situation, when \eqref{3.2}, \eqref{3.3} has no complex eigenvalues. When the complex eigenvalues are not present, then the line \eqref{3.11}, on which the eigenvalue $\eta=0$ exists or the line \eqref{33}, on which $\lambda=0$, divides the coordinate plane into two parts. In one of these parts, to which the coordinate original point ($\gamma_1=0$, $\gamma_2=0$) belongs, all the eigenvalues are positive. In the remaining part one negative eigenvalue exists, all the others are positive. It follows from three statements. First, when $\gamma_1=0$, $\gamma_2=0$, then all eigenvalues are real and positive; second, $\eta=0$ is a simple non-multiple eigenvalue, and third, eigenvalues of the matrix are continuous function with respect to all matrix entries. All these three statements also remain true in the presence of complex eigenvalue. The line crossing the points $(\gamma_1^0, 0)$ and $(0, \gamma_2^0)$ still separates regions of convergence and non-convergence of iterative methods. Again, the region of convergence may shrink because of two reasons. First, some eigenvalues, for which $\mathrm{Re} \lambda_{kl}<0$ may arise. Second, as parameters $\gamma_1$ and $\gamma_2$ change, the positive eigenvalue may continuously become negative, passing not the value $\lambda_{kl}=0$, but the value, for which $\mathrm{Re} \lambda_{kl}=0$, $\mathrm{Im} \lambda_{kl}\neq 0$. Although, during the numerical experiment this situation was not observed, it was successfully modeled with another type of nonlocal conditions. \section{Numerical experiments}\label{s:5} Numerical experiments are performed to illustrate the theoretical results. We consider examples where the exact solutions of \eqref{2.1}--\eqref{2.4} are explicitly known by suitable choice of $f(x,y)$. In the first part of the numerical experiments we calculated the parameters $\gamma_1^0$ and $\gamma_2^0$, $\gamma_1^*$ and $\gamma_2^*$ characterizing the region of convergence of the iterative method. Recall that one-dimensional eigenvalue problem \eqref{3.2}, \eqref{3.3} with $\gamma_1^0$ has the eigenvalue $\eta=0$. The two-dimensional eigenvalue problem \eqref{2.13}--\eqref{2.16} with the value $\gamma_1^*$ (when $\gamma_2=0$) has the eigenvalue $\lambda=0$. If $\gamma_2=0$, $\gamma_1>\gamma_1^*$, the iterative method \eqref{32} diverges. We consider several choices of $p(x)$. \smallskip \noindent \textbf{Case 1.} \[ p(x)=\frac{1}{1-ax}, \quad p'(x)>0. \] In Table \ref{tab:1} the approximate values of $\gamma_1^0$ and $\gamma_1^*$, which are critical for convergence of iterative method, are presented for different values of parameter $a$. \renewcommand{\arraystretch}{1.3} \begin{table}[ht] \caption{Values of $\gamma_1^0$, $\gamma_1^*$ for increasing function $p(x)=1/(1-ax)$; $\gamma_2=0$, $\eta=0$.}\label{tab:1} \begin{center} \begin{tabular}{@ {\hspace{3pt}}lllllll}%{D{.}{.}{1.9}} \hline $a$ & 0 & 0.3 & 0.5 & 0.7 & 0.9 & 0.95 \\ \hline $\gamma_1^0$ & 2 & 2.13 & 2.25 & 2.44 & 2.75 & 2.86 \\ $\gamma_1^*$ & 3.42 & 3.44 & 3.48 & 3.57 & 3.76 & 3.84 \\ \hline \end{tabular} \end{center} \end{table} \smallskip \noindent\textbf{Case 2.} \[ p(x)=\frac{1}{x^2-ax+b}, \quad p'(x)>0. \] The numerical results are provided in Table \ref{tab:2}. In both of one-dimensional and two-dimensional eigenvalue problems the spectrum of the problem is much more sensitive to the change of parameters $a$ and $b$ comparing with Case 1. \begin{table}[ht] \caption{The values of $\gamma_1^0$, $\gamma_1^*$ for increasing function $p(x)=1/(x^2-ax+b)$; $\gamma_2=0$, $\eta=0$.}\label{tab:2} \begin{center} \begin{tabular}{@ {\hspace{3pt}}lllll}%{D{.}{.}{1.9}} \hline $a$ & 2.1 & 3 & 4 & 10 \\ $b$ & 1.1 & 2.05 & 3.05 & 9.05 \\ \hline $\gamma_1^0$ & 3.83 & 3.21 & 3.13 & 3.04 \\ $\gamma_1^*$ & 4.61 & 4.88 & 5.59 & 9.20 \\ \hline \end{tabular} \end{center} \end{table} \smallskip \noindent\textbf{Case 3. } $p(x)=1+bx$. Here the sign of $p'(x)$ depends on the sign of $b$. The numerical results presented in Table \ref{tab:3} show again that the statement in Corollary \ref{col:2} in a quantitative sense strongly depends on the function $p(x)$. \begin{table}[ht] \caption{The values of $\gamma_1^*$, $\gamma_2^*$ for function $p(x)=1+bx$.} \label{tab:3} \begin{center} \begin{tabular}{@ {\hspace{3pt}}lllll}%{D{.}{.}{1.9}} \hline $b$ & 5 & 0.5 & -0.5 & -0.95 \\ $p'(x)$ & $p'(x)>0$ & $p'(x)>0$ & $p'(x)<0$ & $p'(x)<0$ \\ \hline $\gamma_1^*$ & 3.59 & 3.44 & 3.46 & 3.42 \\ $\gamma_2^*$ & 1.95 & 2.96 & 4.33 & 10.78 \\ \hline \end{tabular} \end{center} \end{table} Note that $p(x)=1$, $q(y)=1$ imply $\gamma_1^*=\gamma_2^*\approx 3.42$. In all the cases of the numerical experiment (Tables \ref{tab:1} --\ref{tab:3}), we observed, that $\gamma_1^*> 3.42$ if $p'(x)>0$ and $\gamma_2^*> 3.42$ if $p'(x)<0$. However, it is not a theoretical statement, but practically reliable. The solution of the eigenvalue problem is influenced not only by the monotonicity of the function $p(x)$, but also by its absolute value. In Table \ref{tab:4} we provided the values of $\gamma_1^*$, $\gamma_2^*$, when \[ p(x)=\frac{c}{1+5x}, \] where $c$ varies. Function $p(x)$ is decreasing for $c>0$. However, the values $\gamma_1^0$, $\gamma_2^0$, i.e. the preconditions for the existence of zero eigenvalue for one-dimensional problem does not depend on $c$. This could be observed from the expression \eqref{3.12} for $\beta$. With any value of $c>0$, $\gamma_1^0\approx 1.615$, $\gamma_2^0\approx 2.625$. \begin{table}[ht] \caption{The values of $\gamma_1^*$, $\gamma_2^*$ for decreasing function $p(x)=c/(1+5x)$, $p'(x)<0$. }\label{tab:4} \begin{center} \begin{tabular}{@ {\hspace{3pt}}llllll}%{D{.}{.}{1.9}} \hline $c$ & 0.05 & 0.5 & 1 & 5& 20 \\ \hline $\gamma_1^*$ & 14.90 & 5.30 & 3.94 & 2.26& 1.79 \\ $\gamma_2^*$ & 31.41 & 10.57 & 7.50 & 3.88 & 2.96 \\ \hline \end{tabular} \end{center} \end{table} In the second part of numerical experiments the problem \eqref{2.5}--\eqref{2.8} was solved using the iterative method \eqref{33}. As it was mentioned before, the convergence depends on one generalized parameter \begin{equation}\label{36a} \Tilde{\gamma}= \frac{\gamma_1}{\gamma_1^*}+\frac{\gamma_2}{\gamma_2^*}. \end{equation} When $\Tilde{\gamma}>1$, the iterative method diverges due to existence of negative eigenvalue of $\mathbf{A}$. The role of the condition $\Tilde{\gamma}<1$ is quite obvious in the case of only one nonlocal condition, i.e. $\gamma_1=0$ or $\gamma_2=0$. This situation is typical for many practical problems \cite{re41, re42, re43}. \begin{table}[ht] \caption{The convergence of iterative method for different functions $p(x)$.} \label{tab:5} \begin{center} \begin{tabular}{@ {\hspace{3pt}}lllcc}%{D{.}{.}{1.9}} \hline $\gamma_1$ & $\gamma_2$ & $p(x)=1$ & $p(x)=1/(1-0.5x)$ & $p(x)=1-0.5x$ \\ & & & $p'(x)>0$ & $p'(x)<0$ \\ \hline 3 & 0 & conv. & conv. $(\gamma_1<\gamma_1^*\approx 3.48)$ & \\ % & & & $(\gamma_1<\gamma_1^*\approx 3.48)$ & \\ 0 & 3 & conv. & div. \ \ $(\gamma_2>\gamma_2^*\approx 2.73)$ & \\ %& & & $(\gamma_2>\gamma_2^*\approx 2.73)$ & \\ 0 & 4 & div. & & conv. $(\gamma_2<\gamma_2^*\approx 4.43)$ \\ %& & & & $(\gamma_2<\gamma_2^*\approx 4.43)$ & \\ 4 & 0 & div. & & div. \ \ $(\gamma_1>\gamma_1^*\approx 3.43)$ \\ %& & & & $(\gamma_1>\gamma_1^*\approx 3.43)$ & \\ \hline \end{tabular} \end{center} \end{table} \begin{table}[ht] \caption{The convergence of iterative method, depending of condition $\Tilde{\gamma}<1$, for different functions $p(x)$.}\label{tab:6} \begin{center} \begin{tabular}{@ {\hspace{3pt}}rrlcc}%{D{.}{.}{1.9}} \hline $\gamma_1$ & $\gamma_2$ & $p(x)=1$ & $p(x)=1/(1-0.5x)$ & $p(x)=1-0.5x$ \\ & & & $p'(x)>0$ & $p'(x)<0$ \\ \hline 2 & 2 & div. & div. \ \ \ \ $(\Tilde{\gamma}\approx 1.31>1)$ & div. $(\Tilde{\gamma}\approx 1.03>1)$ \\ -2 & -2 & conv. &conv. $(\Tilde{\gamma}\approx -1.31<1)$ & conv. $(\Tilde{\gamma}\approx -1.03<1)$ \\ \hline \end{tabular} \end{center} \end{table} In Table \ref{tab:5} the convergence of the iterative method is presented. These results fully correspond to the theoretical investigations (see also Figure \ref{fig1}). Table \ref{tab:6} is composed in analogous way. Tables \ref{tab:7} and \ref{tab:8} complement the results of numerical experiment, presented in Tables \ref{tab:5} and \ref{tab:6}. In these tables the errors of the solution \[ \varepsilon_h=\max_{i,j}|u_{ij}^n-u_{ij}^*| \] are provided, where $u_{ij}^n$ is the approximate solution of the system of difference equations, and $u_{ij}^*$ is the exact solution of the differential problem in the point $(x_i, y_j)$. We should admit, that all functions and coefficients in the differential equation \eqref{2.1} and boundary conditions \eqref{2.2}--\eqref{2.4} were choose so that $u(x, y)=1+exp(x+y)$ would be the solution of problem \eqref{2.1}--\eqref{2.4}. It is follows from Table \ref{tab:7}, that an error depends very little on the function $p(x)$ and it starts to grow [Table \ref{tab:8}], when the point $(\gamma_1, \gamma_2)$ comes closer to the line \eqref{36}. In this case the least positive eigenvalues tends to the zero. \begin{table}[ht] \caption{The values of error $\varepsilon_h$ for different functions $p(x)$.} \label{tab:7} \begin{center} \begin{tabular}{@ {\hspace{3pt}}llllll}%{D{.}{.}{1.9}} \hline $\gamma_1$ & $\gamma_2$ & $h$ & $p(x)=1$ & $p(x)=1/(1-0.5x)$ & $p(x)=1-0.5x$ \\ \hline 3 & 0 & $2^{-5}$ & 0.00762 & 0.00795 & 0.00733 \\ 3 & 0 & $2^{-6}$ & 0.00189 & 0.00197 & 0.00182 \\ 3 & 0 & $2^{-7}$ & 0.00047 & 0.00049 & 0.00045 \\ \hline \end{tabular} \end{center} \end{table} \begin{table}[ht] \caption{The values of error $\varepsilon_h$ for $\gamma_1=\gamma_2$, when $\Tilde{\gamma}\rightarrow 1$, $p(x)=1-0.5x$. }\label{tab:8} \begin{center} \begin{tabular}{@ {\hspace{3pt}}lrlllll}%{D{.}{.}{1.9}} \hline &$\gamma_1:$ & 1.8 & 1.9 & 1.93 & 1.94 \\ $h $& $\Tilde{\gamma}:$ &0.94 & 0.99 & 1.00 & 1.01 \\ \hline $2^{-5}$ & & 0.00795 & 0.0347 & 0.438 & div. \\ $2^{-6}$ & & 0.00197 & 0.0080 & 0.053 & div. \\ $2^{-7}$ & & 0.00049 & 0.0020 & 0.012 & div. \\ \hline \end{tabular} \end{center} \end{table} \subsection*{Acknowledgments} This research was partially supported by the Research Council of Lithuania (grant No. MIP-047/2014). \begin{thebibliography}{10} \bibitem{re8} A. Ashyralyev, E.Ozturk; \newblock On a difference scheme of second order of accuracy for the bitsadze-samarskii type nonlocal boundary-value problem. \newblock {\em Bound. Value Probl.}, 14(12):19 pp., 2014. \bibitem{re6} G. Avalishvili, M. Avalishvili, D. Gordeziani; \newblock On a nonlocal problem with integral boundary conditions for a multidimensional elliptic equation. \newblock {\em Appl. Math. Lett.}, 24:566--571, 2011. \bibitem{re34} B. I. Bandyrskii, V. L. Makarov; \newblock Sufficient conditions for the eigenvalue of the operator $-{\mathrm{d}}^2/{\mathrm{d}} x^2+q(x)$ with the {I}onkin--{S}amarskii conditions to be real valued. \newblock {\em Comput. Math. Math. Phys.}, 40(12):1715--1728, 2000. \bibitem{re10} G.~Berikelashvili, N.~Khomeriki; \newblock On the convergence rate of a difference solution of the poisson equation with fully nonlocal constraints. \newblock {\em Nonlinear Anal. Model. Control}, 19(3):367--381, 2014. \bibitem{re41} A. F. Chudnovskij; \newblock {\em Teplofizika {Pochv}}. \newblock Nauka, Moscow, 1976. \newblock (in Russian). \bibitem{re40} R.~\v{C}iupaila, M.~Sapagovas, O.~\v{S}tikonien{\.e}; \newblock Numerical solution of nonlinear elliptic equation with nonlocal condition. \newblock {\em Nonlinear Anal. Model. Control}, 18(4):412--426, 2013. \bibitem{re43a} M.~Dehghan; \newblock Efficient techniques for the second-order parabolic equation subject to nonlocal specifications. \newblock {\em Appl. Num. Math.}, 52(1):39--62, 2005. \bibitem{re24} A.~Elsaid, S.~Helal, A.~El-Sayed; \newblock The eigenvalue problem for elliptic partial differential equation with multi-point nonlocal conditions. \newblock {\em J. Appl. Anal. Comp.}, 5(1):146--158, 2015. \bibitem{re36} J.~Gao, D.~Sun, M.~Zhang; \newblock Structure of eigenvalues of multi-point boundary value problems. \newblock {\em Adv. Difference Equ.}, 2010, Article ID 381932, 24 pp., 2010. \bibitem{re1} D.~Gordeziani; \newblock {\em Solution Methods for a class of nonlocal boundary value problems}. \newblock Tbilisi, 1981. \newblock (in Russian). \bibitem{re26} A.V. Gulin; \newblock Stability of nonlocal difference schemes in a subspace. \newblock {\em Differ. Equ.}, 48(7):940--949, 2012. \bibitem{re33} J.~Henderson, S.K. Ntouyas; \newblock Positive solutions for systems of $n$th order three-point nonlocal boundary value problems. \newblock {\em Electron. J. Qual. Theory Differ. Equ.}, 2007(18):1--12, 2007. \bibitem{re4} V.A. Il'in, E.I. Moiseev; \newblock Two-dimensional nonlocal boundary value problem for {P}oisson's operator in differential and difference variants. \newblock {\em Mat. Model.}, 2:132--156, 1990. \newblock (in Russian). \bibitem{re31} G.~Infante; \newblock Eigenvalues of some nonlocal boundary-value problem. \newblock {\em Proc. Edinb. Math. Soc., II Ser.}, 46:75--86, 2003. \bibitem{re28} F.~Ivanauskas, T.~Me\v{s}kauskas, M.~Sapagovas; \newblock Stability of difference schemes for two-dimensional parabolic equations with nonlocal boundary conditions. \newblock {\em Appl. Math. Comput.}, 215(7):2716--2732, 2009. \bibitem{re30} F. F. Ivanauskas, Yu. A. Novitski, M. P. Sapagovas; \newblock On the stability of an explicit difference scheme for hyperbolic equations with nonlocal boundary conditions. \newblock {\em Differ. Equ.}, 49(7):849--856, 2013. \bibitem{re29} J.~Jachimavi\v(c)ien{\.e}, M.~Sapagovas, A.~\v{S}tikonas, O.~\v{S}tikonien{\.e}; \newblock On the stability of explicit finite difference schemes for a pseudoparabolic equation with nonlocal conditions. \newblock {\em Nonlinear Anal. Model. Control}, 19(2):225--240, 2014. \bibitem{re27} G.~Kalna, S.~McKee; \newblock The thermostat problem with a nonlocal nonlinear boundary conditions. \newblock {\em IMA J. Appl. Math.}, 69:437--462, 2004. \bibitem{re47} A. I. Kozhanov, L. S. Pul'kina; \newblock On the solvability of boundary value problems with a nonlocal boundary condition of integral form for multidimentional hyperbolic equations. \newblock {\em Differ. Equ.}, 42(9):1233--1246, 2006. \bibitem{re32} R.~Ma; \newblock A survey on nonlocal boundary value problems. \newblock {\em Appl. Math. E-Notes}, 7:257--279, 2007. \bibitem{re13} J.~Mart\'{\i}n-Vaquero; \newblock Polynomial-based mean weighted residuals methods for elliptic problems with nonlocal boundary conditions in the rectangle. \newblock {\em Nonlinear Anal. Model. Control}, 19(3):448--459, 2014. \bibitem{re44} J.~Martin-Vaquero, J.~Vigo-Aguiar; \newblock On the numerical solution of the heat conduction equations subject to nonlocal conditions. \newblock {\em Appl. Num. Math.}, 59(10):2507--2514, 2009. \bibitem{re43} A. M. Nakhushev; \newblock {\em The {E}quations of {M}athematical {B}iology}. \newblock Moscow, 1995. \newblock (in Russian). \bibitem{re9} C.~Nie, S.~Shu, H.~Yu, Q.~An; \newblock A high order composite scheme for the second order elliptic problem with nonlocal boundary and its fast algorithms. \newblock {\em Appl. Math. Comp.}, 227:212--221, 2014. \bibitem{re12} C.~Nie, H.~Yu; \newblock Some error estimates on the finite element approximation for two-dimensional elliptic problem with nonlocal boundary. \newblock {\em Appl. Numer. Math.}, 68:31--38, 2013. \bibitem{re21} S.~Pe{\v{c}}iulyt{\.{e}}, O.~\v{S}tikonien{\.{e}}, A.~\v{S}tikonas; \newblock Investigation of negative critical points of the characteristic function for problems with nonlocal boundary conditions. \newblock {\em Nonlinear Anal. Model. Control}, 13(4):467--490, 2008. \bibitem{re46} L. S. Pul'kina; \newblock Solution to nonlocal problems of pseudohyperbolic equations. \newblock {\em Electron. J. Differ. Equ.}, 2012(116):1--9, 2012. \bibitem{re14} S.~Sajavi\v(c)ius; \newblock Radial basis function method for a multidimensional linear elliptic equation with nonlocal boundary conditions. \newblock {\em Comp. Math. Appl.}, 67(7):1407--1420, 2014. \bibitem{re3} M.~Sapagovas; \newblock Numerical methods for two-dimensional problem with nonlocal conditions. \newblock {\em Differ. Equ.}, 20(7):1258--1266, 1984. \bibitem{re15} M.~Sapagovas; \newblock The solution of difference equations, arising from differential problem with integral condition. \newblock In G.I. Marchuk, editor, {\em Comput. Process and Systems}, number~6, pages 245--253, Moscow, Nauka, 1988. \newblock (In Russian). \bibitem{re16} M.~Sapagovas; \newblock The eigenvalues of some problem with a nonlocal condition. \newblock {\em Differ. Equ.}, 38(7):1020--1026, 2002. \bibitem{re5} M.~Sapagovas; \newblock Difference method of increased order of accuracy for the {P}oisson equation with nonlocal conditions. \newblock {\em Differ. Equ.}, 44(7):1018--1028, 2008. \bibitem{re2} M.~Sapagovas, R.~\v{C}iegis; \newblock On some boundary problems with nonlocal conditions. \newblock {\em Differ. Equ.}, 23(7):1268--1274, 1987. \bibitem{re37} M.~Sapagovas, R.~\v{C}iupaila, \v{Z}.~Jok{\v s}ien{\.e}; \newblock The eigenvalue problem for a one-dimensional differential operator with a variable coefficient and nonlocal integral conditions. \newblock {\em Lith. Math. J.}, 54(3):345--355, 2014. \bibitem{re23} M.~Sapagovas, R.~\v{C}iupaila, \v{Z}.~Jok{\v s}ien{\. e}, T.~Me\v{s}kauskas; \newblock Computational experiment for stability analysis of difference schemes with nonlocal conditions. \newblock {\em Informatica}, 24(2):275--290, 2013. \bibitem{re48} M.~Sapagovas, T.~Me{\v s}kauskas, F.~Ivanauskas; \newblock Numerical spectral analysis of a difference operator with non-local boundary conditions. \newblock {\em Appl. Math. Comput.}, 218(14):7515--7527, 2012. \bibitem{re20} M.~Sapagovas, A.~\v{S}tikonas; \newblock {On the structure of the spectrum of a differential operator with a nonlocal condition}. \newblock {\em Differ. Equ.}, 41(7):961--969, 2005. \bibitem{re19} M.~Sapagovas, A.~\v{S}tikonas, O.~\v{S}tikonien{\.e}; \newblock{Alternating Direction Method for the {P}oisson Equation with Variable Weight Coefficients in an Integral Condition}. \newblock{\em Differ. Equ.}, 47(8):1163--1174, 2011. \bibitem{re17} M.~Sapagovas, O.~\v{S}tikonien{\.e}; \newblock A fourth-order alternating-direction method for difference schemes with nonlocal condition. \newblock {\em Lith. Math. J.}, 49(3):309--317, 2009. \bibitem{re39} M.~Sapagovas, O.~\v{S}tikonien{\.e}; \newblock Alternating-direction method for a mildly nonlinear elliptic equation with nonlocal integral conditions. \newblock {\em Nonlinear Anal. Model. Control}, 16(2):220--230, 2011. \bibitem{re45} K.~Sch\"{u}gerl; \newblock {\em Bioreaction {E}ngineering: {R}eactions {I}nvolving {M}icroorganisms and Cells, vol.1}. \newblock Wiley, N.Y., 1987. \bibitem{re22} A.~\v{S}tikonas; \newblock Investigation of characteristic curve for {S}turm--{L}iouville problem with nonlocal boundary conditions on torus. \newblock {\em Math. Model. Anal.}, 16(1):1--22, 2011. \bibitem{re25} A.~\v{S}tikonas; \newblock A survey on stationary problems with nonlocal boundary conditions. \newblock {\em Nonlinear Anal. Model. Control}, 19(3):301--334, 2014. \bibitem{re18} O.~{\v{S}}tikonien{\.e}; \newblock Numerical investigation of fourth-order alternating direction method for {P}oisson equation with integral conditions. \newblock In V.~Kleiza, editor, {\em Proceeding of Intern. Conf. Differ. Equations and their Applic. DETA-2009}, pages 139--146, Kaunas, University of Technology, 2009. \bibitem{re38} O.~\v{S}tikonien{\.e}, M.~Sapagovas, R.~\v{C}iupaila; \newblock On iterative methods for some elliptic equations with nonlocal conditions. \newblock {\em Nonlinear Anal. Model. Control}, 19(3):517--535, 2014. \bibitem{re42} V. A. Vodachova; \newblock A boundary value problem with {N}akhushev nonlocal condition for a certain pseudoparabolic moisture-transfer equation. \newblock {\em Differ. Equ.}, 18:280--285, 1982. \newblock (in Russian). \bibitem{re7} E. A. Volkov, A. A. Dosiyev, S.~C. Buranay; \newblock On the solution of a nonlocal problem. \newblock {\em Comput. Math. Appl.}, 66(3):330--338, 2013. \bibitem{re11} Y.~Wang; \newblock Solutions to nonlinear elliptic equations with a nonlocal boundary condition. \newblock {\em Electron. J. Differ. Equ.}, pages 1--16, 2002. \end{thebibliography} \end{document}