# Numerical Method for Solving Nonhomogeneous Backward Heat Conduction Problem.

1. Introduction

The heat conduction equation is a kind of very important time-dependent parabolic partial differential equation; it describes the distribution of heat or temperature in a given region over time and is widely used in diverse scientific fields, such as the study of Brownian motion [1], to solve the Black-Scholes partial differential equation [2] and the research of chemical diffusion. Many works have been done to study the heat conduction problem [3-5]. The purpose of this article is to numerically solve nonhomogeneous backward heat conduction problem. This problem is one kind of inverse problem, also called final value problem or time inverse problem.

In many engineering and application areas we need to reconstruct the unknown initial heat energy or temperature from the final measured one; it is a typical inverse problem which is related to the initial boundary value problems in heat conduction. As we all know, there are many kinds of inverse problems in mathematical physics, such as boundary inverse problems [6, 7], coefficient inverse problems [8, 9], and evolutionary inverse problems (or time inverse problem) [10, 11]. The backward heat conduction problem is a time inverse problem, in which the initial conditions are unknown; instead the final data are observable.

In the sense of Hadamard [12], the inverse problem is always considered as a class of classically ill-posed problem, which means that a small error in the input data may cause enormous error to the result. The backward heat conduction problem is a typically ill-posed problem, because of the measurement error of the input data. Most standard numerical method can not achieve a good result in solving an ill-posed problem, due to the ill-posedness of the problem and the ill-conditioning of the discretized matrix, so, it is very important to find a stable and efficient numerical algorithm to solve ill-posed problems.

The backward heat conduction problem arises in many physical and chemical fields and has attracted the attention of many researchers. The improperly posed backward heat conduction problem has been considered by L. Payne [13], and Miranker [14] has given a research of the uniqueness conditions for the backward heat conduction problem. There are also many works that have been done to numerically solve the ill-posed backward heat conduction problem. In [15], Fourier regularization method is used to solve one-dimensional backward heat problem, H. Han [16] has considered this problem using boundary element method, and also some other numerical methods for solving backward heat conduction problem have been given in many works, such as finite difference method [17, 18], iterative boundary element method [19, 20], the method of fundamental solution [21], etc. Most of these methods are used to solve one-dimensional homogeneous situations; there are a few papers on the nonhomogeneous case in higher dimensional space. Scientists M. Denche and A. Abdessemed [22] gave extensions of the quasi-boundary methods to the nonhomogeneous case. Paper [23] regularized the two-dimensional nonhomogeneous backward heat problem by perturbing the final value. M. Li, T. S. Jiang, and Y. C. Hon have solved nonhomogeneous situations by using a new meshless method based on radial basis functions [11].

In this paper we consider the nonhomogeneous backward heat conduction problem in two- or three-dimensional spatial domains. We propose a new numerical method to solve this severely ill-posed problem; coupled with the likewise Crank-Nicolson scheme and a new variable, the backward problem is transformed to a nonhomogeneous Helmholtz type problem; by solving this nonhomogeneous Helmholtz type problem and submitting back to the new variable, we can get the solution of the initial condition.

The structure of the paper is organized as follows: In Section 2, we briefly introduce the formulation of the nonhomogeneous backward heat conduction problem. In Section 3, we introduce the numerical method and apply this method on the inverse problem. The results of several two- and three-dimensional numerical experiments are presented to illustrate the stability and accuracy of the proposed method in Section 4. Section 5 is dedicated to a brief conclusion. Finally, some references are introduced at the end.

2. The Formulation of the Problem

We consider the nonhomogeneous backward heat conduction problem in a bounded and connected domain [bar.[OMEGA]] [subset] [R.sup.d] (d = 2,3) with the boundary r = [partial derivative][OMEGA] in the following,

[partial derivative]/[partial derivative]t u(x, t) = [DELTA]u(x,t) + f(x, t), x [member of] [bar.[OMEGA]], t [member of] (0, T), (1)

with the final temperature condition,

u (x, T) = g(x), x [member of] [OMEGA], (2)

and the Dirichlet boundary conditions,

u (x, t) = h(x, t), x [member of] [GAMMA], t [member of] (0, T), (3)

where f(x, t), g(x), and h(x, t) are known functions, [DELTA] is Laplace Operator, and u(x, t) satisfies the nonhomogeneous heat conduction equation (1). The problem which we want to solve in this paper is the backward heat conduction problem; we will determine the temperature before some particular time T from the known data (2) at time T and the boundary conditions (3).

3. The Computational Algorithm

To solve the backward problem (1)-(3) numerically, we introduce the grid of time,

[[omega].sup.t] = {[t.sup.n] | [t.sup.n] = [eta], n = 0, 1, ..., N, N[tau] = T], (4)

Using the notation [u.sup.n](x) = u(x, [t.sup.n]) and a similar format as Crank-Nicolson (C-N) scheme to approximate the heat equation (1) with second order at the time moment [t.sup.n+1/2] = ([t.sup.n+1] + [t.sup.n])/2, we have

[u.sup.n+1] - [u.sup.n]/[tau] = 1/2([DELTA][u.sup.n+1] + [DELTA][u.sup.n]) + f(x, [u.sup.n+1/2]), (5)

where n = 0,1, ..., N - 1. Here, it is worth mentioning that in (5) we did not use the C-N approximation scheme, but a scheme has a similar format as C-N scheme, which we called the likewise C-N scheme [24].

By using a new intermediate variable w, which is given as follows,

[??] = 1/2 ([u.sup.n+1] + [u.sup.n], (6)

(5) can be transformed in the following new form,

[mathematical expression not reproducible]. (7)

By simplifying (7), we can get

[mathematical expression not reproducible]. (8)

Notice that (8) is a nonhomogeneous Helmholtz equation with m as solution; it can be written as

[mathematical expression not reproducible], (9)

where

[k.sup.2] = 2/[tau],

F(x, t) = 2/[tau] [u.sup.n+1] - f(x, [t.sup.n+1/2]), (10)

with the boundary conditions,

[??] = 1/2 ([h.sup.n+1](x) + [h.sup.n](x)), x [member of] [GAMMA], (11)

where [h.sup.n](x) = h(x, [t.sup.n]) is known from the Dirichlet boundary conditions (3).

To get the solution u from the nonhomogeneous Helmholtz equation (9) with boundary conditions (11), the finite element method is used. We multiply by a test function v, which vanishes on the boundary, and apply integration by parts, getting the variational problem with n= 0,1, ...,N -1 as follows,

[mathematical expression not reproducible] (12)

where [k.sup.2], F(x, t) are given in (10). We can obtain the value of [??] by solving the nonhomogeneous Helmholtz equation (8) (or (9)); then, the solution [u.sup.n] can be computed from (6). After repeating the above steps N times, we can get the initial temperature at n =0.

The Helmholtz equation often appears in the study of physical problems; it is caused by the propagation of time harmonics and applied in many science and engineering fields, such as acoustic, electromagnetic science, and geophysical problems. Many experts have studied solving Helmholtz equation by finite element method (FEM); F. Ihlenburg and I. Babulka studied the finite element solution of the Helmholtz equation with high wave number by using the h-version [25] and h-p version [26] of the FEM. The Least-Squares stabilization of finite element computation for the Helmholtz equation is considered by I. Harari and F. Magoules [27]. Y. Wong and G. Li [28] consider a new finite difference scheme for solving the Helmholtz equation at any wavenumber, etc.

From [25-28] and the references it can be found that, due to the so-called pollution effect [29], to compute approximate solutions of the Helmholtz equation for high wave numbers by the standard finite element method (FEM) is unreliable, so, in (9), the wave number should be not high. Suppose the numerical discrete grid size is h; in order to ensure an accurate numerical solution, the condition [k.sup.2]h < 1 must be enforced and it is necessary to require kh to be small. A conclusion that the stability constant does not depend on k if [k.sup.2]h is bounded using FEM to solve the Helmholtz equation was given in [25, 26]. Also, a convergence theorem is stated in [30] under the assumption that [k.sup.2]h is sufficiently small. As a result of the fact that the wave number k in the Helmholtz equation (9) is related to the time step [tau], in order to ensure stability and accuracy, the time step [tau] should be selected carefully; it can not be chosen too small.

4. Numerical Examples

In this section we give several numerical examples in both 2D and 3D to demonstrate the effectiveness and stability of the new computation method. In the computation, we impose the noise to the final temperature condition. We set

[g.sub.[delta]](x) = g(x) x (1 + [delta] x rand n), (13)

where g(x) denotes the exact final temperature condition, [delta] is the tolerated noise level, and rand n is Gaussian random number with mean 0 and variance 1.

Two kinds of errors are used to compare the accuracy of the numerical solutions with the exact solutions. The maximum error (Maxerror) and the root-mean-square error (RMSE) are defined as follows,

[mathematical expression not reproducible], (14)

where M is the total number of testing nodes within the domain and [[??].sub.i]; denotes the approximate solution at the ith node; [u.sub.i] is the exact one.

Example 1. In this example we consider the two-dimensional nonhomogeneous heat conduction equation (1), with the Dirichlet boundary conditions,

u (x, t) = (y sin ([pi]x) + x cos ([pi]y)) cos t,

x = (x, y) [member of] [GAMMA], t [member of] (0, T) (15)

the final conditions,

u(x, T) = (y sin ([pi]x) + x cos ([pi]y)) cos T, x [member of] [OMEGA], (16)

and f(x,t) = [[pi].sup.2](y sin([pi]x) + x cos(ny)) cos t - (y sin([pi]x) + x cos([pi]y)) sin t.

The computational domain [bar.[OMEGA]] is a given rectangle in ([bar.x], [bar.y]) space such that 0 < a [less than or equal to] [bar.x] [less than or equal to] b and c [less than or equal to] [bar.y] [less than or equal to] d, using the map

x = [bar.x] cos ([theta][bar.y]),

y = [bar.x] sin ([theta][bar.y]), (17)

to take a point in the rectangular ([bar.x], [bar.y]) geometry and to map it to a point (x, y) in a big hollow cylinder.

In the computation we choose the mesh 41 x 81 on the rectangular and using the same map to put the mesh on the hollow cylinder, the mesh on the computational domain with a = c = 0.5, b = 1.0, d = 1.5, and [theta] = [pi]/2 is shown in Figure 1.

The Maxerror and RMSE for T = 1.0, 2.0, and 3.0 are shown in Table 1 with noise level [delta] = [10.sup.-3] and various time steps [tau]. It can be seen from Table 1 that the proposed method performs well for solving the backward problem.

The figure of numerical initial solutions with T = 1.0, [tau] = 0.25, and [delta] = [10.sup.-3] is presented in Figure 2.

For further investigations, we give a research of the results with a disturbed [delta] of the final time conditions; the figure of RMSE with different noisy level [delta] with [tau] = 0.25 is present in Figures 3 and 4 for T = 1.0 and T = 2.0, respectively.

We also discuss the results under different grids; the Maxerror and RMSE for T = 1.0, 2.0 and 3.0 are shown in Table 2 with various grids with noise level [delta] = [10.sup.-3] and time steps [tau] = 0.25.

To illustrate the generality, we choose different variables for the computational domain. The figures of the exact and numerical initial solutions are given in Figures 5-7 with different variables of the computational domain. In this computation we choose the mesh 81 x 161, [tau] = 0.2 and [delta] = [10.sup.-3], T = 1.

Example 2. We also consider the 2D nonhomogeneous backward heat conduction problem in this example; the known conditions are presented as follows,

h (x, i) = sin ([pi]xy) exp (-t),

g (x) = sin ([pi]xy) exp (-T),

and f(x, t) = ([[pi].sup.2][y.sup.2] + [[pi].sup.2][x.sup.2] - 1) sin([pi]ry) exp(-t). (18)

In this example, the calculation area is a square of the center at origin; side length is 1. In the computation using the mesh 51 x 51, the computational grid is shown in Figure 8. The Maxerror and RMSE for T = 1.0, 2.0, 3.0 are shown in Table [delta] with [delta] = [10.sup.-3] and various [tau].

The errors of RMSE with different [delta] are also given in Figure 9. In this figure [tau] = 0.2, and the mesh is the same as Figure 8.

The same with Example 1, we also consider the results with different meshes; the Maxerror and RMSE for T = 1.0, 2.0 and 3.0 are shown in Table 4 with various meshes with noise level [delta] = [10.sup.-3] and time steps r = 0.2.

The figure of numerical solutions m(x, y, 0) with T = 1, [tau] = 0.2, and [delta] = [10.sup.-3] is presented in Figure 10; the computational mesh is 51 x 51.

We also consider a slightly more complicated computational area, one unit square to remove a small rectangle and two small circles. In this case, the computation is a little more complicated; we define the unite square as [[OMEGA].sub.1] and the small rectangle as [[OMEGA].sub.2],

[[OMEGA].sub.1] = {(x, y) | [less than or equal to] x, y [less than or equal to] 1}

[[OMEGA].sub.2] = {(x, y) | 0.3 < x < 0.4, 0.2 < y < 0.5}, (19)

the two small circles as

[[OMEGA].sub.3] = {(x, y) | [(x - 0.1).sup.2] + [(y - 0.9).sup.2] < 0.052},

[[OMEGA].sub.4] = {(x, y) | [(x - 0.7).sup.2] + [(y - 0.6).sup.2] < 0.12}, (20)

and the computational domain is defined as

[bar.[OMEGA]] = [bar.[[OMEGA]].sub.1] - [[OMEGA].sub.2] - [[OMEGA].sub.3] - [[OMEGA].sub.4]. (21)

The figure of the domain with the computational mesh can be seen in Figure 11 and the figure of solutions u(x, y, 0) ia also given in Figure 12 with T = 3, [tau] = 0.2.

Example 3. In this example we consider the inverse problem in three dimensions with the Dirichlet boundary conditions,

h (x, t) = sin X sin y sin z cos t,

x = (x, y, z) [member of] [GAMMA], t [member of] (0, T), (22)

the final conditions,

g(x) = sin x sin y sin z cos T, x = (x, y, z) [member of] [OMEGA], (23)

and f(x, t) = sin x sin y sin z(3 cos t - sin t); the analytical solution is

u (x, t) = sin x sin y sin z cos t. (24)

In this example we choose the 3D computational domain defined as a cube without a cylinder insider. The cube is defined as

[[OMEGA].sub.1] = {(x, y, z) | -0.5 [less than or equal to] x, y [less than or equal to] 0.5, -1 [less than or equal to] z [less than or equal to] 1}, (25)

the cylinder is defined as

[[OMEGA].sub.1] = j(x, y, z) | [(4x).sup.2] + [(4y).sup.2] <1, -1 [less than or equal to] z [less than or equal to] 1}, (26)

and the computational domain is [bar.[OMEGA]] = [[OMEGA].sub.1] - [[OMEGA].sub.2]. The domain with the computational grid is shown in Figure 13.

The Maxerror and RMSE for T = 1.0, 3.0, 5.0 are shown in Table 5 with noise level [delta] = [10.sup.-3] and various [tau].

The figure of numerical initial solutions with T = 1, [tau] = 0.2 with noisy level [delta] = [10.sup.-3] can be seen in Figure 14.

Example 4. In this example we also consider the inverse problem with the same conditions in Example 3, but the computational domain is a sphere with the center at origin and the diameter is 1. The domain with the computational grid is shown in Figure 15.

The Maxerror and RMSE for T = 1.0,2.0,3.0 are shown in Table 6 with noise level [delta] = [10.sup.-3] and various [tau].

The figure of numerical initial solutions with T = 3, [tau] = 0.2 with noisy level [delta] = [10.sup.-3] can be seen in Figure 16.

5. Conclusion

In this paper, we proposed a new numerical scheme to solve nonhomogeneous backward heat conduction problem with two- and three-dimensional computational domain. In our method, using likewise Crank Nicolson scheme and introducing a new intermediate variable, transforming the backward problem to a nonhomogeneous Helmholtz problem, and solving the Helmholtz problem using finite element method, we get the initial temperature. This new method is easy to calculate and from the numerical results presented in the previous it can be seen that this new method is effective for solving such problems.

It is worth to note that, due to the limitations of the standard finite element method in solving Helmholtz equations with high wave numbers, to solve the Helmholtz problem (9) by FEM, the time step r should be carefully chosen, in order to ensure the wave number in (9) is not high, and for high dimension problem, dense grids need to be generated to ensure accuracy, resulting in increased computational time. Many works have been done in solving Helmholtz problem with high wave numbers [25-28]; these methods can be used instead of the standard FEM to solve the Helmholtz problem (9), so that these problems will be somewhat improved.

https://doi.org/10.1155/2018/1868921

Data Availability

All the data used to support the findings of this study are included within the article.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

The authors would like to to thank professors P. N. Vabishchevish and V. I. Vasilev for their important guidance and advice in the paper writing. This work is supported by RFBR (project 17-01-00732).

References

[1] J. Fan and L. Wang, "Review of heat conduction in nanofluids," Journal of Heat Transfer, vol. 133, no. 4, pp. 801-815, 2011.

[2] A. Paliathanasis, R. M. Morris, and P. G. L. Leach, "Lie symmetries of (1+2) nonautonomous evolution equations in financial mathematics," Mathematics, vol. 4, no. 2, pp. 1-14, 2016.

[3] M. A. Bartoshevich, "A heat-conduction problem," Journal of Engineering Physics, vol. 28, no. 2, pp. 240-244, 1975.

[4] O. M. Alifanov, Inverse Heat Transfer Problems, International Series in Heat and Mass Transfer, Springer-Verlag, 1994.

[5] C.-H. Huang and S.-P. Wang, "A three-dimensional inverse heat conduction problem in estimating surface heat flux by conjugate gradient method," International Journal of Heat and Mass Transfer, vol. 42, no. 18, pp. 3387-3403, 1999.

[6] N. D. Aparicio and M. K. Pidcock, "The boundary inverse problem for the Laplace equation in two dimensions," Inverse Problems, vol. 12, no. 5, pp. 565-577, 1996.

[7] V. V. Vasilev, M. V. Vasilyeva, and A. M. Kardashevskyc, "The numerical solution of the boundary inverse problem for a parabolic equation," American Institute of Physics, vol. 1773, no. 1, pp. 100010(1)-100010(7), 2016.

[8] L. Su, P. N. Vabishchevish, and V. I. Vasilev, "The inverse problem of the simultaneous determination of the right-hand side and the lowest coefficients in parabolic equations," in Numerical Analysis and Its Applications, Lecture Notes in Computer Science, pp. 633-639, 2017.

[9] P. N. Vabishchevich and V. I. Vasilev, "Computational algorithms for solving the coefficient inverse problem for parabolic equations," Inverse Problems in Science and Engineering, vol. 24, no. 1, pp. 42-59, 2016.

[10] A. A. Samarskii and P. N. Vabishchevich, Numerical Methods for Solving Inverse Problems of Mathematical Physics, Inverse and Ill-Posed Problems Series, Walter de Gruyter, 2007.

[11] M. Li, T. Jiang, and Y. C. Hon, "A meshless method based on RBFs method for nonhomogeneous backward heat conduction problem," Engineering Analysis with Boundary Elements, vol. 34, no. 9, pp. 785-792, 2010.

[12] J. Hadamard, Lectures on Cauchy Problems in Linear Partial Differential Equations, New Haven Yale University Press, 1923.

[13] L. Payne, Improperly Posed Problems in Partial Differential Equations, Society for Industrial and Applied Mathematics, 1975.

[14] W. L. Miranker, "A well posed problem for the backward heat equation," Proceedings of the American Mathematical Society, vol. 12, no. 2, pp. 243-247, 1961.

[15] C.-L. Fu, X.-T. Xiong, and Z. Qian, "Fourier regularization for a backward heat equation," Journal of Mathematical Analysis and Applications, vol. 331, no. 1, pp. 472-480, 2007.

[16] H. Han, D. B. Ingham, and Y. Yuan, "The boundary element method for the solution of the backward heat conduction equation," Journal of Computational Physics, vol. 116, no. 2, pp. 292-299, 1995.

[17] K. Iijima, "Numerical solution of backward heat conduction problems by a high order lattice-free finite difference method," Journal of the Chinese Institute of Engineers, vol. 27, no. 4, pp. 611-620, 2004.

[18] A. Shidfar and A. Zakeri, "A numerical technique for backward inverse heat conduction problems in one-dimensional space," Applied Mathematics and Computation, vol. 171, no. 2, pp. 1016-1024, 2005.

[19] D. Lesnic, L. Elliott, and D. B. Ingham, "An iterative boundary element method for solving the backward heat conduction problem using an elliptic approximation," Inverse Problems in Science and Engineering, vol. 6, no. 4, pp. 255-279, 1998.

[20] N. S. Mera, L. Elliott, D. B. Ingham, and D. Lesnic, "Iterative boundary element method for solving the one-dimensional backward heat conduction problem," International Journal of Heat and Mass Transfer, vol. 44, no. 10, pp. 1937-1946, 2001.

[21] Y. C. Hon and M. Li, "A discrepancy principle for the source points location in using the MFS for solving the BHCP," International Journal of Computational Methods, vol. 6, no. 2, pp. 181-197, 2009.

[22] M. Denche and A. Abdessemed, "On regularization and error estimates for non-homogeneous backward Cauchy problem," Arab Journal of Mathematical Sciences, vol. 18, no. 2, pp. 149-164, 2012.

[23] N. H. Tuan and D. D. Trong, "A new regularized method for two dimensional nonhomogeneous backward heat problem," Applied Mathematics and Computation, vol. 215, no. 3, pp. 873-880, 2009.

[24] S. Y. Reutskiy, C. S. Chen, and H. Y. Tian, "A boundary mesh-less method using Chebyshev interpolation and trigonometric basis function for solving heat conduction problems," International Journal for Numerical Methods in Engineering, vol. 74, no. 10, pp. 1621-1644, 2008.

[25] F. Ihlenburg and I. Babuska, "Finite element solution of the Helmholtz equation with high wave number Part I: The h-version of the FEM," Computers & Mathematics with Applications, vol. 30, no. 9, pp. 9-37, 1995.

[26] F. Ihlenburg and I. Babuka, "Finite element solution of the Helmholtz equation with high wave number part II: the h-p Version of the FEM," SIAM Journal on Numerical Analysis, vol. 34, no. 1, pp. 315-358, 1997.

[27] I. Harari and F. Magoules, "Numerical investigations of stabilized finite element computations for acoustics," Wave Motion, vol. 39, no. 4, pp. 339-349, 2004.

[28] Y. S. Wong and G. Li, "Exact finite difference schemes for solving Helmholtz equation at any wavenumber," International Journal of Numerical Analysis & Modeling- Series B, vol. 2, no. 1, pp. 91-108, 2011.

[29] I. Babuka and S. A. Sauter, "Is the pollution effect of the fem avoidable for the Helmholtz equation considering high wave," SIAM Journal on Numerical Analysis, vol. 34, no. 1, pp. 2392-2423, 1997.

[30] A. Bayliss, C. I. Goldstein, and E. Turkel, "On accuracy conditions for the numerical computation of waves," Journal of Computational Physics, vol. 59, no. 3, pp. 396-404, 1985.

LingDe Su (iD) (1) and TongSong Jiang (iD) (2)

(1) Institute of Mathematics and Information Science, North-Eastern Federal University, Russia

(2) Department of Mathematics, Heze University, Shandong, China

Correspondence should be addressed to TongSong Jiang; jiangtongsong@sina.com

Received 5 June 2018; Revised 14 October 2018; Accepted 15 November 2018; Published 2 December 2018

Guest Editor: Ovidiu Bagdasar

Caption: FIGURE 1: The domain with computational grid.

Caption: FIGURE 2: The numerical solutions u(x,0).

Caption: FIGURE 3: The errors of RMSE with different [delta], T = 1.0.

Caption: FIGURE 4: The errors of RMSE with different [delta], T = 2.0.

Caption: FIGURE 5: The figure with a = 1.0, b = 1.5, c = 0, d = 1.0, and [theta] = [pi]/4.

Caption: FIGURE 6: The figure with a = 1.0, b = 1.5, c = 1.0, d = 2.0, and [theta] = [pi]/2.

Caption: FIGURE 7: The figure with a = 0.5, b = 1.0, c = -0.5, d = 0.5, and [theta] = [pi]/4.

Caption: FIGURE 8: The square domain with computational grid.

Caption: FIGURE 9: The errors of RMSE with different [delta], T = 1.0 and T = 3.0.

Caption: FIGURE 10: The numerical solutions of u(x, y, 0).

Caption: FIGURE 11: The domain with computational grid.

Caption: FIGURE 12: The figures of exact and numerical solutions of m(x, y, 0).

Caption: FIGURE 13: The 3D domain with computational grid.

Caption: FIGURE 14: The numerical solutions w(x, y, 0) with T = 1, [tau] = 0.2.

Caption: FIGURE 15: The sphere domain with computational grid.

Caption: FIGURE 16: The numerical solutions m(x, y, 0) with T = 3, [tau] = 0.2.
```TABLE 1: Maxerror and RMSE for T = 1.0, 2.0 and 3.0 with different
[tau].

[tau] = 0.2            [tau] = 0.25

T = 1.0
Maxerror     4.951 x [10.sup.-4]    3.714 x [10.sup.-4]
RMSE         2.230 x [10.sup.-4]    2.358 x [10.sup.-4]

T = 2.0
Maxerror     2.393 x [10.sup.-3]    2.041 x [10.sup.-3]
RMSE         1.172 x [10.sup.-3]    1.016 x [10.sup.-3]

T = 3.0
Maxerror     2.336 x [10.sup.-3]    1.620 x [10.sup.-3]
RMSE         1.046 x [10.sup.-3]    1.766 x [10.sup.-3]

[tau] = 0.5            [tau] = 1.0

T = 1.0
Maxerror     8.721 x [10.sup.-4]    1.635 x [10.sup.-2]
RMSE         4.734 x [10.sup.-4]    8.604 x [10.sup.-3]

T = 2.0
Maxerror     4.100 x [10.sup.-3]    1.569 x [10.sup.-2]
RMSE         2.133 x [10.sup.-3]    8.278 X [10.sup.-3]

T = 3.0
Maxerror     6.557 x [10.sup.-3]    3.554 x [10.sup.-3]
RMSE         3.320 x [10.sup.-3]    8.416 x [10.sup.-4]

TABLE 2: Maxerror and RMSE for T = 1.0,2.0,3.0
with different grids.

21 x 41                41 x 81

T = 1.0
Maxerror     3.256 x [10.sup.-4]    2.205 x [10.sup.-4]
RMSE         6.786 x [10.sup.-5]    1.358 x [10.sup.-4]

T = 2.0
Maxerror     1.362 x [10.sup.-3]    1.665 x [10.sup.-3]
RMSE         2.888 x [10.sup.-4]    8.172 x [10.sup.-4]

T = 3.0
Maxerror     2.061 x [10.sup.-3]    5.661 x [10.sup.-3]
RMSE         5.832 x [10.sup.-4]    2.751 x [10.sup.-3]

81 x 161              161 x 321

T = 1.0
Maxerror     1.466 x [10.sup.-4]    3.418 x [10.sup.-4]
RMSE         8.334 x [10.sup.-5]    2.276 x [10.sup.-4]

T = 2.0
Maxerror     2.013 x [10.sup.-3]    1.982 x [10.sup.-3]
RMSE         1.020 x [10.sup.-3]    1.009 x [10.sup.-3]

T = 3.0
Maxerror     3.301 x [10.sup.-3]    7.441 x [10.sup.-3]
RMSE         1.565 x [10.sup.-3]    3.705 x [10.sup.-3]

TABLE 3: Maxerror and RMSE for T = 1, 2, 3 with different [tau].

t = 1.0                t = 0.5

T = 1.0
Maxerror     1.301 x [10.sup.-3]    1.374 x [10.sup.-4]
RMSE         5.907 x [10.sup.-4]    5.727 x [10.sup.-5]

T = 2.0
Maxerror     7.725 x [10.sup.-4]    1.513 x [10.sup.-4]
RMSE         3.463 x [10.sup.-4]    6.095 x [10.sup.-5]

T = 3.0
Maxerror     9.877 x [10.sup.-4]    2.187 x [10.sup.-4]
RMSE         4.474 x [10.sup.-4]    9.301 x [10.sup.-5]

t = 0.25               t = 0.2

T = 1.0
Maxerror     2.504 x [10.sup.-4]    4.221 x [10.sup.-4]
RMSE         1.409 x [10.sup.-4]    2.280 x [10.sup.-4]

T = 2.0
Maxerror     6.912 x [10.sup.-5]    2.041 x [10.sup.-4]
RMSE         3.407 x [10.sup.-5]    1.015 x [10.sup.-4]

T = 3.0
Maxerror     6.443 x [10.sup.-5]    4.832 x [10.sup.-4]
RMSE         1.911 x [10.sup.-5]    2.457 x [10.sup.-4]

TABLE 4: Maxerror and RMSE for T = 1.0, 2.0, 3.0 with different grids.

21 x 21                51 x 51

T = 1.0
Maxerror     4.767 x [10.sup.-4]    4.114 x [10.sup.-4]
RMSE         1.409 x [10.sup.-4]    2.221 x [10.sup.-4]

T = 2.0
Maxerror     2.644 x [10.sup.-4]    2.124 x [10.sup.-4]
RMSE         4.982 x KT5            1.058 x [10.sup.-4]

T = 3.0
Maxerror     5.485 x [10.sup.-4]    4.497 x [10.sup.-4]
RMSE         2.807 x [10.sup.-4]    2.288 x [10.sup.-4]

101 x 101              201 x 201

T = 1.0
Maxerror     9.978 x [10.sup.-5]    1744 x [10.sup.-4]
RMSE         4.938 x [10.sup.-5]    9.225 x [10.sup.-5]

T = 2.0
Maxerror     3.899 x [10.sup.-4]    3.138 x [10.sup.-4]
RMSE         1.986 x [10.sup.-4]    1.621 x [10.sup.-4]

T = 3.0
Maxerror     5.738 x [10.sup.-4]    6.567 x [10.sup.-4]
RMSE         2.925 x [10.sup.-4]    3.353 x [10.sup.-4]

TABLE 5: Maxerror and RMSE for T = 1,2,5 with different t.

[tau] = 1.0            [tau] = 0.5

T = 1.0
Maxerror     4.726 x [10.sup.-4]    2.014 x [10.sup.-3]
RMSE         5.653 x [10.sup.-5]    1.499 x [10.sup.-4]

T = 3.0
Maxerror     1.807 x [10.sup.-3]    1.161 x [10.sup.-4]
RMSE         1.983 x [10.sup.-4]    3.789 x [10.sup.-5]

T = 5.0
Maxerror     6.437 x [10.sup.-4]    1.125 x [10.sup.-3]
RMSE         6.681 x [10.sup.-5]    1.225 x [10.sup.-4]

[tau] = 0.25           [tau] = 0.2

T = 1.0
Maxerror     4.399 x [10.sup.-4]    1.409 x [10.sup.-3]
RMSE         4.827 x [10.sup.-5]    1.489 x [10.sup.-4]

T = 3.0
Maxerror     1.717 x [10.sup.-3]    2.961 x [10.sup.-4]
RMSE         1.974 x [10.sup.-4]    6.815 x [10.sup.-5]

T = 5.0
Maxerror     6.596 x [10.sup.-4]    1.206 x [10.sup.-3]
RMSE         7.013 x [10.sup.-5]    1.629 x [10.sup.-4]

TABLE 6: Maxerror and RMSE for T = 1,2,3 with different r.

[tau] = 1/2            [tau] = 1/3

T = 1.0
Maxerror     3.765 x [10.sup.-5]    1.119 x [10.sup.-4]
RMSE         7.204 x [10.sup.-6]    1.780 x [10.sup.-5]

T = 2.0
Maxerror     1.021 x [10.sup.-4]    1.028 x [10.sup.-4]
RMSE         1.742 x [10.sup.-5]    1.645 x [10.sup.-5]

T = 3.0
Maxerror     1.431 x [10.sup.-4]    1.364 x [10.sup.-5]
RMSE         2.497 x [10.sup.-5]    3.921 x [10.sup.-6]

[tau] = 1/4            [tau] = 1/5

T = 1.0
Maxerror     3.521 x [10.sup.-5]    1.119 x [10.sup.-4]
RMSE         6.112 x [10.sup.-6]    1.766 x [10.sup.-5]

T = 2.0
Maxerror     1.025 x [10.sup.-4]    1.031 x [10.sup.-4]
RMSE         1.658 x [10.sup.-5]    1.638 x [10.sup.-5]

T = 3.0
Maxerror     1.463 x [10.sup.-4]    3.104 x [10.sup.-4]
RMSE         2.324 x [10.sup.-5]    1.172 x [10.sup.-5]
```