High-Order Accurate and High-Speed Calculation System of 1D Laplace and Poisson Equations Using the Interpolation Finite Difference Method

Among the methods of the numerical analysis of the physical phenomena of the continuum, the finite difference method (FDM) is the first examined method and has been established as a full numerical calculation system over the regular domain. However, there is a general perception that generality in numerical calculations cannot be expected over complex irregular domains. As using the FDM, the development of computational methods that are applicable over any irregular domain is considered to be a very important contemporary problem. In the FDM, there is a marked characteristic that the theory developed by the (spatial) one-dimensional (1D) problem is naturally applied to the 2D and 3D problems. The calculation method is called the interpolation FDM (IFDM). In this paper, attention is paid to 1D Laplace and Poisson equations, and the whole image of the IFDM using the algebraic polynomial interpolation method (APIM), the IFDM-APIM, is described. Based on the Lagrange interpolation function, the spatial difference schemes from 2nd order to 10th order including odd order are calculated and defined instantaneously over equally/unequally spaced grid points, then, high-order accurate and high-speed computations become possible.Among the methods of the numerical analysis of the physical phenomena of the continuum, the finite difference method (FDM) is the first examined method and has been established as a full numerical calculation system over the regular domain. However, there is a general perception that generality in numerical calculations cannot be expected over complex irregular domains. As using the FDM, the development of computational methods that are applicable over any irregular domain is considered to be a very important contemporary problem. In the FDM, there is a marked characteristic that the theory developed by the (spatial) one-dimensional (1D) problem is naturally applied to the 2D and 3D problems. The calculation method is called the interpolation FDM (IFDM). In this paper, attention is paid to 1D Laplace and Poisson equations, and the whole image of the IFDM using the algebraic polynomial interpolation method (APIM), the IFDM-APIM, is described. Based on the Lagrange interpolation function, the spatial differ...


Introduction
The FDM is a major numerical calculation method mainly used in the field of fluid dynamics. The common recognition of the FDM up to around 2000 was as follows.
"FDM is the oldest method for numerical solutions of partial differential equations (PDEs). On structured grids, FDM is very simple and effective. It is especially easy to obtain higher-order schemes on regular grids. A disadvantage of FDM is that certain conservation laws are not enforced unless special care is taken. Furthermore, the restriction to simple geometries is a significant disadvantage if complex flows are specified" [1]. www.videleaf.com Hoffmann and Chiang [2] issued a comprehensive guide on the FDM and stated the problems when applying the FDM to numerical analysis in arbitrary domains as follows.
"Simple rectangular boundaries rarely occur in nature; indeed, most boundaries are irregular. Irregular boundaries create tremendous difficulties in implementation of the boundary conditions. Various schemes are available to treat irregular boundaries. One may consider using variable step sizes in the neighborhood of the irregular boundaries. This procedure is cumbersome and most cases, very inefficient." [Ref. 2, p.94] Many ideas have been proposed to overcome the aforementioned difficulties with the FDM. The ghost cell method [3], which has been attempted in particular in the field of fluid dynamics, is one noteworthy example. Fukuchi's research [4][5][6][7][8] may be considered an attempt to exceed the conventional limits of the FDM. The IFDM was proposed as a general system subsuming the conventional FDM. Regarding the handling of spatial differences, the IFDM is categorized into the (i) algebraic polynomial interpolation method (APIM) and (ii) boundary polynomial interpolation method (BPIM) [6]. In the IFDM-BPIM calculation system, a method is adopted in which a dummy element is set at the boundary and a dummy value is obtained by successive interpolation in the calculation process so as to match the boundary condition. This method is in common with the idea of the ghost cell method [9][10][11]. From this perspective, the verification of the second-order accurate difference has been investigated with respect to its generality, and the effectiveness of numerical calculations over arbitrary domains has been confirmed [7,8]. There are two cases where all of the Dirichlet boundary conditions are present and the Dirichlet-non-Dirichlet mixed conditions exist as the boundary conditions. In any case of non-Dirichlet conditions in 2D problems, by successively setting the local curved surface at the boundary in the calculation process, the boundary condition comes down to the Dirichlet boundary condition. The mixed boundary conditions method is also incorporated into the IFDM-APIM calculation system, but since the subject is the 1D calculation method, it is not mentioned in this paper. The IFDMwww.videleaf.com BPIM is applied in the case of equidistant grid points, but the IFDM-APIM can be used freely over regular and irregular grid points. Although its basic idea was previously reported [6], ,, uuu u f x y z x y z As ( ) , it reduces to the 3D Laplace equation. (The variable u is indicated in a conventional notation specific to each of the analysis fields, but in this paper, u is uniformly used.) The ultimate goal is to numerically analyze Eq. (1) over arbitrary domains, but in order to do so, we must first consider the 1D problem as comprehensively as possible.
In this paper, we will compare and verify basic conventional FDM numerical analysis methods that numerically analyze 1D Laplace and Poisson equations. The following seven calculation methods are systematically verified: (i) the forward time explicit (FTEX) method, (ii) the time marching successive displacement (TMSD) method, (iii) the successive over-relaxation (SOR) method, (iv) the direct method using the Gauss elimination method (GEM), (v) the direct method using the band diagonal matrix algorithm (BDMA), (vi) the implicit Crank-Nicolson method (CNM), and (vii) the full implicit method (FIM). In section II, focusing on the second-order accurate difference, we comprehensively examine the various calculation methods by IFDM-APIM, and confirm the convergence characteristics. In section III, based on the Lagrange interpolation function, we formulate a method to calculate high-order difference schemes instantly. At the same time, it is confirmed that various calculation methods based on the second-order accurate difference naturally develop to the high-order difference www.videleaf.com calculation methods. In section IV, we numerically analyze four patterns of Poisson's equation using direct method with BDMA and investigate numerical analysis accuracy when using the second-order to tenth-order difference equations.

The Laplace Equation and its Finite Difference Equation
In the paper, the Laplace equation is generally expressed as ( ) , where is the dimension of the domain. As previously described, to progressively understand the simple equation compared to the complicated one, it is better to examine the calculation methods in the order of the 1D Laplace equation ( ) , 2D Laplace equation ( ) , and 3D Laplace equation ( ) . This paper comprehensively considers the numerical calculation method by the IFDM-APIM, limited to 1D problems.
PDEs are classified into (i) parabolic PDEs, (ii) elliptic PDEs, and (iii) hyperbolic PDEs. This paper considers the first two. Conventionally, both numerical methods have each developed as independent numerical calculation methods, but from the perspective that the elliptic PDE is considered the limiting form of the time marching of parabolic PDEs, there is no distinction between them in numerical calculations. Therefore, in this paper, the basic equation for the numerical analysis of 1D Laplace equations is expressed as follows.
Under a fixed boundary condition, as , it becomes ⁄ , then ( ) . The time differential term of the parabolic PDE has a physical meaning. However, the time differential term for solving the elliptic PDE is introduced as a means of numerical analysis. Let the time differential term in this case be defined as the pseudo-acceleration term and call this method the pseudo-acceleration method. 6 www.videleaf.com The method of solving parabolic PDEs is divided into an explicit method and an implicit method. At first, the forward timecentered space (FTCS) scheme expressed by the finite difference equation (FDE) below is adopted as a representation of the explicit method.
The time difference in the FTCS scheme is defined by the following equation [2,5].
In this case, the criteria time difference is equal to the reciprocal of the absolute value of the coefficient of in the FDE of Eq. (3). This is evidence to define in a high-order accurate difference scheme. The stable condition for which the numerical calculation does not diverge is to specify the acceleration factor in the range of . In the case of the second-order accurate difference scheme, the maximum value is theoretically derived by the discrete perturbation stability analysis or the Von Neumann stability analysis [Ref. 2, pp. 113-141].
The variable depends on the time and place. Let us denote this as ( ). The fundamental equation of the 1D thermal conduction problem is expressed by Eq. (2) due to deriving dimensionless equation [12]. In Figure 1(a), the problem of a socalled cooling rod is the problem of setting the temperature of A' and B' to zero, giving the initial condition ( ) and obtaining the time change ( ) of the temperature. This general solution and the solution when the initial condition is in Figure  1(a) are also known. Here, for the sake of explanation, let us treat the problem of the cooling rod as a heating rod problem. A boundary condition of ( ) ( ) is given to the end points of A and B, and an initial condition is provided as 7 www.videleaf.com shown; the initial condition at the center point C is ( ) . Figure 1(a) shows the numerical calculation result using the FTCS. The exact solution to this problem may be expressed as follows using a Fourier series.
where ( | ̅̅̅̅ |) . (We attempted to calculate the case of the initial value of zero between points A and B, but a significant deviation occurs between the numerical solution and the theoretical solution at the initial stage of the calculation. This is thought to be due to the nature of the Fourier series. Therefore, the above initial condition is adopted.) Consideration starts from the case of a regular grid disposition. Figure 1(a) is a calculation example in which the number of divisions is ( ). The numbers of grid points are defined as shown. There is a clear deviation between the numerical solution (dot and interpolation curve) and the theoretical solution (dashed line). The calculation accuracy needs to be evaluated separately for calculation accuracy of the time progression and the calculation accuracy of the convergence value. The error occurring in Figure 1(a) is the error with respect to time progression. However, regarding the convergence value, it is meaningless to compare and verify the accuracy of the spatial difference because the 1D problem here is a problem of the intrinsic-error-zero calculation [6], as described in detail later. This treatment has the advantage that the convergence characteristics of several numerical analysis methods can be clearly compared. Furthermore, if the number of divisions is , the result is strange as shown in Figure 1(b). This is a result of setting the singular point related to the initial condition at the center point C. Usually, since such a singular initial condition is not set and there may be no problem, it shows one characteristic of the numerical calculation by an explicit method. Even if in Figure 1(c), there is a calculation point at the singular point, so unnaturalness remains. However, in the case of Figure 1(d), , the calculation result is sufficiently reliable. By choosing the appropriate 8 www.videleaf.com number of calculation points and time differences, we obtain reasonable calculation results.  . When evaluating precise errors, is adopted. Let us consider the stability of numerical calculation. In Eq. (4), is a reasonable definition in any case? is interpreted as the limit value at which the numerical calculation diverges if it is set larger than the value. This is theoretically determined by slightly disturbing the system in the equilibrium (stationary) state and confirming whether the disturbance decreases; however, it is possible to systematically change the value of to check and determine the maximum value without divergence. The results are shown in Table II, changes depending on whether the number of divisions is an odd number or an even number and converges to 1 as increases. In the method of deriving theoretically [Ref. 2, pp. 113-141], the conclusion comes from the fact that the influence of finite division cannot be evaluated. Therefore, even in the case of adopting a high-order accurate scheme, should be approximately 50 to obtain as a limit value by www.videleaf.com  We can formally develop the discussion into the numerical calculation of the 2D Laplace equation, ( ) . This paper assesses PDEs with higher-order accuracy. Therefore, let us represent the FDE more generally as follows.
In ( ), expresses the order of the algebraic polynomial interpolation, and in the high-order FDE, not necessarily the centered space FDE, m = 3, 5, etc., are also allowed (the meaning is described later). In the case of the quadratic interpolation , the following FTCS scheme is obtained over equidistant grid points. www When is adopted, Eq. (9) is reduced to the scheme of the following 2D Jacobi iteration method [Ref.
It is easy to extend the above definitions to 3D FDE. Equations (6) and (11) guarantee the fastest convergence based on the premise of simultaneous displacement of variables in the calculation process. Let us change Eq. (3) as follows.
This is defined the time marching successive displacement (TMSD) scheme [5]. The definition of the TMSD method is inseparably close to the definition of the FTCS scheme. The concept of the TMSD should be applied to parabolic and elliptic PDEs in general, and it will be able to quickly grasp the convergence value. In this scheme, Eq. (4) is changed as follows.
As ( ) ⁄ , it becomes: This is the 1D Gauss-Seidel iteration method [Ref. 2, pp. 162-164]. The calculation result when is shown in Figure  2(a). In the TMSD calculation, since only the convergence value is concerned, the initial value between points A and B is defined as zero. The convergence characteristics are shown in Table I (II), including the cases of other divisions. The Gauss-Seidel iteration method has meaning as an iteration calculation method of elliptic PDE, but it is a special case of Eq. (12). To execute the calculation of Eq. (12), it is sufficient to successively change the variable in the FTCS calculation. Let us use ( ) ⁄ in Eq. (12), then we obtain the following results.
Replacing the acceleration factor of the TMSD, , with the relaxation parameter, , it becomes the expression of the successive over-relaxation (SOR) method [Ref. 2, pp. 166-167). Therefore, and are completely equivalent. The equivalence of the TMSD calculation and the SOR calculation is limited to calculations based on the second-order accurate difference, which does not hold in the calculation with a higher-order difference, as described in detail later. In the case of the 2D Laplace equation, the optimum relaxation parameter that www Therefore, even when the domain is a cube and the division number is , it also reduces to a one-dimensional expression. The author confirmed that this relation holds in a 4D problem by numerical calculation. From this, the general form of in Eq. (16) is expressed by the following equation.    where is the dimension of the domain. www.videleaf.com

Even if
in a 2D problem, it is convenient if can be defined by reducing to an equivalent 1D calculation. Let us confirm this with Eq. (16-1). It holds that ⁄ ⁄ . When θ is small, the third and subsequent terms on the right side are ignored. Then Eq. (16-1) is expressed as follows.    Therefore, the following equation holds.
is the harmonic mean of and . When and , , then ( ⁄ ) gives the result, . On the other hand, Eq. (16) and (16-1) give the result . The OAF in multidimensional TMSD calculation leads to the equivalent OAF of the 1D TMSD calculation. As will be described later, the concept of the OAF is also effective when using higher-order differences in the TMSD calculation. There, the theoretical solution is unknown and is obtained by numerical calculation. Then it is convenient if it is possible to use the OAF of the 1D problem for the OAF in the 2D and 3D problems. Equation  is the theoretical solution in the case of a rectangular cross-section, but also in the case of an arbitrary section, it can be applied accordingly by adopting the maximum value for the division number of each and .
As in the case of the SOR method, let us define 1 as under-relaxation and ( ) as over-relaxation. There is an OAF, , that guarantees the fastest convergence in the range of over-relaxation. It has been confirmed that the following holds [5].
(i) The convergence process at each calculation point progresses while affecting each other, then the state of convergence progress can be confirmed by specifying one calculation point near the center of the domain. www.videleaf.com (ii) In the process of convergence, there are two cases where the calculated sequence is (I) converging monotonously and (II) converging while oscillating.
(iii) At that time, the OAF, , is a limiting pattern of case (I), at which monotonicity is maintained.
Therefore, the range of the over-relaxation is divided into (I) (the range of monotonous convergence) and (II) (the range of oscillation convergence). Suppose that it is a monotonous increase in (I) the range of monotonous convergence. If it is assumed that , the calculated sequence monotonously increases to the maximum value exceeding the convergence value and follows the process of converging with oscillation. At this time, the oscillation amplitude decreases monotonously. If the convergence range is designated beforehand, there are acceleration factors arriving at this range with the smallest calculation steps. This is the allowable oscillation OAF, expressed as . As the convergence range , it becomes .  Table  I (III) shows the calculation results including cases of other divisions. As the number of divisions increases, the advantage of using increases. In both Table I (I) (FTCS) and Table I (II) (TMSD), the convergence speed will be expressed as . c is a proportional constant inherent to each problem. This relationship always holds if the acceleration factor is a fixed value. However, in the method that adopts , depending on , it becomes . It is ideal to adopt whenever possible.
Figure 2(c) shows the calculation result of the TMSD with . Once the calculated sequence exceeds the convergence value, the oscillation is decreases. If the convergence range is specified as 0.8 to 1.2, it reaches the specified range in three steps, so is one allowable www.videleaf.com oscillation solution of the acceleration factor under this condition. However, ( )=8 and ( )=34 are larger than those in Figure 2 . In this case, the maximum amplitude occurs at ( ), but the definition of the convergence-confirmation grid point defined previously is not changed. The reason is that even if the convergence check is executed at all of the calculation points, there is almost no change in the convergence index evaluation. For example, Figure  2 The difference between a one-point confirmation and an all-point confirmation will not be a problem. If we adopt , it becomes ( )=33, ( )=108. From this, it is understood that even if the approximate value is adopted for the acceleration factor, it does not become a problem.
Currently we deal with the ideal cases where we can use the theoretical solution of the OAF, but strictly the application of the theoretical solution is restricted to the regular domain. Furthermore, there is no known theoretical solution in the higher-order difference FDE. In relation to this, Fukuchi [5] proposed an algorithm to obtain the OAF by numerical calculation and confirmed the agreement with the theoretical solution. In that case, by adopting the "bisection method", it was possible to numerically obtain a strict OAF with the specified number of digits. However, it is cumbersome to set a full-fledged routine. The precision of two decimal places is sufficient. Simple calculation by trial is performed as follows. To obtain an allowable oscillation OAF, we change the previously defined convergence definition . This corresponds as follows. 1 1 Since is different in the calculation points, it is expressed as . Assuming ( ) , the result is consistent with the calculation using the FTCS scheme. This result is generalized in cases of unequally spaced grid points. This method is also effective in the TMSD scheme. This conclusion is obtained as a result of imparting the pseudoacceleration term to an elliptic PDE. The above theory evolves into the derivation method of higher-order difference scheme. The polynomial interpolation is not limited to an algebraic polynomial. Any interpolation equation in which the derivative is defined could be used.
Let us refer to the method of directly obtaining the convergence value. In Figure 1, the following relationship holds at 4 calculation points, .
Depending on the Dirichlet boundary condition, and the solution is . If the number of calculation points is small, the solution can be found instantly using the Gauss elimination method (GEM). However, as the number of calculation points increases, the computation time increases substantially. Equation (18) is a tridiagonal matrix, and when using the tridiagonal matrix algorithm (TDMA) [Ref. 13, pp. 50-54], the calculation time is shortened. In this paper, the expression for using the TDMA is as follows.
Matrix is obtained by placing the diagonal elements of matrix in the center column of matrix B. x is undefined and irrelevant to the calculation. Using the TDMA, the calculation is completed by the operation of matrix , and the storage capacity is saved. www.videleaf.com When the number of elements is small, there is little difference between the computation time of Eq. (18) using the GEM and that of Eq. (19) using the TDMA; however, as the number of elements increases, a decisive difference occurs in the computation time. Table III shows the computation time CT(s) using the author's PC. Although there are some shakes in these computation times, they are effective as reference values. The calculation results using the TMSD with the OAF (TMSD-OPT) are also shown regarding . The convergence time with ( ), 2 holds. In the GEM calculation, the evaluation of the calculation time is , which agrees with our knowledge. However, in the TDMA calculation, the calculation time is 0.5× seconds or less with . The calculation results of and using the TDMA are shown. The estimation of the calculation time is 1 as expected. (Such an examination is unnecessary as far as it is limited to a 1D analysis, but it is not special in a 3D analysis that the number of elements exceeds 10000.) In the GEM calculations, we found that a large calculation time is required as the number of elements increases. In 2D and 3D problems, Eq. (18) becomes a diagonal dominant large sparse matrix. Then we have to consider an effective calculation method to replace the GEM. The conjugate gradient method (CGM) is one method of calculating this problem. The convergence speed using this method is almost equivalent to the TMSD [6]. However, since it is limited to the case where the coefficient matrix is symmetrical, it cannot be applied to higher-order differences. SOR is also commonly used. When the diagonal components are equal, let be the expression obtained by transforming u=f such that the diagonal component is 1.
Rearranging this, the following equation is obtained.
The term for setting the boundary condition shifts to the right side. For each calculation step, simultaneous equations are solved using the TDMA. With reference to Eq. (19), the expression defining the whole of Eq. (20) is expressed as . is produced at the beginning of the calculation and is used in common in each calculation step. When β = 1, it becomes the full implicit method, and for β = 0.5, the Crank-Nicolson method.   Figure 1 has been resolved. The implicit method is also very useful when only the convergence value is a problem, then we use full implicit method, β = 1. In this method, no matter how much is taken, the calculated sequence does not oscillate and uniformly converges. In other words, the Crank-Nicolson method is conditionally stable, whereas the full implicit method is unconditionally stable [Ref.
2, pp. 66-68]. Figure 3(b) shows the calculation result of the full implicit solution with . There is a discrepancy between the numerical solution and the analytical solution, but this is not a problem because only the convergence value is concerned. If we use , it becomes ( )=1, ( )=1. The full implicit scheme is expressed as follows with β = 1 in Eq. (21).
Dividing both sides of this expression by and setting results in the direct method as shown in Eq. (18). Therefore, the full implicit solution has significance as a semi-direct method in 1D problems. We studied the case of the "regular" domain. The physical phenomenon is irrelevant to the coordinate system. When this is applied to a 1D analysis, the same result (with little difference) must be obtained using an arbitrary coordinate movement. Figure 4(a) shows the result obtained by shifting the coordinates in Figure 1(a) to the left by . It is defined as a regular element disposition [4]. The reason to investigate this problem is that it is necessary to consider regular element disposition for examining arbitrary domains in 2D and 3D problems. Although the calculation result is strange, it is the same as the phenomenon occurring in Figure 1(b), and the influence of a singular point in the initial condition setting appears. Here, only the convergence of the numerical calculation is a problem. The same problem is calculated using in Figure 4(b). In Figure 4(a), the numerical calculation usually www.videleaf.com diverges using . However, it does not diverge because we changed Eq. (7) as follows [5].   [14]. With reference to Figure 4(b), it is defined as follows.
As , it results in a usual second-order accurate central-difference scheme and becomes ( ) ⁄ . The criteria time difference is ( ) ⁄ ⁄ . As described previously, using Eq. (7), the calculation diverges when ( ) ⁄ is adopted in Figure 4(a). At , When adopting the reciprocal number for the criteria time difference, it becomes ( ) ⁄ and gives a stable convergence. However, assuming that takes an arbitrary value in the range of , a problem arises. If it is assumed to be , it becomes ( ) ⁄ , and we have to use a very small time interval. The first way to avoid this problem is to adopt . This process is adopted in the ghost fluid method 9,10 . As a second method, it is possible to avoid the problem by adjusting the time difference in this part to ( ⁄ ). Equation (23) means that the latter method was adopted. To eliminate the influence of the initial condition singularity in Figure 4(a), we adopt as shown in Figure 4(c). The numerical solution and the theoretical solution are in good agreement. Using 0.8 in Figure 4(c), in the movement ( ) of the coordinate, regardless of how is specified, the numerical solution and the theoretical solution seem to be perfectly matched only by changing the calculation points. The definition of Eq. (23) seems to be inevitable that it is not simply a way of suppressing divergence but must be specified like this. As previously mentioned, if we choose the grid disposition with the regulation or adopt the Crank-Nicolson method, this problem can be completely avoided. In the TMSD calculation, only the convergence value is a problem, so the above method is www.videleaf.com completely effective. Figure 4(d) is a calculation example corresponding to Figure 2(b). The OAF, , of can be applied accordingly. (Calculating using a numerical stability analysis will result in different results, but not very different.) In this case, the calculation using Eq. (23) is equivalent to the calculation by SOR. This can be easily confirmed by numerical calculations.

Derivation of High-Order Finite Differences Based on Lagrange Interpolation Function
Conventionally, the difference scheme was derived via Taylor series expansion [15]. Let us assess the second-order accurate scheme and the SW scheme. The following two formulas are obtained using the Taylor series expansion.
This equation reduces to Eq. (27) under the condition of , but when , the fourth term on the righthand side remains. If it is defined that the highest-order error term shows the accuracy of the difference, this is the first-order accuracy. Equation (27) has been conventionally referred to as a second-order accurate difference, which is a special form of Eq. (30). Both are inseparable from quadratic interpolation curves as described below. For this reason, we use the term "second-order difference" to include both. Although the above method has the advantage that it is possible to evaluate the accuracy of the difference, it becomes difficult as it gradually shifts to a higherorder difference and handling irregularly spaced grid points. At this time, from the perspective that it is sufficient to calculate only the difference scheme, more efficient derivation methods are necessary. Consider the following Lagrange interpolation function passing through three points. 2 3 1  3  1  2  2  1  2  3  1  2  1  3  2  1  2  3  3  1  3 The notation of ( ) indicates that this is a quadratic expression. It is obvious that this is equivalent to a quadratic equation going through the three points ( ) ( ) ( ). The first term on the right side is at and zero at and . The second and third terms can be explained in the same way. Let us express each term as follows.
The same result can be obtained by arbitrarily setting the coordinate zero point, but in this paper, it is defined as . The first derivative changes depending on the position. The second derivative is constant regardless of the position. This is a characteristic of a quadratic curve. Equation (34) is equivalent to the SW scheme. Figure 5 shows the calculation result of the Lagrange interpolation function calculated using the difference coefficient matrix, defined as and . The above method is expanded to the interpolation expression of an arbitrary order, and instantaneously the derivatives, of which the orders are equal to or less than the order of interpolation www.videleaf.com equation, are derived at the designated position. This process can be determined instantaneously using a calculator. As can be seen from the derivation process, this method is applied irrespective of the calculation point disposition (equidistant or not). Using Lagrange's polynomial interpolation directly, Figure 5 can be drawn, but it is difficult to define the difference schemes. Numerically, the same results can be obtained from n-th order polynomials; however, it is impossible to define the difference schemes. The advantage of obtaining is that the difference schemes are instantaneously determined, and the explicit/implicit methods become possible using high-order difference schemes. As previously mentioned, the method of deriving the difference scheme using the Taylor series expansion can evaluate the accuracy, but it lacks maneuverability. Although it is not possible to directly evaluate the accuracy of the difference scheme in the above method, it is possible to evaluate the accuracy by comparing the results of numerical calculations with some analytical solutions [6]. The above method was published as the author's original method; however, Klein and Berrut [16] previously suggested a method to derive finite difference formulas over an equidistant grid interval using Lagrange's polynomial interpolants. Although there seems to be a difference in the algorithms between the author's method and their method, it was confirmed that their results of finite difference formulas can be obtained from the author's method as well. www.videleaf.com  Table IV shows the results of the central difference of at equidistant intervals ( ). In the central difference in Table IV, the center point number is expressed as zero.   in this table are shown in finite decimal numbers, but they are instantaneously and exactly given in the process of numerical calculation analysis).
As previously mentioned, in this paper, a local coordinate system with at the left end point is used. For this reason, the second-order (accurate) difference is expressed as shown in Table V-P2. For second-order differences, the interpolation function is quadratic and the evaluation of is only the center point, but in fourth-order differences, the evaluation of there are three points including the central difference, except for two end points, and it changes depending on the position. They are expressed as: Including the other orders, the general expression is as follows.
Now it is assumed that =constant in Eq. (36). Then ( ) , meaning ∑ . Table V shows the calculation results. In the central difference, in Eq. (36), ⁄ ( ) is the central term. Table V demonstrates that the maximum value of the absolute value of each coefficient is | |. In second-order differences, the criteria time difference is set to ( ) ⁄ in Eqs. (4) and (13). There may be instances that commonly use this standard. However, in this paper, it is used as follows. It is equivalent to define ( ) ⁄ , and in ( ) to define | |. Therefore, the criteria time difference is generally defined as | |. The m-th order difference scheme is not limited to the central difference. The expression of Eq. (36) holds even if m = 3, 5, , 9. Table VI shows the calculation results. For odd-order differences, in the case of equidistant grid points, there are two occurrences of the maximum value of the coefficient (usually one place in case of unequal spacing). In odd-order differences, the criteria time difference is determined in the same manner as in even-order differences. Table VI-P3 demonstrates the third-order difference retracted to the secondorder difference. That is, at , ( ) is determined irrespective of .
This is the case of equidistant grid points, not at unequally spaced grid points. The validity of this result can be easily confirmed by numerical calculations and is also derived from Lagrange interpolation function. Tables VI-P5, P7, and P9 demonstrate that the retraction occurred in the pair of central portions. www.videleaf.com Figure 6(a) shows the result calculated by the forward time explicit method using the fifth-order difference scheme, ( ) , called FTEX(5). From Table VI-P5, it becomes ( ) ⁄ , and is used. The calculation accuracy with respect to time progression is equivalent to the case in Figure 1. Calculating this using the Crank-Nicolson implicit method improves the calculation accuracy even if the acceleration factor . Figure 6(b) shows the result calculated with β = 0.8, which is highly accurate, but the rationale for specifying β = 0.8 is unclear. When adopting the acceleration factor in the full implicit method (β = 1), the calculation becomes equivalent to the direct method and converges in one set calculation. Similar results are obtained for the seventh-order difference and the ninth-order difference. The TMSD calculation can also be applied and convergence accelerates. An important fact arises here. For the second-order accurate difference, the TMSD calculation and the SOR calculation is completely equivalent. In this case, however, they deviate. Figure 6(c) shows the calculation result via the TMSD. The convergence indexes are ( ) . Figure 6(d) shows the calculation result by the SOR. The convergence indexes are ( ) 42. The TMSD www.videleaf.com calculation converges faster than the SOR calculation. In the actual calculation, we do not use an odd-order difference alone (described in detail later). However, in the multidimensional problem and the higher-order difference numerical analysis, it is possible to improve the calculation accuracy by partially using the odd-order difference at the corner portion of the irregular domain.  Figure 1(c) using the fourth-order difference. Figure 7(a) adopts three fourth-order polynomials , and . This is a sparse setting of polynomials [6]. For example, the interpolation of is used for the calculation at each point of and . Equations (35-1), (35-2), and (35-3) are used at each calculation point, respectively. In this setting, if the number of divisions is even, the polynomial setting is symmetrical, but if it is an odd number, symmetry cannot be always maintained. A total of 7 polynomials are used in Figure 7(b). This is a dense www.videleaf.com setting of polynomials [6]. In this setting, the symmetry is maintained irrespective of the even number and odd number of divisions. This is called the standard fourth-order API scheme and defined as SAPI (4). The setting methods of polynomials at are innumerable. It is also possible to set even-order and odd-order mixed polynomials and to set only odd-order polynomials. The corresponding numerical solution is obtained if the polynomials cross one another. However, if ( ) , the symmetrical polynomial setting is uniquely determined only due to specifying the value of m. The highestorder polynomial considered in the paper is up to the tenth order. In Figure 7(b), in ( ) 6, the calculation results seem to visually coincide strictly with analytical solutions. A slight deviation occurs at m = 8, and at 10, it is reduced to as described previously. Its time progressing solution is not very accurate.  (6). The TMSD calculation is superior to the "SUR" calculation; in particular, in the SUR calculation in Figure 7(f), divergence occurs unless it is set to under-relaxation ( ), and the convergence speed decelerates. When using the direct method, the symmetry of the matrix confirmed by the secondorder difference is not retained. The calculation equation for , SAPI (4)  , which includes the TDMA as a special case, the computation time is terminated at 7/3 of the calculation in the case of the second-order difference. When there are many calculation points, a decisive difference occurs in the computation time for both. Also, in the calculation of the implicit method, a 100 × 7 band matrix must be adopted. The bandwidth of ( ) is generally 2 m-1 and the number of the centre position is m when m is an even number. The above shows the handling of higher-order differences in regular domains, but it can manage exactly the same in irregular domains. In the calculation of FTEX and the TMSD and in the calculation of both calculation end points over an irregular domain, time difference width adjustment is performed as shown in Eq. (23). In the implicit solution, such an adjustment is unnecessary.
Finally, we consider the problems of unequally spaced grid points that will be needed in the future. Although the discussion thus far can be applied to any unequally spaced grid point, singularities must be considered in the FDM. Singularities in the initial condition were examined in this paper, but singularities in the boundary shape may be a problem. The need to set the difference interval smaller for a specific area occurs in multidimensional problems. In preparation for such a situation, we will confirm the problem of unequally spaced grid points in www.videleaf.com 1D problems. Figure 8(a) is an example calculated using FTEX with ( ) and ⁄ in the range using ( ). The criteria time difference is defined from and the central difference in ( ), becoming , and as , it becomes . Changing the subdivision range to 0.3 ≤ x ≤ 0.6 does not cause any problems. The result is shown as The upper limit of where the calculation does not diverge is . The calculation can be done using all ( ) , but changes individually. In this case, using the implicit Crank-Nicolson method (CNM), the range of the subdivision point is set to 0.3 ≤ x ≤ 0.6, and the result is calculated by ( ) as shown in Figure 8(b). Compared with Figure 8(a), the number of calculations is 1/20. Using in the full implicit method, the convergence value can be obtained by one set of calculations. Figure 8(c) shows the TMSD calculation result using SAPI (6); the convergence index is , NC(3)=1910. If we use , the calculation diverges; therefore, it is . This is characteristic of the TMSD calculation using high-order differences. Figure 8(d) demonstrates the SUR calculation result using SAPI(6); convergence index is , ( )=2559. In this case, if we use , the calculation diverges. The TMSD calculation is superior to the SUR calculation.  (6) , ( )=2559.

Numerical Calculation Accuracy Regarding 1D Poisson Equations
The convergence of the 1D problem has been examined thus far, but here we assess the accuracy of the convergence value. Since the convergence process is not a problem, the numerical calculation is performed via the direct method with the BDMA. The Laplace equation ( ) has only a few variations, therefore, the Poisson equation ( ) is used. All of the calculation methods studied in the Laplace equation are also effective in the Poisson equation. We will examine SAPI(m), but we cannot visually distinguish the difference in the results using . However, there is a large difference in the calculation accuracy. Figure 9(I) demonstrates a case of that results in the previously discussed problem. When there is almost no digit-loss error, m = 2, the average of the absolute value of the error (the absolute error) is the smallest ( ). The results including other cases are shown in Table VII(I). In Table  VII, the area surrounded by the thick frame is the part of the intrinsic-error-zero range. Theoretically, all are the intrinsicerror-zero calculations; however, we can clearly confirm the digit loss error in the higher-order difference calculations. As a whole, the lower-order difference and the smaller number of divisions ensure the calculation accuracy. Fukuchi [6] confirmed that high-order differences do not necessarily guarantee high accuracy. First, the problem of Runge's phenomenon [16,17] occurring using polynomial interpolants can be considered a reason. This phenomenon will inevitably cause an error in the evaluation of the derivative coefficient. Second, when using higher-order differences, the digit-loss error becomes noticeable. The numerical calculation is performed in double precision, but to improve digit-loss error of the higher-order difference, it is necessary to calculate in quadruple precision. www.videleaf.com  Figure 9(II) shows a case where the analytical solution is . Let us first consider the case where is a positive integer. The calculation example is the case of . In SAPI(m), the minimum value of m is 2. If , it becomes the intrinsicerror-zero calculation for all SAPI(m). However, if and m = 2, it is the intrinsic-error-zero calculation. This is because in the equally spaced central difference equation, Eq. (27), there is no third-order difference term in the error terms. The calculation results are shown in Table VII(II) for k=3. There is an error arising from digit-loss error of higher-order differences, so the overall trend is similar to Table VII(I). In Table VII (II) for k = 4, m ≥ 4 gives the intrinsic-error-zero calculation. In the calculation of SAPI(2), it can be confirmed that it is the secondorder accuracy (as the difference width is 1/10, and the error becomes 1/100), whereas, in the range m ≥ 4, the digit-loss error becomes dominant. In Table VII(II) for k = 5, the range of www.videleaf.com intrinsic-error-zero is , because we have to use differences other than the central difference. If the analytical solution is an n-th algebraic polynomial, the Taylor series expansion has finite terms. In this case, if it is at least m ≥ n in SAPI(m), it becomes the intrinsic-error-zero calculation. However, in all other cases, intrinsic errors always occur. The case where is shown in Table VII(II)k = 2.5. Figure 9(III) and Table VII(III) for k = 4 show a case where the analytical solution is ( ). The error criterion in this case is ( ). Figure 9(IV) and Table VII(IV) for k = 8 demonstrate a case where the analytical solution is ( ). The calculation error of SAPI(m) exhibits a complicated appearance based on the combination of , and digit loss errors.
Below we describe the basic features of error under double precision calculation.
(i) As 10 and m = 2, for second-order differences, the error is less than 1%. In this case, the digit-loss error is minimized, and the error shrinks in proportion to ⁄ .
(ii) Although SAPI (8) and SAPI (10) are effective when the division number is small, as the number of divisions increases, the precision tends to decrease due to digit-loss errors. It should not be used with 100.
(iii) In conclusion, SAPI(4) or SAPI(6) can be used for higherorder differences, and by setting 1000, the calculation accuracy ( ) or less can be expected. This means that the calculation error will be within . In this case, the error does not result in shrinking in proportion to ⁄ . The influence of the digit-loss error and Runge's phenomenon is involved in a complicated way.

Conclusion and Discussion
Starting from the usual second-order accurate difference, there are two methods to realize high-precision calculation with FDM. One is called compact FDM, and there are many reports on its usage [18,19]. The other is the method described in this paper. www.videleaf.com Both are important, but the latter has not yet been systematically studied. The reason is that there is no simple algorithm for deriving higher-order differences and managing the boundary is difficult. This paper shows the possibility to fully solve these problems. The contents of this paper are summarized as follows.
(1) The solution to calculating higher-order differences is the difference coefficient matrix defined in this paper, and it is obtained from Vandermonde matrix, as . From this, difference schemes of each class can be obtained instantaneously. The results calculated using this method are completely equivalent to the conventional results in equally spaced central differences. It is also applied at unequally spaced grid points and has meaning as a generalized scheme of the SW scheme.
(2) The solution to elliptic PDEs is a convergent solution when the boundary is fixed in parabolic PDEs. From this perspective, there is no need to distinguish between their numerical calculation methods.
(3) In the second-order difference, the simplest scheme of parabolic PDEs in explicit methods is the FTCS scheme, assuming the simultaneous displacement of the variables. When these are successively displaced, it becomes a calculation system including all of the conventional iterative calculation methods of elliptic PDEs. This is called the TMSD scheme. This method evolves into a numerical analysis method of 2D and 3D problems.
(4) In second-order accurate FDM calculations, calculations using the TMSD and SOR are completely equivalent. However, they are not equivalent in calculations using a higher-order difference. Thus, it was confirmed that the TMSD calculation is faster in convergence.
(5) Although there is a direct method for numerically solving elliptic PDEs, the calculation using the GEM requires much time. When using the TDMA, the computation time is surprisingly shortened. www.videleaf.com Then, the numerical calculation of elliptic PDEs can be performed using the full implicit method. By increasing the time difference (acceleration factor), it always converges with several iterations (usually less than 10 times), and high-accurate and high-speed calculations over arbitrary domains become possible as in the 1D problem. The details will be reported separately in the near future.