Application of the Adams-Störmer Method to Higher-Order Equations
Abstract Generated abstract
This paper develops a direct extension of the Adams-Störmer numerical integration method for higher order ordinary differential equations, especially equations of the form y to the mth derivative equals f(x,y). Instead of first reducing such equations to systems of first or second order equations, the method constructs difference tables for the function and for the highest derivative, then derives recurrence formulas for advancing the solution using lower diagonal differences. The paper proves that the resulting multipliers are independent of the step index, gives systems and formulas for computing them, and lists practical cases for low orders. Two worked examples illustrate the procedure and indicate that it can reduce computational work while maintaining good numerical agreement with known or tabulated solutions.
Full Text
APPLICATION OF THE ADAMS–STÖRMER METHOD TO HIGHER-ORDER EQUATIONS
A. Kh. Pinsky
Equation of the form \(y^{(m)}=f(x,y)\)
§ 1. The Adams–Störmer method, when applied to higher-order equations, requires a preliminary transformation of the given equation in order to represent it in the form of a system of equations, the order of each of which does not exceed two. Thus, the equation
\[ y'''=f(x,y) \]
by the substitution \(y'=z\) is reduced to a system of equations, one of which is of the second order.
The system is solved by the Adams–Störmer method; however, the appearance of a new variable, as well as the use of two formulas, introduces a superfluous element into the process of solving the original equation. There will be more such “superfluous elements” the higher the order of the equation, if it is reduced to a system. The method of numerical integration of differential equations proposed here makes it possible to solve the given equation either without any preliminary transformations, or by first reducing it to a system that is most conveniently solved.
The method is based on using differences of the values of the derivative of highest order, arranged along the lower oblique row, i.e., on using the same elements as in the Adams–Störmer method. Thus, the proposed method is a further development of the Adams–Störmer method.
Let us first consider the equation
\[ y^{(m)}=f(x,y) \tag{1} \]
Let the particular integral of this equation corresponding to certain initial conditions \((x_0,y_0,y'_0,\ldots,y_0^{(m-1)})\) be
\[ y=F(x). \]
Suppose further that we know \((n+1)\) values of this integral
\[ y_0,y_1,y_2,\ldots,y_n \qquad (n\geq m), \]
corresponding to the values of the argument \(x\):
\[ x_0,\ x_1=x_0+h,\ x_2=x_0+2h,\ldots,x_n=x_0+nh, \]
and it is required to find \(y_{n+1}\).
To solve this problem we proceed as follows. From the \((n+1)\) values of the function available to us, we form a table of differences up to order \(m\) inclusive.
Introducing the notation
\[
y_n=F(x_0+nh),\qquad y_{n+1}-y_n=\delta y_n
\]
\[
\delta y_{n+1}-\delta y_n=\delta^2y_n\quad\text{and so on,}
\]
we shall present this table in the following form (Table 1). Further, by the definition of differences, we have
Table No. 1.
| \(y\) | \(\delta y\) | \(\delta^2 y\) | \(\delta^3 y\) | \(\cdots\) | \(\delta^{m-1}y\) | \(\delta^m y\) |
|---|---|---|---|---|---|---|
| \(y_0\) | \(\delta y_0\) | |||||
| \(y_1\) | \(\delta y_1\) | \(\delta^2y_0\) | ||||
| \(y_2\) | \(\delta y_2\) | \(\delta^2y_1\) | \(\delta^3y_{n-3}\) | |||
| \(y_3\) | — | — | ||||
| — | — | — | — | |||
| — | — | — | — | \(\delta^m y_0\) | ||
| — | — | \(\delta^2y_{m-2}\) | — | — | — | — |
| \(y_m\) | \(\delta y_{m-1}\) | — | ||||
| \(y_{m+1}\) | \(\delta y_m\) | \(\delta^2y_m\) | — | |||
| — | — | — | — | — | ||
| — | — | — | — | — | — | |
| — | — | — | — | \(\delta^m y_{n-m}\) | ||
| \(y_{n-3}\) | \(\delta y_{n-3}\) | \(\delta^{m-1}y_{\,n-(m-1)}\) | ||||
| \(y_{n-2}\) | \(\delta y_{n-2}\) | \(\delta^2y_{n-3}\) | \(\delta^3y_{n-3}\) | |||
| \(y_{n-1}\) | \(\delta y_{n-1}\) | \(\delta^2y_{n-2}\) | ||||
| \(y_n\) |
\[ y_{n+1}=y_n+\delta y_n \]
\[ \delta y_n=\delta y_{n-1}+\delta^2y_{n-1} \]
\[ \delta^2y_{n-1}=\delta^2y_{n-2}+\delta^3y_{n-2} \tag{2} \]
\[ \ldots\ldots\ldots\ldots \]
\[ \ldots\ldots\ldots\ldots \]
\[ \delta^{m-2}y_{\,n-(m-3)} = \delta^{m-2}y_{\,n-(m-2)} + \delta^{m-1}y_{\,n-(m-2)} \]
\[ \delta^{m-1}y_{\,n-(m-2)} = \delta^{m-1}y_{\,n-(m-1)} + \delta^m y_{\,n-(m-1)} \]
Adding term by term, we shall have
\[ y_{n+1}=y_n+\delta y_{n-1}+\delta^2 y_{n-2}+\cdots+\delta^{m-1}y_{n-(m-1)}+\delta^m y_{n-(m-1)} \tag{2'} \]
On the right-hand side of equality \((2')\) all terms except the last are known, since they are located in the lower diagonal row of the table of differences; consequently, if we found the term
\[ \delta^m y_{n-(m-1)}, \]
then the value \(y_{n+1}\) would be found by addition.
Let us now suppose that, for the same values of the argument \(x\) for which the values of the function are known, the values of the derivative of order \(m\) are also known:
\[ y_0^{(m)}, y_1^{(m)}, y_2^{(m)}, \ldots \ldots \ldots y_n^{(m)} . \]
Multiply each of these values by \(h^m\) and form from them a table of differences, assuming that their difference of order \(k\) will be constant. Introducing the notation
\[ h^m f(x_i,y_i)=\Delta_0^m y_i, \]
\[ \Delta_0^m y_{i+1}-\Delta_0^m y_i=\Delta_1^m y_i, \]
\[ \Delta_1^m y_{i+1}-\Delta_1^m y_i=\Delta_2^m y_i \]
and so on,
we shall present this table in the following form (Table 2).* As in the Adams–Störmer method, we shall seek the value
\[ \delta^m y_{n-(m-1)} \]
through the differences of the lower diagonal row of the second table, by means of the equality
\[ \delta^m y_{n-(m-1)}=a_{0m}y_n+a_{1m}\Delta_1^m y_{n-1}+a_{2m}\Delta_2^m y_{n-2}+\ldots \]
\[ +\ldots \ldots +a_{im}\Delta_i^m y_{n-i}+\ldots \ldots +a_{km}\Delta^m y_{n-k} \tag{4} \]
where \(a_{im}\) are certain multipliers.
To determine the multipliers \(a_{im}\), we transform equality \((4)\). Putting
\[ x=x_0+nh, \]
where \(n\) is a new variable, and expanding the integral function
\[ y=F(x) \]
in a Taylor series, we shall have
\[ F(x_0+nh)=F(x_0)+\frac{nh}{1}F'(x_0)+\frac{n^2h^2}{2!}F''(x_0)+\frac{n^3h^3}{3!}F'''(x_0)+\cdots+ \]
\[ +\cdots+\frac{n^m h^m}{m!}F^{(m)}(x_0)+\cdots+\frac{n^{m+k}h^{m+k}}{(m+k)!}F^{(m+k)}(x_0)+\cdots \tag{5} \]
Putting, in general,
\[ y^{(i)}h^i=\Delta_0^i y_0 \qquad i=1,2,3,\ldots m,\ldots (m+k),\ldots \]
we rewrite \((5)\) in the form:
\[ y_n=y_0+\frac{n}{1}\Delta_0^2 y_0+\frac{n^2}{2!}\Delta_0^2 y_0+\frac{n^3}{3!}\Delta_0^3 y_0+\ldots+ \]
\[ +\ldots+\frac{n^m}{m!}\Delta_0^m y_0+\ldots \ldots+\frac{n^{m+k}}{(m+k)!}\Delta_0^{m+k}y_0+\ldots \tag{6} \]
* Table 2, see p. 11.
From the last equality, substituting for \(n\) the values \(1,2,3,\ldots,n,(n+1)\) and taking differences up to the \(m\)-th order inclusive, we find an expression for the left-hand side of (4) in the following form:
\[
\delta^{m}y_{n-(m-1)}=
\frac{A_m}{m!}\Delta_0^m y_0+
\frac{A_{m+1}}{(m+1)!}\Delta_0^{m+1}y_0+
\frac{A_{m+2}}{(m+2)!}\Delta_0^{m+2}y_0+\cdots+
\]
\[
+\cdots+
\frac{A_{m+i}}{(m+i)!}\Delta_0^{m+i}y_0+\cdots+
\frac{A_{m+k}}{(m+k)!}\Delta_0^{m+k}y_0
\tag{7}
\]
where
\[ A_{m+i}=(n+1)^{m+i}-mn^{m+i} +\frac{m(m-1)}{2!}(n-1)^{m+i}+\cdots\pm[n-(m-1)]^{m+i}. \]
The right-hand side of (4) is transformed analogously. Differentiating (6) \(m\) times with respect to \(n\) and multiplying by \(h^m\), we shall have
\[ \Delta_0^m y'_n=\Delta_0^m y_0+ \frac{n}{1}\Delta_0^{m+2}y_0+ \frac{n^2}{2}\Delta_0^{m+2}y_0+\cdots+ \frac{n^k}{k!}\Delta_0^{m+k}y_0 \tag{8} \]
Hence, substituting for \(n\) the values \(1,2,3,\ldots,n\) and taking differences up to the \(k\)-th order inclusive, we obtain
\[
\Delta_0^m y_n=\Delta_0^m y_0\frac{n}{1}
+\Delta_0^{m+1}y_0+
\frac{n^2}{2}\Delta_0^{m+2}y_0+
\frac{n^3}{3!}\Delta_0^{m+3}y_0+\cdots+
\]
\[
+\cdots+
\frac{n^i}{i!}\Delta_0^{m+i}y_0+\cdots+
\frac{n^k}{k!}\Delta_0^{m+k}y_0
\]
\[
\Delta_1^m y_{n-1}=
\frac{A_{11}}{1!}\Delta_0^{m+1}y_0+
\frac{A_{12}}{2!}\Delta^{m+1}y_0+
\frac{A_{13}}{3!}\Delta_0^{m+3}y_0+\cdots+
\]
\[
+\cdots+
\frac{A_{1i}}{i!}\Delta_0^{m+i}y_0+\cdots+
\frac{A_{1k}}{k!}\Delta_0^{m+k}y_0
\]
\[
\Delta_2^m y_{n-2}=
\frac{A_{22}}{2!}\Delta_0^{m+2}y_0+
\frac{A_{23}}{3!}\Delta^{m+3}y_0+
\frac{A_{24}}{4!}\Delta_0^{m+4}y_0+
\]
\[
+\cdots+
\frac{A_{2i}}{i!}\Delta_0^{m+i}y_0+\cdots+
\frac{A_{2k}}{k!}\Delta_0^{m+k}y_0.
\]
\[
\Delta_3^m y_{n-3}=
\frac{A_{33}}{3!}\Delta_0^{m+3}y_0+
\frac{A_{34}}{3!}\Delta_0^{m+4}y_0+\cdots+
\]
\[
+\cdots+
\frac{A_{3i}}{i!}\Delta_0^{m+1}y_0+
\frac{A_{3k}}{k!}\Delta_0^{m+k}y_0
\tag{9}
\]
\[ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot \]
\[ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot \]
\[ \Delta_i^m y_{n-i}= \frac{A_{ii}}{i!}\Delta_0^{m+i}y_0+ \frac{A_{i,i+1}}{(i+1)!}\Delta_0^{m+i+1}y_0+\cdots+ \frac{A_{ik}}{k!}\Delta_0^{m+k}y_0 \]
\[ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot \]
\[ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot \]
\[ \Delta_k^m y_{n-k}= \frac{A_{kk}}{k!}\Delta_0^{m+k}y_0 \]
where
\[
A_{ik}=n^k-i(n-1)^k+
\frac{i(i-1)}{2!}(n-2)^k-
\frac{i(i-1)(i-2)}{3!}(n-3)^k+\cdots+
\]
\[
+\cdots\pm(n-i)^k,\qquad i=1,2,3,\ldots,k.
\]
and, according to the notation introduced,
\[ \Delta_0^{m+i} y_0 = y_0^{(m+i)} h^{m+i}. \]
Multiplying the first of equalities (9) by \(\alpha_{0m}\), the second by \(\alpha_{1m}\), and so on, and substituting into the right-hand side of (4). Replacing the left-hand side of the same equality (4) by the right-hand side of (7) and comparing the coefficients of identical constants
\[ \Delta_0^i y_0 \]
in the left- and right-hand sides, we obtain the following system of equations relating the multipliers \(\alpha_{im}\):
\[ \alpha_{0m} A_{00}=\frac{A_m}{m!} \]
\[ A_{01}\alpha_{0m}+A_{11}\alpha_{1m}=\frac{A_{m+1}}{(m+1)!} \]
\[ \frac{1}{2!}\left(A_{02}\alpha_{0m}+A_{12}\alpha_{1m}+A_{22}\alpha_{2m}\right)=\frac{A_{m+2}}{(m+2)!} \]
\[ \frac{1}{3!}\left(A_{03}\alpha_{0m}+A_{13}\alpha_{1m}+A_{23}\alpha_{2m}+A_{33}\alpha_{3m}\right)=\frac{A_{m+3}}{(m+3)!} \]
\[ \cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots \]
\[ \cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots \]
\[ \frac{1}{i!}\left(A_{0i}\alpha_{0m}+A_{1i}\alpha_{1m}+A_{2i}\alpha_{2m}+\cdots+A_{ii}\alpha_{im}\right)+\frac{A_{m+i}}{(m+i)!} \]
\[ \cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots \]
\[ \cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots \]
\[ \frac{1}{(k-1)!}A_{0,k-1}\alpha_{0m}+A_{1,k-1}\alpha_{1m}+\cdots+A_{i,k-1}\alpha_{im}+\cdots+A_{k-1,k-1}\alpha_{k-1,m} =\frac{A_{m+k-1}}{(m+k-1)!} \]
\[ \frac{1}{k!}\left(A_{0k}\alpha_{0m}+A_{1k}\alpha_{1m}+A_{2k}\alpha_{2m}+\cdots\cdots\cdots\cdots\cdots+A_{kk}\alpha_{km}\right) =\frac{A_{m+k}}{(m+k)!} \]
where, for symmetry, \(A_{0i}=n^i\) has been put.
§ 2. In equations (10) the coefficients
\[ A_{ik}=n^k-i(n-1)^k+\frac{i(i-1)}{2}(n-2)^k+\cdots\pm(n-i)^k \tag{11} \]
and
\[ A_{m+i}=(n+1)^{m+i}-mn^{m+i}+\frac{m(m-1)}{2}(n-1)^{m+i}+\cdots\pm[n-(m-1)]^{m+i} \]
are functions of \(n\), and at first glance it seems that the multipliers \(\alpha_{im}\) are likewise functions of \(n\).
Let us prove that the multipliers \(\alpha_{im}\) do not depend on \(n\).
First of all, note that from (11) the following relation follows:
\[ \frac{d}{dn}\left(\frac{A_{ik}}{k!}\right) = \frac{ k\left[n^{k-1}-i(n-1)^{k-1}+\frac{i(i-1)}{2}(n-2)^{k-1}+\cdots\pm(n-i)^{k-1}\right] }{k!} = \frac{A_{i,k-1}}{(k-1)!} \]
and
\[ \frac{d}{dn}\left(\frac{A_{m+i}}{(m+i)!}\right) = \frac{A_{m-i-1}}{(m+i-1)!}. \]
We also note that if \(i=k\), then
\[ -\frac{d}{dn}\left(\frac{A_{ik}}{k!}\right)=0, \]
because the expression \(A_{ik}\) is the difference of order \(i\) of the function \(n^k\), and if \(i=k\), then \(A_{ik}=k!\).
Likewise \(A_m=m!\), whence it follows that \(a_{0m}=1\).
Let us now put \(a_{im}=\varphi_i(n)\), where \(\varphi_i(n)\) is some function of \(n\), and substitute into the last two equations of system (10). We have:
\[ \frac{1}{(k-1)!}\left[ A_{0,k-1}\varphi_0(n)+A_{1,k-1}\varphi_1(n)+\cdots +A_{i,k-1}\varphi_i(n)+\cdots+ \right. \]
\[ \left. +\cdots+A_{k-1,k-1}\varphi_{k-1}(n) \right] = \frac{A_{m+k-1}}{(m+k-1)!} \]
\[ \frac{1}{k!}\left[ A_{0,k}\varphi_0(n)+A_{0,k}\varphi_1(n)+\cdots+ A_{k-1}\varphi_{k-1}(n)+A_{k,k}\varphi_k(n) \right] = \frac{A_{m+k}}{(m+k)!}. \]
Differentiating the second equality with respect to \(n\) and subtracting the upper one, we shall have, taking into account the remarks made above:
\[ A_{0k}\varphi_0'(n)+A_{1k}\varphi_1'(n)+\cdots+ A_{ik}\varphi_i'(n)+A_{kk}\varphi_k'(n)=0, \]
whence, successively setting \(k=1,2,3\ldots k\), we find
\[ \varphi_i'(n)=0, \]
therefore
\[ \varphi_i(n)=C. \]
The last relation makes it possible to simplify system (10) considerably by a choice of \(n\). Thus, putting \(n=0\) and taking into account that, for \(n=0\), \(A_{0i}=0\), we rewrite system (10) in the following form:
\[ C_{11}a_{1m}=\frac{C_{m+1}}{(m+1)!} \]
\[ \frac{1}{2!}(C_{12}a_{1m}+C_{22}a_{2m})= \frac{C_{m+2}}{(m+2)!} \]
\[ \frac{1}{3!}(C_{13}a_{1m}+C_{23}a_{2m}+C_{33}a_{3m})= \frac{C_{m+3}}{(m+3)!} \tag{12} \]
\[ \frac{1}{4!}(C_{14}a_{1m}+C_{24}a_{2m}+C_{34}a_{3m}C_{44}a_{4m})= \frac{C_{m+4}}{(m+4)} \]
\[ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot \]
\[ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot \]
\[ \frac{1}{i!}(C_{1i}a_{1m}+C_{2i}a_{2m}+\cdots+C_{ii}a_{im})= \frac{C_{m+i}}{(m+i)!} \]
\[ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot \]
\[ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot\ \cdot \]
\[ \frac{1}{k!}(C_{1k}a_{1m}+C_{2k}a_{2m}+\cdots+C_{ik}a_{im}+\cdots+C_{kk}a_{kk})= \frac{C_{m+k}}{(m+k)!}, \]
where
\[ C_{ik}=A_{ik}\big|_{n=0} = (-1)^{i+k} \left[ i^k-i(i-1)^k+\frac{i(i-1)}{2}(i-2)^k+\cdots \pm i \right]. \tag{12'} \]
and
\[ C_{m+i}=A_{m+i}\big|_{n=0}=1+(-1)^{i+m}\left[ \frac{m(m-1)}{2!}\,1^{m+i} -\frac{m(m-1)(m-2)}{3!}\,2^{m+i} +\cdots\pm (m-1)^{m+i} \right] \]
The solution of system (12) in general form, from the formal point of view, presents no difficulties. Stopping at the determination of the factor \(a_{im}\), we shall have
\[ a_{ij}=\frac{\Delta'}{\Delta} \tag{13}, \]
where
\[ \Delta= \left| \begin{array}{ccccc} C'' & 0 & 0\ldots 0 & 0\\ C_{12} & C_{22} & 0\ldots 0 & 0\\ \ldots & \ldots & \ldots & \ldots & \ldots\\ \ldots & \ldots & \ldots & \ldots & \ldots\\ C_{1j} & C_{2j} & \ldots & \ldots & \ldots\\ \ldots & \ldots & \ldots & \ldots & \ldots\\ \ldots & \ldots & \ldots & \ldots & \ldots\\ C_{1,i-1} & C_{2,i-1} & \ldots & C_{i-1,i-1} & 0\\ C_{1,i} & C_{2,i} & \ldots & C_{i-1,i} & C_{ii} \end{array} \right| \]
is the determinant of the system, and \(\Delta'\) is the determinant obtained by replacing the elements of the \(i\)-th column of the determinant \(\Delta\) by the corresponding terms of the right-hand side of system (12).
Expanding \(\Delta'\) with respect to the elements of the \(i\)-th column, we represent the expression \(a_{im}\) in the following form:
\[
a_{im}=\frac{1}{\Delta}\left[
\frac{C_{m+1}i!}{(m+i)!}\Delta_{ii}
-\frac{C_{m+i-1}(i-1)!}{(m+i-1)!}\Delta_{i-1,i}
+\frac{C_{m+i-2}(i-2)!}{(m+i-2)!}\Delta_{i-2,i}
+\mp\cdots \ldots
\right.
\]
\[
\left.
+\cdots+(-1)^{i+j}\frac{C_{m+i}j!}{(m+j)!}\Delta_{ji}
+\cdots+(-1)^{i+1}\frac{C_{m+1}}{(m+1)!}\Delta_{1i}
\right]
\tag{14},
\]
where \(\Delta_{ji}\) is the minor corresponding to the element of the \(i\)-th column and \(j\)-th row of the determinant \(\Delta\).
We would obtain the same expression by successive substitution.
§ 3. For the arithmetic calculation of the factors \(a_{im}\), formula (14) is inconvenient. The calculation can be carried out considerably more quickly by solving system (12) for particular values of \(m\), especially since the left-hand side of the system does not depend on the order of the equation.
Let us note that the coefficients \(C_{ik}\) have the following properties:
\[ \sum_{i=1}^{i=k} C_{ih}=1 \]
\[ C_{ik}+i\left(C_{i,k-1}-C_{i-1,k-1}\right)=0 \]
\[ \sum_{i=k}^{i=k}\frac{C_{ik}}{i+1}=B_k \quad \text{for even } k \tag{14'} \]
\[ \sum_{i=1}^{i=1}\frac{C_{ik}}{i+1}=0 \quad \text{for odd } k, \]
where \(B_k\) is the \(k\)-th Bernoulli number.
One of the properties of (14′), namely the second, makes it possible to find the coefficients \(C_{ik}\) very quickly.
Let us consider several particular cases of the system (12). Having computed the values \(C_{ik}\) for \(k=1,2\ldots 7\) and \(i=1,2\ldots k\), and putting \(m=1\), we rewrite (12) in the form
\[ \alpha_{11}=\frac{1}{2}, \]
\[ -\alpha_{11}+2\alpha_{21}=\frac{1}{3} \]
\[ \alpha_{11}-6\alpha_{21}+6\alpha_{31}=\frac{1}{4} \]
\[ -\alpha_{11}+14\alpha_{21}-36\alpha_{31}+24\alpha_{41}=\frac{1}{5} \]
\[ \alpha_{11}-30\alpha_{21}+150\alpha_{31}-240\alpha_{41}+120\alpha_{51}=\frac{1}{6} \]
\[ -\alpha_{11}+62\alpha_{21}-540\alpha_{31}+1560\alpha_{41}-1800\alpha_{51}+720\alpha_{61}=\frac{1}{7} \]
\[ \alpha_{11}-126\alpha_{21}+1806\alpha_{31}-8400\alpha_{41}+16800\alpha_{51}-15120\alpha_{61}+5040\alpha_{71}=\frac{1}{8} \]
whence we find
\[ \alpha_{11}=\frac{1}{2},\quad \alpha_{21}=\frac{5}{12},\quad \alpha_{31}=\frac{3}{8},\quad \alpha_{41}=\frac{251}{720}, \]
\[ \alpha_{51}=\frac{95}{288},\quad \alpha_{61}=\frac{19087}{60480},\quad \alpha_{71}=\frac{5257}{17280}. \]
Substituting in (4) for \(m=1\) and remembering that \(\alpha_{0m}=1\), we shall have
\[ \delta y_n=\Delta'_0 y_n+\frac{1}{2}\Delta'_1 y_{n-1} +\frac{5}{12}\Delta'_2 y_{n-2} +\frac{3}{8}\Delta'_3 y_{n-3} +\frac{251}{720}\Delta'_4 y_{n-4} \]
\[ +\frac{95}{288}\Delta'_5 y_{n-5} +\frac{19087}{60480}\Delta'_6 y_{n-6} +\frac{5257}{17280}\Delta'_7 y_{n-7} \tag{16} \]
This is the Adams formula.
For \(m=2\) the system is rewritten as follows:
\[ \alpha_{11}=0 \]
\[ -\alpha_{12}+2\alpha_{22}=\frac{1}{6} \]
\[ \alpha_{12}-6\alpha_{22}+6\alpha_{32}=0 \]
\[ -\alpha_{12}+14\alpha_{22}-36\alpha_{32}+24\alpha_{42}=\frac{1}{12} \]
\[ \alpha_{12}-30\alpha_{22}+150\alpha_{32}-240\alpha_{42}+120\alpha_{52}=0 \]
\[ -\alpha_{12}+62\alpha_{22}-540\alpha_{32}+1560\alpha_{42}-1800\alpha_{52}+720\alpha_{62}=\frac{1}{28} \]
\[ \alpha_{12}-126\alpha_{22}+1806\alpha_{32}-8400\alpha_{42}+16800\alpha_{52}-15120\alpha_{62}+5040\alpha_{72}=0 \]
whence
\[ \alpha_{12}=0,\quad \alpha_{12}=\frac{1}{12},\quad \alpha_{32}=\frac{1}{12},\quad \alpha_{42}=\frac{19}{240}, \]
\[ \alpha_{52}=\frac{3}{40},\quad \alpha_{62}=\frac{863}{12096},\quad \alpha_{72}=\frac{275}{4032}. \]
and equality (14) takes the form
\[ \delta^2 y_{n-1} = \Delta_0^2 y_n + \frac{1}{12}\Delta_2^2 y_{n-2} + \frac{1}{12}\Delta_3^2 y_{n-3} + \frac{19}{240}\Delta_4^2 y_{n-4} + \]
\[ + \frac{3}{40}\Delta_5^2 y_{n-5} + \frac{863}{12096}\Delta_6^2 y_{n-6} + \frac{275}{4032}\Delta_7^2 y_{n-7} \tag{18} \]
The last formula is known as Störmer’s formula.
For \(m=3\), system (12) is rewritten as follows:
\[ \alpha_{13}=-\frac{1}{2} \]
\[ -\alpha_{13}+2\alpha_{23}=\frac{1}{2} \]
\[ \alpha_{13}-6\alpha_{23}+6\alpha_{33}=-\frac{1}{2} \]
\[ -\alpha_{13}+14\alpha_{23}-36\alpha_{33}+24\alpha_{4}=\frac{3}{5} \]
\[ \alpha_{13}-30\alpha_{23}+150\alpha_{33}-240\alpha_{43}+120\alpha_{53}=-\frac{3}{4} \]
\[ -\alpha_{12}+62\alpha_{23}-540\alpha_{33}+1560\alpha_{43}-1800\alpha_{53}+720\alpha_{63} = \frac{85}{84} \]
\[ \alpha_{13}-126\alpha_{23}+1806\alpha_{33}-8400\alpha_{43}+16800\alpha_{53} -15120\alpha_{63}+5040\alpha_{73} = -\frac{17}{12} \tag{19} \]
and the values of \(\alpha_{im}\) will be as follows:
\[ \alpha_{13}=\frac{1}{2}, \qquad \alpha_{23}=0, \qquad \alpha_{33}=0, \qquad \alpha_{43}=\frac{1}{240} \]
\[ \alpha_{53}=-\frac{1}{160}, \qquad \alpha_{63}=\frac{221}{30240}, \qquad \alpha_{73}=-\frac{95}{12096}. \]
Substituting into (4) for \(m=3\), we obtain
\[ \delta^3 y_{n-2} = \Delta_0^3 y_n - \frac{1}{2}\Delta_1^3 y_{n-1} + \frac{1}{240}\Delta_4^3 y_{n-4} + \frac{1}{160}\Delta_5^3 y_{n-5} + \]
\[ + \frac{221}{30240}\Delta_6^3 y_{n-6} +\cdots . \tag{20} \]
Without going into details, we give the result of solving system (12) for \(m=4\).
In this case
\[ \alpha_{14}=-1, \qquad \alpha_{24}=\frac{1}{2}, \qquad \alpha_{34}=0, \qquad \alpha_{44}=-\frac{1}{720} \]
and equality (4) takes the form:
\[ \delta^4 y_{n-3} = \Delta_0^4 y_n - \Delta_1^4 y_{n-1} + \frac{1}{6}\Delta_2^4 y_{n-2} - \frac{1}{720}\Delta_4^4 y_{n-4} + \tag{21} \]
§ 4. There is no need to stop at further particular values of \(m\). For every equation of the form
\[ y^{(m)}=f(x,y) \]
the system (12) and formula (4) give the corresponding formula for its solution, which is obtained by the same device by which we found formulas (16), (18), (20), and (21).
Let us dwell briefly on the error which we admit, assuming that the differences of the \(k\)-th order of the values
\[ i^{(m)} h^m = \Delta_0^m y_i \]
are constant. Accordingly, we expanded the function \(y\) in a Taylor series, restricting ourselves to the term of order \((m+k)\) with respect to \(h\); thus, the first of the omitted terms in the expansion of \(y\) into a series is
\[ \frac{n^{m+k+1} h^{m+k+1}}{(m+k+1)!}\,F^{(m+k+1)}(x_0) \]
Consequently, in expression (7), which determines the sought difference, the first of the omitted terms will be
\[ \frac{A_{m+k+1}}{(m+k+1)!}\,\Delta_0^{m+k+1} y_0 \]
and the corresponding term in (4) will be
\[ a_{k+1,m}\Delta_{k+1}^{m} y_{n-(k+1)} \]
But the last of equalities (9) gives
\[ \Delta_k^m y_{n-k}=\Delta_0^{m+k}y_0=y_0^{(m+k)}h^{m+k} \]
Thus, formula (4) is valid up to terms of order \((m+k)\) with respect to \(h\). A more detailed investigation of the error of formula (4) is a rather serious problem, on which there is no possibility of dwelling here.
§ 5. The equation \(y^{(m)}=f(x,y,y^{(p)})\)
For every equation of the form
\[ y^{(m)}=f(x,y,y^{(p)}) \tag{22} \]
equality (4) is unsuitable, since we do not know how to determine the derivative \(y^{(p)}\), on which the value of the derivative \(y^{(m)}\) depends.
Equation (2) will be solved if to formula (4), which determines \(y\), we add a formula determining the derivative \(y^{(p)}\). We shall now derive such a formula.
We shall assume that all the data which we used previously are known to us, i.e. that the \((n+1)\) values of the function \(y\) and of its derivative of order \(m\) are known, which in turn presupposes that the \(n+1\) values of the derivative \(y^{(p)}\) are also known to us:
\[ y_0^{(p)},\ y_1^{(p)},\ldots\ldots\ldots y_n^{(p)}, \tag{22} \]
since, not knowing these values, we cannot find the values of \(y^{(m)}\).
Thus, at our disposal there will be three groups of values:
\[ \begin{aligned} 1)&\quad y_0,\quad y_1,\quad y_2,\ldots\ldots y^n\\ 2)&\quad y_0^{(m)},\ y_1^{(m)},\ y_2^{(m)},\ldots\ldots y_n^{(m)}\\ 3)&\quad y_0^{(p)},\ y_1^{(p)},\ y_2^{(p)},\ldots\ldots y_n^{(p)} \end{aligned} \tag{22''} \]
which we shall use in the following way.
From the first group of values we shall form a table of differences up to order \(m\) inclusive. This will be our table No. 1, which we have already used. From the second group of values we shall form a table of differences up to order \(k\) inclusive. For convenience of reasoning, as before, we shall use the notation that we introduced, so that
the numerical value of the product of the derivative of order \(i\) by \(h^i\) will be denoted by the symbol
\[ \Delta_0^i y \]
regardless of the form of expression of the derivative.
Thus, the table composed from the second group of values will remain in the form in which it was clothed at the beginning, i.e., this will be Table No. 2.
Table No. 2
| \(\Delta_0^m y_0\) | |||||||
| \(\Delta_0^m y_1\) | \(\Delta_1^m y_0\) | ||||||
| \(\Delta_0^m y_2\) | \(\Delta_1^m y_1\) | \(\Delta_2^m y_0\) | |||||
| \(\Delta_0^m y_3\) | \(\Delta_1^m y_2\) | \(\Delta_2^m y_1\) | \(\Delta_3^m y_{n-3}\) | ||||
| — | |||||||
| — | |||||||
| \(\Delta_k^m y_0\) | |||||||
| — | — | ||||||
| \(\Delta_2^m y_{k-2}\) | — | — | |||||
| \(\Delta_0^m y_k\) | \(\Delta_1^m y_{k-1}\) | — | |||||
| \(\Delta_0^m y_{k+1}\) | \(\Delta_1^m y_k\) | \(\Delta_2^m y_k\) | — | ||||
| — | — | — | |||||
| — | — | ||||||
| — | \(\Delta_k^m y_{n-k}\) | ||||||
| \(\Delta_0^m y_{n-3}\) | — | ||||||
| \(\Delta_0^m y_{n-2}\) | \(\Delta_1^m y_{n-3}\) | \(\Delta_2^m y_{n-3}\) | — | ||||
| \(\Delta_0^m y_{n-1}\) | \(\Delta_1^m y_{n-2}\) | \(\Delta_2^m y_{n-2}\) | \(\Delta_3^m y_{n-3}\) | ||||
| \(\Delta_0^m y_n\) | \(\Delta_1^m y_{n-1}\) |
From the third group of values we shall compose a table of differences up to and including the \((m-p)\)-th order. This table will differ from the table of values of the function itself only in that, where the sign \(m\) stood, it is replaced by \((m-p)\), and instead of \(y\) there will stand \(y^{(p)}\).
The last column of this table will be \((m-p+1)\), and the last term of the lower oblique row will be
\[ \delta^{m-p} y_{\,n-(m-p)} \]
Suppose now that we need to find not only \(y_{n+1}\), but also \(y_{n+1}^{(p)}\). We know how to find \(y_{n+1}\). To determine \(y_{n+1}^{(p)}\), we turn to formulas (2) and \((2')\). Replacing the sign \(m\) by the sign \((m-p)\) and putting \(y^{(p)}\) in place of \(y\), we arrive at the conclusion that the value \(y_{n+1}^{(p)}\) is determined by filling in the lower diagonal row of the third table of differences; for this, in the expression
\[ y_{n+1}^{(p)}=y_n^{(p)}+\delta y_{n-1}^{(p)}+\delta^2 y_{n-2}^{(p)}+\cdots+\delta^{m-p-1}y_{n-(m-p-1)}^{(p)}+\delta^{m-p}y_{n-(m-p)1-} \]
it is necessary to know the last term, i.e. the term
\[ \delta^{m-p}y_{n-(m-p-1)}^{(p)} . \]
We shall determine this term through the differences of the lower diagonal row of the second table by the following equality:
\[ \delta^{m-p}y_{n-(m-p-1)}^{(p)} = \beta_{0p}\Delta_0^m y_n +\beta_{1p}\Delta_1^m y_{n-1} +\beta_{2p}\Delta_2^m y_{n-2} + \]
\[ +\cdots+\beta_{ip}\Delta_i^m y_{n-i}+\cdots+\beta_{kp}\Delta_k^m y_{n-1} \tag{23} \]
where \(\beta_{ip}\) are certain, as yet undetermined, multipliers. Equality (23) is no different from (4), and in order to determine the multipliers \(\beta_{ip}\) it is necessary to perform the same operations as with (4).
Assuming, as before, that a particular integral of the equation
\[ y^{(m)}=f(x,y,y^{(p)}) \]
will be
\[ y=F(x) \]
and expanding the derivative
\[ y_n^{(p)}=F^{(p)}(x_0+nh) \]
in a Taylor series, we shall have, taking into account the notation introduced by us:
\[ y_n^{(p)}=\frac{1}{h^p}\left[ \Delta_0^p y_0 +\frac{n}{1}\Delta_0^{p+1}y_0 +\frac{n^2}{2!}\Delta_0^{p+2}y_0 + \right. \]
\[ \left. +\frac{n^3}{3!}\Delta_0^{p+1}y_0+\cdots+ \frac{n^{m-p}}{(m-p)!}\Delta_0^m y_0+\cdots+ \frac{n^{m-p+k}}{(m+p+k)!}\Delta_0^{m+k}y_0 \right] \tag{24} \]
Changing \(n\) from \(1\) to \((n+1)\) and taking differences up to order \((m-p)\), we find
\[ \delta^{m-p}y_{n-(m-p-1)} = \frac{1}{h^p}\left[ \frac{A_{m-p}}{(m-p)!}\Delta_0^m y_0 + \frac{A_{m-p+1}}{(m-p+1)!}\Delta_0^{m+1}y_0 + \right. \]
\[ \left. + \frac{A_{m-p+2}}{(m-p+2)!}\Delta_0^{m+2}y_0+\cdots+ \frac{A_{m-p+i}}{(m-p+i)!}\Delta_0^{m+i}y_0+ \right. \]
\[ \left. +\cdots+ \frac{A_{m-p+k}}{(m-p+k)!}\Delta_0^{m+k}y_0 \right] \tag{25} \]
where
\[ A_{m-p+i}=(n+1)^{m-p+i}-(m-p)n^{m-p+i}+ \]
\[ +\frac{(m-p)(m-p-1)}{2!}(n-1)^{m-p+i}+\cdots+ \left[n-(m-1)\right]^{m-p+i}. \]
Transforming the right-hand side of (23) will bring us to equalities (9). We multiply the first of equalities (9) by \(\beta_{0p}\), the second by \(\beta_{ip}\), etc., and substitute into the right-hand side of (23). Replacing the left-hand side of (23) by the right-hand side of (25) and comparing the coefficient at the constants
\[ \Delta_0^i y_0 \]
in the left- and right-hand sides, we obtain the following system, determining the multipliers \(\beta_{ip}\):
\[ \begin{aligned} A_{00}\beta_{0p}&=\frac{A_{m-p}}{(m-p)!h^p}\\ A_{01}\beta_{0p}+A_{11}\beta_{1p}&=\frac{A_{m-p+1}}{(m-p+1)!h^p}\\ \frac{1}{2}\left(A_{02}\beta_{0p}+A_{12}\beta_{1p}+A_{22}\beta_{2p}\right)&=\frac{A_{m-p+2}}{(m-p+2)!h^p}\\ \frac{1}{3}\left(A_{03}\beta_{0p}+A_{13}\beta_{1p}+A_{23}\beta_{2p}+A_{33}\beta_{3p}\right)&=\frac{A_{m-p+3}}{(m-p+3)!h^p}\\ &\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\\ &\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\\ \frac{1}{i!}\left(A_{0i}\beta_{0p}+A_{1i}\beta_{1p}+\cdots\cdots\cdots+A_{ii}\beta_{ip}\right)&=\frac{A_{m-p+i}}{(m-p+i)!h^p}\\ &\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\\ &\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\\ \frac{1}{k!}\left(A_{0k}\beta_{0p}+A_{1k}\beta_{1p}+\cdots\cdots\cdots+A_{ik}\beta_{ip}+\cdots\cdots+A_{kk}\beta_{kp}\right)&=\frac{A_{m-p+k}}{(m-p+k)!h^p} \end{aligned} \tag{26} \]
§ 6. The obtained system (26) differs from system (10) only in that the factor \(\dfrac{1}{h^p}\) appears on the right-hand side. The coefficients of the unknowns \(\beta_{ip}\) remain the same as in (10), and, consequently, the properties of the coefficients of system (10) carry over entirely to (26); however, in order to solve system (26) there is no need for the operations that we performed when solving system (10).
Indeed, suppose we have the equation
\[ v^{(m-p)}=f(x,y) \]
Solving this equation by the same method as equation (1), we must everywhere replace \(m\) by \(m-p\), and the general form of the equation relating the multipliers \(\alpha_{i,m-p}\) will be
\[ \frac{1}{i!}\left(A_{0i}\alpha_{0,m-p}+A_{12}\alpha_{1,m-p}+A_{2i}\alpha_{2,m-p}+\cdots\cdots+A_{ii}\alpha_{i,m-p}\right) =\frac{A_{m-p+i}}{(m-p+i)!} \]
But from (26) we have
\[ \frac{1}{i!}\left(A_{0i}\beta_{0p}+A_{1i}\beta_{1p}+A_{2i}\beta_{2p}+\cdots\cdots+A_{ii}\beta_{ip}\right) =\frac{A_{m-p+i}}{(m-p+i)!h^p} \]
The last two equations will be identities if
\[ \beta_{ip}=\frac{\alpha_{i,m-p}}{h^p},\quad i=0,1,2\ldots k \tag{27} \]
Thus formula (23) takes the form:
\[ \delta^{m-p}v^{(p)}_{\,n-(m-p-1)} =\frac{1}{h^p}\left[ \alpha_{0,m-p}\Delta_0^m y_n +\alpha_{1,m-p}\Delta_1^m y_{n-1} ++-\+ \right. \]
\[ \left. +\cdots+\alpha_{i,m-p}\Delta_i^m y_{n-i} +\cdots+\alpha_{k,m-p}\Delta_k^m y_{n-k} \right] \]
Relation (27), which connects the coefficients of equalities (4) and (23), makes it possible to solve numerically any equation of the form
\[ y^{(m)}=f(x,y,y,\ldots,y^{(m-1)}), \]
without transforming it into a system, which makes it possible to reduce the number of tables to the number of variables entering into the equation. Hence it follows that the equations
\[ y^{(m)}=f(x,y) \]
and
\[ y'=f(x,y) \]
are equivalent in the difficulty of their solution; however, the formulas derived show that increasing the order of the equation does not complicate, but, on the contrary, simplifies the computation. In this respect formula (20) deserves special attention, where two adjacent coefficients are equal to zero.
Without giving the details, we give a summary of formulas for the first four values of \(m\):
\(m=1\)
\[ \delta y_n=\Delta^1_0 y_n+\frac{1}{2}\Delta^1_2 y_{n-2} +\frac{5}{12}\Delta^1_2 y_{n-2} +\frac{3}{8}\Delta^1_3 y_{n-3}+\cdots . \tag{28} \]
\(m=2,\ p=1\)
\[ \delta^2 y_{n-1}=\Delta^2_0 y_n+\frac{1}{12}\Delta^2_2 y_{n-2} +\frac{1}{12}\Delta^2_3 y_{n-3} +\frac{19}{240}\Delta^2_4 y_{n-4}+\cdots \tag{29} \]
\[ \delta y^{(1)}_n=\frac{1}{h}\left[ \Delta^2_0 y_n+\frac{1}{2}\Delta^2_1 y_{n-1} +\frac{5}{12}\Delta^2_2 y_{n-2} +\frac{3}{8}\Delta^2_2 y_{n-3}+\cdots \right]. \]
\(m=3,\ p=1,2\)
\[ \delta^3 y_{n-2}=\Delta^3_0 y_n-\frac{1}{2}\Delta^3_1 y_{n-1} +\frac{1}{240}\Delta^3_4 y_{n-4}+\cdots . \]
\[ \delta^2 y^{(1)}_{n-1}=\frac{1}{h}\left[ \Delta^3_0 y_n+\frac{1}{12}\Delta^3_0 y_{n-2} +\frac{1}{12}\Delta^3_3 y_{n-3}+\cdots \right] \tag{30} \]
\[ \delta y^{(2)}_n=\frac{1}{h^2}\left[ \Delta^3_0 y_n+\frac{1}{2}\Delta^3_1 y_{n-1} +\Delta^3_2 y_{n-2}+\cdots \right]. \]
\(m=4,\ p=1,2,3\)
\[ \delta^4 y_{n-3}=\Delta^4_0 y_n-\Delta^4_1 y_{n-1} +\frac{1}{6}\Delta^4_2 y_{n-2} -\frac{1}{720}\Delta^4_u y_{n-4} \]
\[ \delta^3 y^{(1)}_{n-2}=\frac{1}{h}\left[ \Delta^4_0 y_n-\frac{1}{2}\Delta^4_1 y_{n-1} +\frac{1}{240}\Delta^4_4 y_{n-4}+\cdots \right] \tag{31} \]
\[ \delta^2 y^{(2)}_{n-1}=\frac{1}{h^2}\left[ \Delta^4_0 y_n+\frac{1}{12}\Delta^4_2 y_{n-2}+\cdots \right] \]
\[ \delta y^{(3)}_n=\frac{1}{h^3}\left[ \Delta^4_0 y_n+\frac{1}{2}\Delta^4_1 y_{n-1} +\frac{5}{12}\Delta^4_2 y_{n-2}+\cdots \right]. \]
This summary can be continued as far as desired, using formulas (14), (12), (23), and relation (27).
Formulas (28), (29), (30), and (31) can be obtained from (23), considering that \(p\) varies not from one, but from zero.
§ 7. Examples.
No. 1. Let the equation be given
\[ y'''=y+\sin x \]
with initial conditions \(x_0=0,\ y_0=1.5,\ y'_0=0.5,\ y''_0=0.5\), and it is required to compile a table of values of the particular integral of this equation satisfying the stated conditions, over the interval from \(x=0\) to \(x=1\). The process of solution coincides in all details with the process of solution by the Adams—Störmer method.
We find the first two values of the function by some method, then compile Tables No. 1 and No. 2. Without dwelling on the question of how to find the first values of the function, in the present case, taking the interval \(h=0.1\), from Taylor’s series we find the values
\[ y_1=1.552756 \quad \text{and} \quad y_2=1.612101 \]
corresponding to the values of the argument \(x\):
\[ x_1=0.1,\quad x_2=0.2 \]
We compile the table according to scheme No. 1 (Table 3).
Next, we find the corresponding values of the derivative \(y'''\). These values will be as follows:
\[ y'''_0=1.50000,\quad y'''_1=1.65258,\quad y'''_2=1.81077 \]
We multiply each of them by \(h^3\) and compile the table according to scheme No. 2 (Table 4).
Table 3.
| n | x | y | \(\delta y\) | \(\delta^2 y\) | \(\delta^3 y\) |
|---|---|---|---|---|---|
| 0 | 0.0 | 1.500000 | 0.052756 | 0.006589 | |
| 1 | 0.1 | 1.552756 | 0.059345 | ||
| 2 | 0.2 | 1.612101 | |||
| 3 | 0.3 |
Table 4.
| n | x | \(\Delta^3_0 y\) | \(\Delta^3_1 y\) |
|---|---|---|---|
| 0 | 0.0 | 0.001500 | 0.000153 |
| 1 | 0.1 | 0.001653 | 158 |
| 2 | 0.2 | 0.001811 | |
| 3 |
Table 4 shows that, if in the computation we limit ourselves to five decimal places, the difference \(\Delta^3_1 y\) may be regarded as constant; otherwise we would have had to compute one more value of the function for \(x_3=0.3\).
To find the subsequent values of the function we apply the first of formulas (30), restricting ourselves in it to the first two terms.
For \(n=2\) we shall have
\[ \delta^3 y_0=\Delta^3_0 y_2-\frac{1}{2}\Delta^3_1 y_1 \]
But from Table 4 we have
\[ \Delta^3_0 y_2=0.001811 \quad \text{and} \quad \Delta^3_1 y_1=0.000158, \]
hence
\[ \delta^3 y_0=0.001811-1/2 \cdot 0.000158=0.001732 \]
We enter the obtained value in the corresponding place of Table No. 3 and add it to the number \(\delta^2 y_0=0.006589\), etc., until we reach the column containing \(y\). Having found \(y_3\), we substitute it into the equation and find \(y'''\). We repeat the process.
The results of the computations are presented in Table 5.
Table 5
| n | x | y | \(\delta y\) | \(\delta^2 y\) | \(\delta^3 y\) | \(\Delta_0^3 y\) | \(\Delta_1^3 y\) |
|---|---|---|---|---|---|---|---|
| 0 | 0,0 | 1,500000 | 0,052756 | 0,006589 | 0,001732 | 0,001500 | 0,000153 |
| 1 | 0,1 | 1,552756 | 0,059345 | 0,008321 | 0,001893 | 0,001653 | 0,000158 |
| 2 | 0,2 | 1,612101 | 0,067666 | 0,010214 | 0,002061 | 0,001811 | 0,000164 |
| 3 | 0,3 | 1,679767 | 0,077880 | 0,0012275 | 0,002237 | 0,001975 | 0,000172 |
| 4 | 0,4 | 1,757647 | 0,090155 | 0,0014512 | 0,002422 | 0,002147 | 0,000180 |
| 5 | 0,5 | 1,847802 | 0,104667 | 0,0016934 | 0,002617 | 0,002327 | 0,000190 |
| 6 | 0,6 | 1,952469 | 0,121601 | 0,0019551 | 0,002825 | 0,002517 | 0,000201 |
| 7 | 0,7 | 2,074070 | 0,141152 | 0,0022376 | 0,003047 | 0,002718 | 0,000215 |
| 8 | 0,8 | 2,215222 | 0,163528 | 0,0025423 | 0,002933 | 0,000229 | |
| 9 | 0,9 | 2,378750 | 0,188951 | 0,003162 | |||
| 10 | 1,0 | 2,567701 |
Knowing that the particular integral satisfying the prescribed conditions is
\[ v=e^x+\frac{1}{2}(\operatorname{Cos}x-\operatorname{Sin}x), \]
it is not difficult to verify that the error in the value of \(y\) corresponding to the value \(x=1\) is approximately equal to two units in the sixth decimal place.
No. 2. As a second example, let us consider the equation
\[ y''+\frac{y'}{x}+y=0 \]
with the initial conditions \(x_0=1, y_0=0,7652,\ y'_0=-0,4401\).
Reducing the equation to the form
\[ y''=-\frac{y'}{x}-y, \]
we take the first two values of the function \(y\) and of its derivative \(y'\) from the table given in Zanden’s book,\(^1\) whence the example was borrowed by us. Next we construct a table of differences according to schemes No. 1 and No. 2, as well as a table
\(^1\) (Zanden: Elements of Applied Analysis)
differences of the derivative \(y'\). Having carried out all these preliminary operations, we apply successively formulas (29), in each of which we take only the first two terms. The result of integration over the interval from \(x_1=1\) to \(x_2=2\) with step \(h=0.1\) is presented in Table 6, in the last column of which we have placed the more accurate values of the function, taken by us from Zanden’s book. The discrepancy nowhere exceeds one unit in the fourth digit after the decimal point.
Table 6
| n | x | y | \(\delta y\) | \(\delta^2 y\) | \(y'\) | \(\delta y'\) | \(\Delta_0^2 y\) | \(\Delta_1^2 y\) | More accurate value of \(y\) |
|---|---|---|---|---|---|---|---|---|---|
| 0 | 1.0 | 0.7652 | −0.0456 | −0.0029 | −0.4401 | −0.0308 | −0.0032,5 | 0.0003,5 | 0.7652 |
| 1 | 1.1 | 0.7196 | −0.0485 | −0.0025 | −0.4709 | −0.0174 | −0.0029 | 0.0004 | 0.7196 |
| 2 | 1.2 | 0.6711 | −0.0510 | −0.0022 | −0.4983 | −0.0274 | −0.0025 | 0.0004 | 0.6711 |
| 3 | 1.3 | 0.6201 | −0.0532 | −0.0018 | −0.5220 | −0.0200 | −0.0022 | 0.0004 | 0.6201 |
| 4 | 1.4 | 0.5669 | −0.0550 | −0.0014 | −0.5420 | −0.0160 | −0.0018 | 0.0004 | 0.5669 |
| 5 | 1.5 | 0.5119 | −0.0564 | −0.0010 | −0.5580 | −0.0120 | −0.0014 | 0.0004 | 0.5118 |
| 1 | 1.6 | 0.4555 | −0.0574 | −0.0006 | −0.5700 | −0.0080 | −0.0010 | 0.0004 | 0.4554 |
| 7 | 1.7 | 0.3980 | −0.0580 | −0.0002 | −0.5780 | −0.0040 | −0.0006 | 0.0004 | 0.3980 |
| 8 | 1.8 | 0.3400 | −0.0582 | +0.0002 | −0.5820 | 0 | −0.0002 | 0 | 0.3400 |
| 9 | 1.9 | 0.2818 | −0.0580 | −0.5820 | +0.0020 | +0.0002 | 0.2818 | ||
| 10 | 2.0 | 0.2238 | −0.5800 | 0.2239 |
It should be noted that if we had reduced the given equation to a system of linear equations by the substitution \(y'=z\) and integrated by Adams’ formula, the volume of the entire work would have doubled.
Examples of equations of higher orders will be given together with an investigation of the error of formula (4).