An Exact Finite-Difference Solution of a Boundary-Value Problem for an Ordinary Linear Differential Equation of Second Order

B. F. Krutoi

(Presented by the scientific seminar of the Department of Engineering and Computational Mathematics)

An ordinary linear differential equation of second order is given,

\[ y''(x)+p(x)y'(x)+q(x)y(x)=f(x) \tag{1} \]

with boundary conditions

\[ \begin{cases} \alpha_0 y(a)+\alpha_1 y'(a)=A,\\ \beta_0 y(b)+\beta_1 y'(b)=B \end{cases} \tag{2} \]

at the points \(a, b\) of the \(Ox\) axis (Fig. 1). The coefficients \(p(x)\), \(q(x)\) and the free term \(f(x)\) of equation (1) will be assumed to be continuous functions of \(x\) on the interval \(a\le x\le b\). Under these conditions it is required to find the solution \(y(x)\) of the boundary-value problem (1), (2), taking into account that then the function \(y(x)\) and its derivatives \(y'(x)\), \(y''(x)\) will be continuous on the interval \(a\le x\le b\).

Fig. 1. Regions of application of expansions (28), (29), and (30) for computing \(y_s'\) and \(y_s''\).

We shall set as our goal to find a numerical solution of the indicated boundary-value problem by determining the values \(y_s=y(x_s)\) at \(n+1\) equally spaced ...

points \(x_s\) of the interval \([a,b]\), including its endpoints \(a=x_0\) and \(b=x_n\). Thus, if we denote by \(h\) the length of one subdivision of the interval \([a,b]\) when it is divided into \(n\) equal parts,

\[ h=\frac{b-a}{n}=\frac{x_n-x_0}{n}, \tag{3} \]

then the corresponding points \(x_s\) of the interval \([a,b]=[x_0,x_n]\) are defined as follows:

\[ x=x_0+sh,\qquad (s=0,1,2,\ldots,n). \tag{4} \]

Suppose now that for the \(n-1\) interior points \(x_s\), \((s=1,2,\ldots,n-1)\), of the interval \([x_0,x_n]\), specified according to (3)—(4), we have written the \(n-1\) equalities obtained by substituting \(x=x_s\) into the differential equation (1). Adding to these two more boundary conditions (2), we obtain a system of \(n+1\) algebraic equations of the first degree

\[ \begin{aligned} &1)\quad y''(x_s)+p(x_s)y'(x_s)+q(x_s)y(x_s)=f(x_s),\\ &\hspace{4.2cm}(s=1,2,\ldots,n-1)\\ &2)\quad \alpha_0y(x_0)+\alpha_1y'(x_0)=A,\\ &\hspace{1.1cm}\beta_0y(x_n)+\beta_1y'(x_n)=B \end{aligned} \tag{5} \]

with \(3(n+1)\) unknowns \(y_s=y(x_s)\), \(y'_s=y'(x_s)\), \(y''_s=y''(x_s)\) at the \(n+1\) points \(x_0,x_1,x_2,\ldots,x_{n-1},x_n\) of the interval \([x_0,x_n]\).

Since in the system of equations (5) the number \(3(n+1)\) of unknowns \(y_s\), \(y'_s\), and \(y''_s\), where \(s=0,1,2,\ldots,n-1,n\), is three times the number \(n+1\) of the equations themselves, from this system the indicated set of unknowns \(y_s\), \(y'_s\), \(y''_s\) cannot be found even approximately unless the method of least squares is used for this purpose.

However, there is a simple way to eliminate from the system (5) the \(2(n+1)\) unknowns \(y'_s\), \(y''_s\), if one takes a sufficiently large number \(n\) of subdivisions of the interval \([x_0,x_n]\) and takes into account that then the derivatives \(y'=y'(x_s)\) and \(y''_s=y''(x_s)\) can be expressed with a high degree of accuracy in terms of the values \(y_x=y(x_x)\), where \(x=0,1,2,\ldots,n-1,n\). We shall show how this is done.

First of all, note that if a point \(x\) of the axis \(Ox\) lies on the interval \([x_0,x_n]\), and hence \(x_0\le x\le x_n\), then this same point \(x\) also satisfies the narrower condition \(x_s\le x\le x_{s+1}=x_s+h\). The latter condition is equivalent to the relations

\[ 1)\ 0\le \xi\le h,\quad 2)\ x=x_s+\xi,\quad 3)\ x_{s+1}=x_s+h, \tag{6} \]

whence it follows that one may then write

\[ 1)\ y(x)=y(x_s+\xi),\quad 2)\ y(x_{s+1})=y(x_s+h). \tag{7} \]

We next introduce the finite-difference operators \(E,\Delta,\delta,\mu,\nabla\), defining them by the equalities [1]:

\[ \begin{aligned} &1)\quad y(x_{s+1})=y(x_s+h)=E^h y(x_s),\qquad y(x_{s-1})=y(x_s-h)=E^{-h}y(x_s),\\ &2)\quad y(x)=y(x_s+\xi)=E^\xi y(x_s),\\ &3)\quad y(x_{s+1})-y(x_s)=y(x_s+h)-y(x_s)=E^h y(x_s)-y(x_s)=(E^h-1)y(x_s),\\ &4)\quad y(x_s)-y(x_{s-1})=y(x_s)-y(x_s-h)=y(x_s)-E^{-h}y(x_s)=(1-E^{-h})y(x_s), \end{aligned} \]

$$ \begin{gathered} 5)\quad y(x_{s-1})-y(x_s)=(E^h-1)y(x_s)=\Delta y(x_s), \tag{8}\\ 6)\quad y(x_{s+1})-y(x_s)=(E^h-1)y(x_s)\delta y\!\left(x_{s+\frac12}\right) =\delta y\!\left(x_s+\frac{h}{2}\right)=\delta E^{\frac{h}{2}}y(x_s),\\ 7)\quad y(x_s)-y(x_{s-1})=(1-E^{-h})y(x_s) =\delta y\!\left(x_{s-\frac12}\right)=\delta y\!\left(x_s-\frac{h}{2}\right) =\delta E^{-\frac{h}{2}}y(x_s),\\ 8)\quad \frac12\left[y\!\left(x_{s+\frac12}\right)y\!\left(x_{s-\frac12}\right)\right] =\frac12\left[y\!\left(x_s+\frac{h}{2}\right)+y\!\left(x_s-\frac{h}{2}\right)\right]=\\ =\frac12\left(E^{\frac{h}{2}}+E^{-\frac{h}{2}}\right)y(x_s)=\mu y(x_s),\\ 9)\quad y(x_s)-y(x_{s-1})=(1-E^{-h})y(x_s)=\nabla y(x_s). \end{gathered} $$

From equalities (8) there follow the following relations among the indicated operators:

$$ \begin{gathered} 1)\quad 1+\Delta+E^h,\qquad 2)\quad (1-\nabla)=E^{-h},\qquad 3)\quad \delta E^{\frac{h}{2}}=(E^h-1)=\Delta,\\ 4)\quad \delta E^{-\frac{h}{2}}=(1-E^{-h})=\nabla,\qquad 5)\quad \delta=E^{\frac{h}{2}}-E^{-\frac{h}{2}},\\ 6)\quad (1+\Delta)(1-\nabla)=1,\qquad 7)\quad \mu=\frac12\left(E^{\frac{h}{2}}+E^{-\frac{h}{2}}\right). \end{gathered} $$

Using the relations found among the operators $E$, $\Delta$, $\nabla$, $\delta$, and $\mu$, we shall write the representations we need for $y(x)=y(x_s+\xi)$ in terms of descending $\Delta^\nu y_s$, central $\delta^\nu y_s$, and ascending $\nabla^\nu y_s$ finite differences, with accuracy up to terms of the 6th order of smallness. For this purpose we first reduce relation (6.1) to the form

$$ \frac{0}{h}=0\leq \frac{\xi}{h}=t\leq \frac{h}{h}=1 $$

and then take $t=\dfrac{\xi}{h}$ as the new variable. Thus, we obtain the following relations for $\xi$ and $t$:

$$ 1)\quad \xi=th,\qquad 2)\quad 0\leq \xi\leq h,\qquad 3)\quad 0\leq t\leq 1. \tag{10} $$

Now we can write, with accuracy up to terms of the 6th order of smallness, that

$$ \begin{aligned} 1)\quad y(x)&=y(x_s+\xi)=E^\xi y(x_s)=E^{th}y(x_s)=(E^h)^t y(x_s)=\\ &=(1+\Delta)^t y_s=\left[1+(t)_1\Delta+(t)_2\Delta^2+(t)_3\Delta^3+(t)_4\Delta^4+(t)_5\Delta^5+(t)_6\Delta^6\right]y_s \end{aligned} $$

and, consequently,

$$ y(x_s+\xi)=y_{s+t}=(1+\Delta)^t y_s \approx \sum_{\nu=0}^{6}(t)_\nu \Delta^\nu y_s, \tag{11} $$

where it is denoted

$$ (t)_\nu=\frac{t(t-1)\cdots[t-(\nu-1)]}{\nu!};\qquad (t)_0=1; \tag{12} $$

\[ \begin{aligned} 2)\quad y(x_s+\xi)=y_{s+t}=y_s'&+\frac{1}{2}t(\delta y_{s+1/2}+\delta y_{s-1/2})+\frac{t^2}{2!}\delta^2y_s+ \\ &+\frac{1}{2}\frac{t(t^2-1)}{3!}(\delta^3y_{s+1/2}+\delta^3y_{s-1/2}) +\frac{t^2(t^2-1)}{4!}\delta^4y_s+ \tag{13}\\ &+\frac{1}{2}\frac{t(t^2-1)(t^2-4)}{5!}(\delta^5y_{s+1/2}+\delta^5y_{s-1/2}) +\frac{t^2(t^2-1)(t^2-4)}{6!}\delta^6y_s=\\ &=y_s+t\mu\delta y_s+\frac{t^2}{2!}\delta^2y_s +\frac{t(t^2-1)}{3!}\mu\delta^3y_s +\frac{t^2(t^2-1)}{4!}\delta^4y_s+\\ &+\frac{t(t^2-1)(t^2-4)}{5!}\mu\delta^5y_s +\frac{t^2(t^2-1)(t^2-4)}{6!}\delta^6y_s; \end{aligned} \]

\[ \begin{aligned} 3)\quad y(x)&=y(x_s-\xi)=E^{-\xi}y(x_s)=(E^{-h})^t y(x_s)=\\ &=y_{s-t}=(1-\nabla)^t y_s=\sum_{\nu=0}^{6}(-1)^\nu (t)_\nu \nabla^\nu y_s, \end{aligned} \tag{14} \]

where the coefficients \((t)\) are determined according to (12).

Let us note that the representations (11) and (13) for \(y(x_s+\xi)=y_{s+t}\) and \(y(x_s-\xi)=y_{s-t}\) are called Newton expansions in descending \(\Delta^\nu y_s\) and ascending \(\nabla^\nu y_s\) finite differences. The derivation of these expansions is obvious.

As for representation (12) for \(y(x_s+\xi)=y_{s+t}\), it is called Stirling’s expansion in central differences \(\delta^\nu y_{s+x}\) and is obtained by a rather complicated route. To obtain this expansion, one must first construct Fraser’s diagram [2] for descending differences \(\Delta^\nu y_x\). Then, using this diagram, we write Gauss expansions along the upper and lower zigzags, after which we recalculate the differences \(\Delta^\nu y_x\) into the differences \(\delta^\nu y_x\), using relation (4.5). Then we obtain the following two Gauss expansions for \(y(x_s+\xi)=y_{s+t}\):

\[ 1)\quad y(x_s+\xi)=y_{s+t}=y_s+(t)_1\delta y_{s-1/2}+(t+1)_2\delta^2y_s+ \]

\[ +(t+1)_3\delta^3y_{s-1/2}+(t+2)_4\delta^4y_s+(t+2)_5\delta^5y_{s-1/2}+(t+3)_6\delta^6y_s; \tag{14} \]

\[ 2)\quad y(x_s+\xi)=y_{s+t}=y_s+(t)_1\delta y_{s+1/2}+(t)_2\delta^2y_s+ \]

\[ +(t+1)_3\delta^3y_{s+1/2}+(t-1)_4\delta^4y_s+(t+2)_5\delta^5y_{s+1/2}+(t+2)_6\delta^6y_s. \tag{15} \]

The half-sum \(\dfrac{1}{2}[(14)+(15)]\) of these Gauss expansions gives, after simple transformations, Stirling’s expansion (12).

Let us return again to the expansions (11)—(13). Our goal is to obtain from them series for computing \(y_s'\) and \(y_s''\), which enter the boundary-value problem (1)—(2), represented in finite-difference form in (5). To construct such series for \(y_s'\) and \(y_s''\), we shall proceed as follows.

The expression for the first derivative \(y_s'\) for all three expansions (11)—(13) will be obtained by the general rule, namely:

\[ 1)\quad y_s'=y'(x_s)=\lim_{\xi\to0}\frac{y(x_s+\xi)-y(x_s)}{\xi} =\lim_{t\to0}\frac{y(x_s+th)-y(x_s)}{th}= \]

\[ =\lim_{t\to 0}\frac{1}{t!}\left[t\Delta+\frac{t(t-1)}{2!}\Delta^2+\frac{t(t-1)(t-2)}{3!}\Delta^3+\cdots+\frac{t(t-1)\cdots(t-5)}{6!}\Delta^6\right]y_s, \]

or finally

\[ y_s'=y'(x_s)=\left(\frac{dy}{dt}\right)_{y=y_s} =\frac{1}{h}\left(\Delta-\frac{1}{2}\Delta^2+\frac{1}{3}\Delta^3-\frac{1}{4}\Delta^4+\frac{1}{5}\Delta^5-\frac{1}{6}\Delta^6\right)y_s=Dy_s. \tag{16} \]

\[ \begin{aligned} 2)\quad y_s'&=y'(x_s)=\lim_{\xi\to 0}\frac{y(x_s+\xi)-y(x_s)}{\xi} =\lim_{t\to 0}\frac{y(x_s+th)-y(x_s)}{th}=\\ &=\lim_{t\to 0}\frac{1}{th}\left[ t\mu\delta y_s+\frac{t^2}{2!}\delta^2 y_s+\frac{t(t^2-1)}{3!}\mu\delta^3 y_s +\frac{t^2(t^2-1)}{4!}\delta^4 y_s+\right.\\ &\qquad\left. +\frac{t(t^2-1)(t^2-4)}{5!}\mu\delta^5 y_s +\frac{t^2(t^2-1)(t^2-4)}{6!}\delta^6 y_s \right] \end{aligned} \]

or finally

\[ y_s'=y'(x_s)=\frac{1}{h}\left(\mu\delta y_s-\frac{1}{6}\mu\delta^3 y_s+\frac{1}{30}\mu\delta^5 y_s\right)=Dy_s. \tag{17} \]

\[ \begin{aligned} 3)\quad y_s'&=y'(x_s)=\lim_{\xi\to 0}\frac{y(x_s-\xi)-y(x_s)}{-\xi} =\lim_{t\to 0}\frac{y(x_s-th)-y(x_s)}{-th}=\\ &=\lim_{t\to 0}\frac{1}{-th}\left[ -t\nabla+\frac{t(t-1)}{2!}\nabla^2-\frac{t(t-1)(t-2)}{3!}\nabla^3+\cdots+ \frac{t(t-1)\cdots(t-5)}{6!}\right]y_s. \end{aligned} \]

or finally

\[ y_s'=y'(x_s)=\left(\frac{dy}{dt}\right)_{y=y_s} =\frac{1}{h}\left(\nabla+\frac{1}{2}\nabla^2+\frac{1}{3}\nabla^3+\frac{1}{4}\nabla^4+\frac{1}{5}\nabla^5+\frac{1}{6}\nabla^6\right)y_s=Dy_s, \tag{18} \]

In order to pass from the expansions (16)—(18) for \(y_s'\) to expansions for \(y_s''\), we introduce the differentiation operator \(D\):

\[ D=\frac{d}{dt}, \tag{19} \]

whence it further follows that

\[ \frac{d^x}{dt^x}=\left(\frac{d}{dt}\right)^x=D^x, \]

and, therefore,

\[ y_s^{(x)}=\left(\frac{d^x y}{dt^x}\right)_{y=y_s}=D^x y_s. \tag{20} \]

Thus, in particular,

\[ y_s''=\left(\frac{d^2y}{dt^2}\right)_{y=y_s}=D^2y_s=(DD)y_s. \tag{21} \]

Proceeding now from the expansions (16)—(18) for \(y'_s=Dy_s\) and using the differentiation operator \(D\), we obtain the following three expressions for \(y''_s\), accurate up to terms of the 6th order of smallness:

\[ 1)\quad y''_s=y''(x_s)=D^2y_s=\frac{1}{h^2}\left(\Delta-\frac{1}{2}\Delta^2+\frac{1}{3}\Delta^3-\frac{1}{4}\Delta^4+\frac{1}{5}\Delta^5-\frac{1}{6}\Delta^6\right)^2 y_s, \]

or, finally,

\[ y''_s=\frac{1}{h^2}\left(\Delta^2-\Delta^3+\frac{11}{12}\Delta^4-\frac{5}{6}\Delta^5+\frac{137}{180}\Delta^6\right)y_s; \tag{21} \]

\[ 2)\quad y''_s=y''(x_s)=D^2y_s=\frac{1}{h^2}\left(\mu\delta-\frac{1}{6}\mu\delta^3+\frac{1}{30}\mu\delta^5\right)^2y_s= \]

\[ =\frac{1}{h^2}\left(\mu^2\delta^2-\frac{1}{3}\mu^2\delta^4+\frac{17}{180}\mu^2\delta^6\right)y_s; \tag{*} \]

\[ 3)\quad y''_s=y''(x_s)=D^2y_s=\frac{1}{h^2}\left(\nabla+\frac{1}{2}\nabla^2+\frac{1}{3}\nabla^3+\frac{1}{4}\nabla^4+\frac{1}{5}\nabla^5+\frac{1}{6}\nabla^6\right)^2 y_s. \]

or, finally,

\[ y''_s=\frac{1}{h^2}\left(\nabla^2+\nabla^3+\frac{11}{12}\nabla^4+\frac{5}{6}\nabla^5+\frac{137}{180}\nabla^6\right)y_s. \tag{22} \]

The finite-difference expressions (21), (22) for \(y''_s\) have final form. The expression (*) for \(y''_s\), however, we transform by expressing \(\mu^2\) in terms of \(\delta^2\) with the help of relations (9.5) and (9.7). Relying on these relations, for \(h=1\) we obtain:

\[ 1)\quad \mu=\frac{1}{2}\left(E^{\frac12}+E^{-\frac12}\right),\qquad 2)\quad \delta=E^{\frac12}-E^{-\frac12},\qquad 3)\quad \delta^2=E+E^{-1}-2, \]

\[ 4)\quad \mu^2=\frac{1}{4}\left(E+E^{-1}+2\right)=\frac{1}{4}\left[(E+E^{-1}-2)+4\right]= \]

\[ =\frac{1}{4}(E+E^{-1}-2)+1=\frac{1}{4}\delta^2+1. \]

Consequently,

\[ \mu^2=\frac{1}{4}\delta^2+1. \tag{23} \]

Substituting further in (*) the factor \(\mu^2\) by its expression (23) in terms of \(\delta^2\), we find one more finite-difference expression for \(y''_s\):

\[ y''_s=y''(x_s)=\frac{1}{h^2}\left[\delta^2-\frac{1}{12}\delta^4+\frac{1}{90}\delta^6\right]y_s. \tag{24} \]

In conclusion, we need the obtained finite-difference expansions (16)—(18) and (21), (24), (22) for \(y'_s\) and \(y''_s\) to be expressed directly in terms of \(y_x\) \((x=0,1,2,\ldots,n)\). This can be done with the help of

relations (9.1)—(9.3) for the operators \(E, \Delta, \nabla, \delta\), which for \(h=1\) give

\[ 1)\ \Delta E-1,\qquad 2)\ \delta=E^{\frac12}-E^{-\frac12},\qquad 3)\ \Delta=1-E^{-1}. \tag{25} \]

From these relations it further follows that

\[ 1)\ \Delta^\nu y_s=(E-1)^\nu y_s,\qquad 2)\ \delta^\nu y_x=\left(E^{\frac{\nu}{2}}-E^{-\frac{\nu}{2}}\right)y_x,\qquad 3)\ \nabla^\nu y_s=(1-E^{-1})y_s \tag{26} \]

If one also takes into account the meaning of the transformation (8.3), performed by the operator \(E\) on \(y_s, y_x\) for \(h=1\),

\[ 1)\ E^\lambda y_s=y_{s+\lambda},\qquad 2)\ E^{\pm\frac{\mu}{2}}y_x=y_{x\pm\frac12\mu},\qquad 3)\ E^{-\lambda}y_s=y_{s-\lambda}, \tag{27} \]

then we obtain the following expressions for \(y_s'\), \(y_s''\) in terms of \(y_x\), \((x=0,1,2,\ldots,n)\), following from (16) and (21), (17) and (24), (18) and (22):

\[ \begin{cases} \displaystyle 1)\ y_s'=\frac{1}{60h}\bigl[-147y_s+360y_{s+1}-450y_{s+2}+400y_{s+3}\\ \displaystyle\hspace{3.2cm} -225y_{s+4}+72y_{s+5}-10y_{s+6}\bigr], \\[1.2em] \displaystyle 2)\ y_s''=\frac{1}{180h^2}\bigl[812y_s-3132y_{s+1}+5265y_{s+2}-5080y_{s+3}\\ \displaystyle\hspace{3.2cm} +2970y_{s+4}-972y_{s+5}+137y_{s+6}\bigr]. \end{cases} \tag{28} \]

\[ \begin{cases} \displaystyle 1)\ y_s'=\frac{1}{60h}\bigl[y_{s+3}-9y_{s+2}+45y_{s+1}-45y_{s-1}+9y_{s-2}-y_{s-3}\bigr], \\[1.2em] \displaystyle 2)\ y_s''=\frac{1}{180h^2}\bigl[2y_{s+3}-27y_{s+2}+270y_{s+1}-490y_s\\ \displaystyle\hspace{3.2cm} +270y_{s-1}-27y_{s-2}+2y_{s-3}\bigr]. \end{cases} \tag{29} \]

\[ \begin{cases} \displaystyle 1)\ y_s'=-\frac{1}{60h}\bigl[147y_s-360y_{s-1}+450y_{s-2}-400y_{s-3}\\ \displaystyle\hspace{3.2cm} +225y_{s-4}-72y_{s-5}+10y_{s-6}\bigr], \\[1.2em] \displaystyle 2)\ y_s''=\frac{1}{180h^2}\bigl[812y_s-3132y_{s-1}+5265y_{s-2}-5080y_{s-3}\\ \displaystyle\hspace{3.2cm} +2970y_{s-4}-972y_{s-5}+137y_{s-6}\bigr]. \end{cases} \tag{30} \]

We apply expansions (28) for calculating \(y_s'\) and \(y_s''\) when \(s=0,1,2\). With the aid of expansions (29) we express \(y_s\) and \(y_s'\) for \(s=3,4,\ldots,n-3\). Finally, expansions (30) serve us for calculating \(y_s'\) and \(y_s''\) when \(s=n-2,n-1,n\) (Fig. 1).

Substitution of expansions (28)—(30) for \(y_s'\) and \(y_s''\) into system (5) of \((n+1)\) equations with \(3(n+1)\) unknowns \(y_s, y_s', y_s''\) transforms it into a system of \((n+1)\) first-degree algebraic equations with \(n+1\) unknowns \(y_x\) \((x=0,1,2,\ldots,n)\):

\[ \left\{ \begin{aligned} 1)\;& \sum_{x=0}^{n} b_{1x} y_x = A,\\ 2)\;& \sum_{x=0}^{n} b_{sx} y_x = 0,\qquad (s=1,2,\ldots,n-1),\\ 3)\;& \sum_{x=n}^{0} b_{nx} y_x = B. \end{aligned} \right. \tag{31} \]

Since expressions (28)—(30) for the derivatives \(y_s'\), \(y_s''\) are given with a high degree of accuracy, the values \(y_s\) will also be obtained, with approximately the same high degree of accuracy, from the solution of the system of equations (31) on a digital computer. For this same reason, the values \(y_s\) from the solution of the system (31) will be found with sufficient accuracy even with a comparatively sparse grid of the corresponding points \(x_s\) on the \(Ox\) axis.

REFERENCES

1. E. D. But. Numerical Methods. Moscow, Fizmatgiz, 1959.
2. E. Whittaker and G. Robinson. Mathematical Processing of the Results of Observations. Moscow–Leningrad, ONTI, 1935.
3. B. P. Demidovich et al. Numerical Methods of Analysis. Moscow, Fizmatgiz, 1963.
4. L. Collatz. Eigenvalue Problems. Moscow, “Nauka,” 1968.