Introducing the integral operator 4.5 (c) 1999 Akio Wakabayashi

4.5 The predictor and corrector.
. . The integral operator obtains the solution of differential equations and the operation is expressed in Eq.(4.35). It must be carried out by use of Eq.(4.31) except for the initial value. If the integrand Y₀' consists of the p values of y₀' to y_p−1', the result ΔY₀ consists of the p+1 values of y₀ to y_p as expressed in Eq.(4.28). Accordingly, the value y_p' at the last point t_p has no effect on the result. By the reason, if the integrand is the function of order p−1, the result comes to the function of order p with no remainder. However, if the integrand is any continuous function, the result with effect of y_p' will have better accuracy than the result without use of y_p', because the integrand is of higher order.
. . The author names the integral operator multiplied by Δ, that is, hD⁻¹ in Eq.(4.28) as the predictor by the reason that it obtains the integral from t₀ to t_p by use of the integrand from t₀ to t_p−1. The usual predictors are obtained by expressing y_p in y₀', y₁', ………, y_p−1' for successive p values as follows.
. . When the matrix D⁻¹ is of order one, that is p=1, it is Euler's rule because the matrix equation is equivalent to Δy₀=h[1]y₀'. When the matrix D⁻¹ is of order two, that is p=2, it is the midpoint rule because of,
Eq4_e1

When the matrix D⁻¹ is of order three, it comes to,
Eq4_e2

When the matrix D⁻¹ is of order four, it is Milne's predictor, that is,

. . Considering the subvector which has not the last component Δ^py₀ of the vector Y₀ integrated by the predictor hD⁻¹, it expresses the integral from t₀ to t_p−1 by use of the integrand of the same range. Hence the author names the operator
hD⁻¹ as the corrector when the matrix D⁻¹ has all zero elements in the last row. Accordingly, the corrector for the predictor of order p is the submatrix of the matrix of order p+1 and is expressed in the rectangular matrix of order (p×(p+1)). The usual correctors are obtained by expressing y_p in y₀', y₁', ………, y_p−1' and y_p' for successive p values as follows.
. . When the matrix D⁻¹ is of order (1×2), it is the corrector for Euler's rule and is the trapezoidal rule because of,
Eq4_e4

When the matrix D⁻¹ is of order (2×3), it is the corrector for the midpoint rule and is Simpson's 1/3 rule because of,
Eq4_e5

Runge-Kutta's method consists of the two predictors, which are Euler's rule and the midpoint rule, and this corrector as shown in Chapter 2 Section 1.4.
. . When the matrix D⁻¹ is of order (3×4), it is Simpson's 3/8 rule, that is,

When the matrix D⁻¹ is of order (4×5), it is Newton-cotes rule using five points, that is,

. . This should be the corrector for Milne's predictor. It is strange idea that the usual predictors should use Simpson's 1/3 rule as the corrector. When we observe microscopic phenomena more precisely, we must increase the magnifying power. In a similar way, we must use the corrector of higher order than the predictor in order to obtain more precise results.
. . Printing out Milne's solutions of the differential equation y' =−ty by every iteration, we shall find that they approach the accurate solutions by once or twice iterations around t=0 but they does not approach the accurate solutions around t=5 even if the corrector is iterated many times. The cause is that it is sufficient to use Simpson's 1/3 rule as corrector around t=0 but it is not sufficient around t=5. It has no effect to use Milne's predictor, because it is sufficient to use the midpoint rule as mentioned in Chapter 2 Section 1 and 2.
. . In Eq.(4.31), the vector ΔY₀ and Y₀' have the equal number of components because the operator is the square matrix of order p. The corrector must append the last component Δ^py₀' to the vector Y₀', because it must obtain the integral from t₀ to t_p by use of the integrand from t₀ to t_p. In the case, the operator becomes the matrix of order p+1 and the vector ΔY₀ has the last component Δ^p+1y₀. It is removed by premultiplying ΔY₀ by Δ^T and the vector Δ^TΔY₀ has the initial value zero as the first element. This operation means subtracting the integral of t to t+h of y' (t) from the integral of t₀ to t+h of y' (t) as expressed in Eq.(4.33) and the matrix equation is expressed in Eq.(4.34). Adding the initial value y₀ to the both sides of Eq.(4.33), the solution y(t) is obtained. It is expressed in Eq.(4.35) and also expressed in following equation.

_______(4.39). . The equation is expressed in detailed matrix equation as follows, denoting the element expressed in Eq.(4.38) by a_n.
Eq4_39a

Premultiplying the both sides by the operator Δ,
Eq4_40

____(4.40). . This matrix multiplied by h is the corrector for the predictor of order p. The matrix is equivalent to the matrix D of order p+1 except the last diagonal is zero. Because all the elements of the last row are zero, it may be expressed in the rectangular matrix of order (p×(p+1)) by omitting the expression of the last row. Eq.(4.40) is equivalent to Eq.(4.39) premultiplied by the operator Δ. Hence,

________(4.41) Accordingly, the corrector which gives ΔY₀ should be denoted by ΔΔ^ThD⁻¹ rigorously speaking but the product ΔΔ^T usually may be omitted because it is equivalent to the unit matrix except the last diagonal element is zero and there is no confusion by omitting it.
. . The integral operator h(ΔD)⁻¹ in Eq.(4.35) is the corrector rigorously speaking, because it obtains the integral from t₀ to t_p by use of the integrand from t₀ to t_p as shown in the detailed expression of Eq.(4.39). However, supposing the last component Δ^py₀' is zero when it is unknown, the operator obtains the integral from t₀ to t_p by use of the integrand from t₀ to t_p−1 and the supposition means the prediction. Hence the author usually uses both of Eq.(4.31) and (4.35) without distinction between the predictor and the corrector.
. . When we will obtain the integral of a given function, it is better to use the corrector of order (p×(p+1)) because it usually gives more precise results than the predictor of order p. When we will obtain the solution of a given differential equation, it is able to use a pair of the predictor and corrector as follows.
. . We first obtain y₀' by the given initial value y₀, which is the only starting value required in case of this calculus. Then we obtain the required precise p solutions by successive use of the predictor of order 1, the corrector of order (1×2), the predictor of order 2, the corrector of order (2×3), ………, the predictor of order p and the corrector of order (p×(p+1)). All the correctors except the last corrector may be omitted because the predictor of order k includes the corrector of order ((k−1)×k).
. . The last corrector must not be repeated many times by use of a constant p, because it has no effect on obtaining more precise result except of twice or three times repeats and has the cause of the unstable result as the usual predictor-corrector methods.

4.6 The double integral operator.
. . The derivative y' (t) is obtained by the differentiation of y(t)−y₀ and the second derivative y"(t) is obtained by the differentiation of y' (t)−y₀'. Hence,

____(4.42) Integrating the both sides from t₀ to t,
Eq4_43

___(4.43) This is the inverse operation of Eq.(4.42) and is unique. If the initial values y₀ and y₀' are given, the function y(t) is the solution of the second differential equation and is expressed in,
Eq4_44

_____(4.44) In order to introduce the double integral operator, Eq.(4.43) is first expressed in the matrix equation using the successive differences of y(t) and y"(t).

. . [Theorem 4.9]. .Eliminating y₀ and y₀' from Eq.(4.43),

_____(4.45) [Proof]. .Denoting the integral variables by x and z,
Eq4_45a

Eliminating y₀ from Eq.(4.43),
Eq4_45b

Then eliminating y₀',
Eq4_45c

[Q.E.D.]

Appendix
Eq4_45d

. . [Theorem 4.10]. .When m≥0, the double integral from t to t+h of Δ^my"(t) is,
Eq4_46

_____(4.46) [Proof]. .It is held for m=0 by Theorem 4.9. Supposing it is held for m=k, it is also held for m=k+1 as follows.
Eq4_46a

[Q.E.D.]

Appendix
Eq4_46b

. . [Theorem 4.11]. .Supposing y(t)∈C^p+1[a, b] and 1≤m≤p+1, there exists a value ξ such that,
Eq4_47

______(4.47) [Proof]. .Δ^my(t)∈C^p+1[a, b] is held because Δ^my(t) consists of the subtractions of y(t+kh). By the mean value theorem,

1). .	if m=1,. .
2)	if m=2,. . where,. .
3)	Supposing the theorem is held for m=n, if m=n+1, . .[Q.E.D.]

. . [Theorem 4.12]. .Supposing y(t)∈C^p+1[a, b], Δ^my(t) is expressed in Eq.(4.48) for the real number r such that t=t₀+rh∈[a, b].
Eq4_48

___(4.48) [Proof]

1)	. .If m=0, Eq.(4.48) is Newton's interpolation formula of y(t) shown in Theorem 3.1.
2)	. .If m=1, it is the difference of the formula as shown in Theorem 3.3. The number of the equidistant points at which the remainder vanishes is decreased to the p points of r=0 to p−1, while Δy(t)∈C^p+1[a, b].
3)	. .If it is held for m=n, by Theorem 4.11 there exist a value ξ such that, This is the remainder of Newton's interpolation formula of p−n degree of the function Δⁿy(t) and vanishes at the p−n+1 points of r=0 to p−n. Hence the formula and shifted formula may be expressed in, These are of the same as the formulas which are obtained by substituting Δⁿy(t) and p−n for y(t) and p of the proof of Theorem 3.3 respectively. . . Considering the function F(u) as the residual of the function Δⁿy(t₀+uh) and its approximate formula with the constant R_q, the function G(u) or F(u+1) is the residual of the function Δⁿy{t₀+(u+1)h} and its approximate formula with the constant R_q instead of R_h. Using that G(u)−F(u) vanishes at the p−n+1 points of u=0 to p−n−1 and r, the following result is obtained in a similar way as the proof of Theorem 3.3. By Theorem 4.11 there exists a value ξ_r such that .. .[Q.E.D.]

. . The step 3 of the above proof may be shown by immediately using Eq.(4.48) substituted the value n for m as follows.
Eq4_48e

where the remainder R is,
Eq4_48f

because there exists a value ξ such that ξ_r<ξ<ξ_r+h by the theorem of intermediate value.
. . The inverse operation of the second differentiation is expressed in Eq.(4.43). When we express it in the matrix equation, we must obtain the successive differences of the left side. If t=t₀, the left side comes to y(t)−y₀−y₀' (t−t₀)=0 but it is not the result of the double integral. If t=t₀+h, it is the result of the double integral but the difference of the two values of the left side has the value independent of the double integral and is not the component of the vector Y₀, because it is expressed in Δy₀−y₀' h. Hence these are not the results of the inverse operation in this operational calculus. The k-th difference of the left side comes to Δ^ky₀, if k≥2. Accordingly, the result of the inverse operation is denoted by Δ²Y₀ and the inverse operation is equivalent to Eq.(4.45).

. . [Theorem 4.13]. .When y"(t)∈C^p−1[a, b], expressing it in Newton's interpolation formula, the successive differences of the results of the double integral is,

___(4.49)
[Proof]. .The relation between y(t) and y"(t) is expressed in Eq.(4.46). Substituting t=t₀ into it, Δ^m+2y₀ is obtained.
Eq4_49a

Expressing the integral variables in x=t₀+rh and z=t₀+qh,
Eq4_50

____(4.50) The integrand Δ^my"(t₀+rh) is expressed in,
Eq4_50a

by changing the notation y and p of Theorem 4.12 for y" and p−2 respectively, because there exists the condition y"(t)∈C^p−1[a, b]. The remainder depends on the (p−1)th derivative of the derivative y"(t), that is, y^(p+1)(t). It is the continuous function of the variable r because the left side of the above equation is continuous. Accordingly, substituting the equation into Eq.(4.50), the proof is completed.

[Q.E.D.]

. . [Theorem 4.14]. .Denoting the double integral in the coefficient of Δ^k+my₀" of Eq.(4.49) by a_k²,
Eq4_51

____(4.51)
[Proof]. .Using Stirling's numbers,
Eq4_51a

When j=0, Stirling's numbers are one for k=0 and zero for the other k values.

[Q.E.D.]

. . Using Stirling's numbers in Table 3.3, the values of Eq.(4.51) are obtained for k=0 to 10 as follows.
Eq4_51b

Hence the relation between Δ²y(t) and y"(t) in Eq.(4.45) is expressed in following matrix equation by arranging the equations in Eq.(4.49) with the m values from zero to p−2.
Eq4_52

__(4.52) where Eq4_52a

. . When the m value is m=p−1, Eq.(4.49) becomes the remainder only.
Eq4_53

_____(4.53) By twice use of the mean value theorem,
Eq4_53a

Hence,________ Eq4_53b

Accordingly, it is not expressed in the matrix equation (4.52).

. . [Theorem 4.15]. .If h→0, the remainder of Eq.(4.52) is h^p+1R→0.
[Proof]. .The diagonal elements h^p+1R_m of the matrix into which the vector h^p+1R is expanded are equivalent to the remainder at t=t₀+mh of the equation which expresses the integrand of Eq.(4.45) in Newton's interpolation formula. Changing the integral variables x, z of Eq.(4.45) into r, q such that x=t₀+rh and z=t₀+qh.

Substituting m into u, Δ²y_m becomes,
Eq4_54

___(4.54) Because of m≤p−2 and r≤q+1≤p, the coefficient of y^(p+1)(ξ_r) is less than or equal to p by Theorem 4.6. Because there exists the value M such that |y^(p+1)(t)|≤M,
Eq4_54a

[Q.E.D.]

. . Accordingly, denoting the matrix of Eq.(4.52) by D⁻², the double integral is expressed in,

_______(4.55) Premultiplying Δ²Y₀ by Δ⁻² mentioned in Theorem 2.12 of Section 2.6,
Eq4_56

___(4.56) This is the operational equation which corresponds to the double integral expressed in Eq.(4.44). Hence the author names the operator h²Δ⁻²D⁻² the double integral operator. This operator requires that the two values of y₀ and Δy₀ are known instead of the two initial values y₀, y₀'. Accordingly, in order to solve the second differential equations by this operator, the initial value y₀' must be transformed to the difference Δy₀.
. . The initial value y₀' is expressed in the first row of the matrix equation of the differentiation Y₀' =h⁻¹DΔY₀, that is,
Eq4_57

_____(4.57) Accordingly, the required value Δy₀ is obtained from this by use of the given initial value y₀' and the components of the vector Δ²Y₀ obtained by Eq.(4.55).
. . The matrix equation (4.55) obtains the results from t₀ to t_p by integrating the integrand from t₀ to t_p−2 as shown in Eq.(4.52). Hence the operator h²D⁻² is the predictor which is the matrix of order p−1 and predicts the two interval from t_p−2 to t_p. On the other hand, the matrix equation (4.56) obtains the results from t₀ to t_p by integrating the integrand from t₀ to t_p, if the integrand may be expressed in the vector with the p+1 components of y₀" to Δ^py₀". The matrix (Δ^T)²D⁻² is the matrix whose first and second rows have all zero and whose third to (p+1)th rows are of the same as the first to (p−1)th rows of the matrix D⁻² of order p+1. Premultiplying the equation (4.56) by Δ²,
Eq4_58

_______(4.58) The operator (ΔΔ^T)²h²D⁻² is the corrector of the double integral and the last two rows are all zero, because the operator (ΔΔ^T)² is equivalent to the unit matrix except the last two diagonal elements are zero. Hence usually omitting the operator (ΔΔ^T)², the corrector of the double integral is expressed in the rectangular matrix of order ((p−1)×(p+1)).

. . [Theorem 4.16]. .The following is always held independently on the remainder.

D⁻²=D⁻¹×D⁻¹_________(4.59) [Proof]. .Considering the function z"(t) which is of (p−2)th degree and equal to the function y"(t) at the p−1 equidistant points of t=t₀+mh for 0≤m≤p−2, where y"(t)∈C^p−1[a, b],
Eq4_59a

The integral and double integral of z"(t) have no remainder because of z^(p+1)(t)≡0. Hence by Eq.(4.32) and Eq.(4.45),
Eq4_59b

Expressing these in the operational equations,
Eq4_59c

The second equation is equivalent to the matrix equation (4.52), because the remainder of Δ²y(t₀+uh) is equal to that of Eq.(4.54) for u=m. Hence,
Eq4_59d

_____________[Q.E.D.]
. . Accordingly, the double integral by h²D⁻² is equivalent to the twice integrals by hD⁻¹. In the latter case, it is able to use the program of the single integral twice and it is not required to transform the initial value y₀' to the difference Δy₀ by Eq.(4.57) because the difference Δy₀ is obtained by the second integral expressed in ΔY₀=hD⁻¹Y₀', using the initial value y₀'. However, it requires the twice calculating times of the former case and there may be the case in which some accuracy of results will be lost by the cause that Theorem 4.16 is not held accurately by the truncation errors of the elements of D⁻¹.

4.7 The n multiple integral operator.
. . The derivative y' (t) is obtained by the differentiation of y(t)−y₀ and the second derivative y"(t) is obtained by the differentiation of y' (t)−y₀'. In a similar way, the n-th derivative is obtained by the differentiation of y⁽ⁿ⁻¹⁾(t)−y₀⁽ⁿ⁻¹⁾. Hence,

___(4.60) Successively integrating the both sides from t₀ to t,
Eq4_60a

____(4.61) This is the inverse operation of Eq.(4.60) and is unique. If the n initial values from y₀ to y₀⁽ⁿ⁻¹⁾ are given, the function y(t) is the solution of the n-th differential equation and is expressed in,
Eq4_62

___(4.62) This is called Taylor's expansion, where the remainder is usually expressed in use of the mean value theorem of the n multiple integral mentioned below. The solution may be also expressed in integrals as follows, transposing the initial values to the right sides of above every integral.
Eq4_63

___(4.63)

. . [Theorem 4.17]. .If y(t)∈C^p[a, b], there exists a value ξ such that,

_____(4.64) [Proof]. .If the value p is one it is held as well known. Supposing it is held for p=n, integrating Eq.(4.63) by use of the mean value theorem of single integral,
Eq4_64a

By the supposition and Eq.(4.62), Eq4_64b

In case of p=n+1, the mean value y⁽ⁿ⁾(ξ) is expressed in,
Eq4_64c

Substituting this into the n multiple integral of Eq.(4.62), it is equivalent to the solution y(t) expressed in the n+1 multiple integral of y⁽ⁿ⁺¹⁾(t). Accordingly, there exists the mean value y⁽ⁿ⁺¹⁾(ξ₁) in case of the n+1 multiple integral of y⁽ⁿ⁺¹⁾(t) and the coefficient φ is n/(n+1) because,
Eq4_64d

Hence,_________ Eq4_64e

. . In order to introduce the n multiple integral operator, Eq.(4.61) is first expressed in the matrix equation using the successive differences of y(t) and y⁽ⁿ⁾(t). Accordingly, the initial values must be all eliminated from the equation because these are not the components of the both vectors of y(t) and y⁽ⁿ⁾(t).

. . [Theorem 4.18]. .Eliminating all the initial values from Eq.(4.61),

________(4.65) [Proof]. .If n=2, it is the case of Theorem 4.9. Supposing it is held in case of n=k, the case of n=k+1 comes to,
Eq4_65a

[Q.E.D.]

. . [Theorem 4.19]. .If the integer m is m≥0, Δ^m+ny(t) is,
Eq4_66

________(4.66) [Proof]. .If m=0, it is the case of Theorem 4.18. Supposing it is held in case of m=k, the case of m=k+1 comes to,
Eq4_66a

[Q.E.D.]

. . [Theorem 4.20]. .When y⁽ⁿ⁾(t)∈C^p−n+1[a, b], expressing it in Newton's interpolation formula with remainder, the following is the successive differences of the n multiple integral from t₀ to t of it, where 0≤m≤p−n.
Eq4_67

____(4.67) [Proof]. .Substituting t₀ into Δ^m+ny(t) of Eq.(4.66), they are Δ^m+ny₀ and expressed in,
Eq4_67a

Changing the integral variable t into the variable r such that t=t₀+rh, the integral range from t₀ to t₀+h is changed into the range of r from 0 to 1 and the integral range from t to t+h is changed into the range of r from r to r+1.
Eq4_68

________(4.68) The integrand is expressed in,
Eq4_68a

by changing the notation y and p of Theorem 4.12 for y⁽ⁿ⁾ and p−n respectively because there exists the condition y⁽ⁿ⁾(t)∈C^p−n+1[a, b]. The remainder depends on the (p−n+1)th derivative of the derivative y⁽ⁿ⁾(t), that is, y^(p+1)(t) and it is a continuous function of r because the left side of the above equation is continuous. Accordingly, substituting the equation into Eq.(4.68), it comes to Eq.(4.67).

[Q.E.D.]

. . [Theorem 4.21]. .When the integers m and n satisfy 0≤m<n,
Eq4_69

_____(4.69) [Proof]. .1). .When m=0, the equation [1] is held because,
Eq4_69a

2). .When m=1, it is also held because,
Eq4_69b

3). .If it is held for m=k, in case of m=k+1,
Eq4_69c

Hence the equation [1] is always held. Denoting the left side of [2] by s_n,
Eq4_69d

. . [Theorem 4.22]. .Denoting the coefficient of Δ^k+my₀⁽ⁿ⁾ of Eq.(4.67) by a_kⁿ,
Eq4_70

_____(4.70) [Proof]. .Expressing the binomial coefficients in use of Stirling's numbers α_j^(k),
Eq4_71

______(4.71) For any integer such as j≥0,
Eq4_71a

Hence, supposing the n multiple integral from r to r+1 of r^j is expressed in,
Eq4_72

_____(4.72) the n+1 multiple integral comes to,
Eq4_72a

Therefore, Eq.(4.72) is always held for the integer n≥1. Accordingly, integrating Eq.(4.71) by the n multiple integral and substituting zero into the integral limit r of the outermost integral, the proof is completed, where a₀ⁿ is one because the integer j only has zero in case of k=0 by Stirling's numbers and the summation of Eq.(4.72) becomes n! in case of r=0 by Theorem 4.21.

[Q.E.D.]

. . By Theorem 4.20 and 4.22, it is able to express the relation between the vectors of y⁽ⁿ⁾(t) and Δⁿy(t) in the following matrix equation, arranging the equations expressed in Eq.(4.67) with the m values such that 0≤m≤p−n.

__(4.73)
where. . Eq4_73a

. . The component Δ^p+1y₀ of the vector of the left side is expressed in the remainder only as follows, because the result of summation of Eq.(4.67) vanishes for m=p−n+1.
Eq4_74

______(4.74) This dose not appear in the matrix equation because it is independent on the components of the vector Y₀⁽ⁿ⁾. There exists a value ξ such that the n multiple integral of the right side comes to y^(p+1)(ξ) by Theorem 4.11 substituted t=t₀. It is also shown by the mean value theorem of the n multiple integral as follows.
Eq4_74a

Accordingly, the value Δ^p+1y₀ vanishes when the functions y(t) and y⁽ⁿ⁾(t) are expressed in the vectors.

. . [Theorem 4.23]. .If both y(t) and g(t) are continuous and g(t) does not change the sign in the interval between t₀ and t, there exists a value ξ such that,

____(4.75) [Proof]. .Supposing g(t)>0 and L≤y(t)≤G,
Eq4_75a

In case of L=G, Eq.(4.75) is evidently held. In case of L≠G, if p=1,
Eq4_75b

Hence,

By the theorem of intermediate value, there exists ξ such that,
Eq4_75d

. .Supposing Eq.(4.75) is held for p=n, in case of p=n+1,
Eq4_75e

Hence,

By the theorem of intermediate value, there exists ξ₁ such that,
Eq4_75g

In case of g(t)<0, the proof is also shown in a similar way.______[Q.E.D.]
. . Denoting the innermost integral of the remainder Δ^mR₀ by f(r), by Theorem 4.17 there exists a value u such that,
Eq4_75h

The function f(u) has an integer j such that it satisfies u<j<u+1 and changes the sign of the binomial coefficient of the integrand. By Theorem 4.23 there exist two values ξ₁, ξ₂ such that,
Eq4_75i

There exist two values u₁, u₂ such that they satisfy u<u₁<j, j<u₂<u+1 and,
Eq4_75j

Accordingly, h^p+1|Δ^mR₀|→0 is satisfied and by Theorem 1.7 the remainder is,
Eq4_75k

. . Therefore denoting the matrix of Eq.(4.73) by D⁻ⁿ, the equation is expressed in,
Eq4_76

______(4.76) This is the operational equation of the n multiple integral expressed in Eq.(4.65). Premultiplying ΔⁿY₀ by Δ⁻ⁿ mentioned in Theorem 2.12 of Section 2.6,
Eq4_77

_____(4.77) This is the operational equation which corresponds to the n multiple integral expressed in Eq.(4.62). Hence the author names the operator hⁿΔ⁻ⁿD⁻ⁿ the n multiple integral operator. This operator requires that the n values from y₀ to Δⁿ⁻¹y₀ are known instead of the n initial values from y₀ to y₀⁽ⁿ⁻¹⁾. Accordingly, in order to solve the n-th differential equations by this operator, the initial values from y₀' to y₀⁽ⁿ⁻¹⁾ must be transformed to the successive differences from Δy₀ to Δⁿ⁻¹y₀.
. . The transformation was the start how the author found this operational calculus. It is expressed in the matrix equation similar to this operational calculus but the matrix is not the operator of this operational calculus because the operand is the vector whose elements are not successive differences but successive derivatives at the initial point. It is the operator interfacing usual numerical calculus with this operational calculus. The details will be mentioned in Chapter 4.
. . The operator hⁿD⁻ⁿ of Eq.(4.76) obtains the n multiple integral from t₀ to t_p by use of the integrand from t₀ to t_p−n because the matrix D⁻ⁿ is of order p−n+1. Hence it is the predictor which predicts the n multiple integral from t_p−n to t_p. When the matrix D⁻ⁿ of Eq.(4.77) is of order p+1, premultiplying Eq.(4.77) by Δⁿ,
Eq4_78

________(4.78) This operator is the corrector of the n multiple integral, because it obtains the n multiple integral from t₀ to t_p by use of the integrand from t₀ to t_p. The operator (ΔΔ^T)ⁿ is equivalent to the unit matrix except the last n diagonal elements are zero. Hence usually omitting the operator, the corrector of the n multiple integral is expressed in the rectangular matrix of order ((p−n+1)×(p+1)).

. . [Theorem 4.24]. .The following is always held independently on the remainder.

_______(4.79) [Proof]. .Considering the function z⁽ⁿ⁾(t) which is of (p−n)th degree and equal to the function y⁽ⁿ⁾(t) at the p−n+1 equidistant points of t=t₀+mh for 0≤m≤p−n, where y⁽ⁿ⁾(t)∈C^p−n+1[a, b],
Eq4_79a

The successive integral of z⁽ⁿ⁾(t) have no remainder because of z^(p+1)(t)≡0. Hence by Theorem 4.18, the n−1 multiple integral of the (n−1)th derivative of z' (t) and the n multiple integral of y⁽ⁿ⁾(t) are expressed in,
Eq4_79b

Calculating the successive differences, these are expressed in the following matrix equations by Theorem 4.19 and 4.20.
Eq4_79d

The second equation is equal to the matrix equation (4.73) together with the remainder. Hence,
Eq4_79e

In a similar way, using Δz⁽ⁿ⁻¹⁾(t₀+rh) and Δⁿy(t₀+rh),
___________ Eq4_79f

________[Q.E.D.]
. . Accordingly, the n multiple integral may be carried out by n times use of the single integral program. In the case, there is no need to transform the initial values from y₀' to y₀⁽ⁿ⁻¹⁾ into the differences from Δy₀ to Δⁿ⁻¹y₀, because the successive integrals are,
Eq4_80

_____(4.80) Hence it is very simple and useful. However, the run time is about n times that of the case using the n multiple integral operator and there may be some cases where the results lose some accuracy by the truncation error of D⁻¹ to the numerical n-th power.

Return to CONTENTS