## Abstract

This study presents an approach for reconstruction of harmonic functions in three dimensions from the finite number of field and surface measurements. The approach, based on the Trefftz method, performs reconstruction as the best fit to the data and provides smoothness of the reconstructed function. Two particular algorithms are proposed; the first one uses specific radial basis functions and the second one is of finite element type. Either of them can be applied to analyse different data types but the latter can handle larger problems. The data types considered in this study also cover direct and inverse boundary value problems. Therefore, the proposed approach is universal and capable of dealing with both well-posed and ill-posed formulations. Examples from steady heat conduction and elastostatics are examined in order to investigate the efficiency of the approach.

## 1. Introduction

Solutions of diverse potential problems can be presented in terms of harmonic functions, i.e. by those satisfying the Laplace equation. Numerical methods for solving direct potential problems are very well developed for well-posed formulations; however, in many applications, these are not possible. Such a wide class of problems is usually referred to as inverse or indirect problems; they require special numerical approaches (see general discussion in Tikhonov & Arsenin (1977)). Inverse problems often involve field measurements in order to compensate for the lack of boundary data, which mathematically leads to ill-posed formulations. This study is focused on the reconstruction of harmonic functions in three dimensions from a finite number of field measurements. Such problems do not possess unique solutions, which makes them, in a certain sense, similar to interpolation problems. For dense regular data, three-dimensional interpolation techniques have been developed (Lekien & Marsden 2005); however, they do not obey any physical equations, which restricts application of these techniques. We do not attempt to fit the measurements exactly but rather find an approximate solution (satisfying governing equations) that provides the best fit for the data. This is in line with modern approaches used for solving inverse potential problems; a brief review of these is presented in this section.

Steady heat conduction is a classical example of potential problems. The direct problem, in this case, requires the determination of temperature at interior points of a domain by prescribed boundary conditions on temperature or heat flux. In many practical situations, boundary conditions cannot be formulated on the entire boundary and have to be complemented by experimental measurements, which lead to an inverse heat conduction problem (IHCP). It should be mentioned that the measurements may be seriously affected by sensors, which leads to loss of accuracy. Alternatively, the surface may be inaccessible for sensors and interior measurements are easier to obtain. A review on IHCP prior to 1985 is given by Beck *et al*. (1985).

Earlier studies on the Cauchy problem for the Laplace equation were carried out by Hadamard (1923), who had introduced the well-known definition for ill-posed mathematical problems. Apart from non-uniqueness, instability of the solution is associated with the ill-posed problems. In numerical implementations, it leads to an ill-conditioned matrix and standard numerical methods cannot achieve satisfactory accuracy and stabilize the solution as a result of a large condition number of the matrix. Several regularization methods have been developed for solving these problems. The most common regularization technique (Tikhonov & Arsenin 1977) introduces a regularization parameter whose optimum can be determined by other techniques, such as the L-curve technique (Lawson & Hanson 1995). Other common regularization methods are based on singular value decomposition (SVD; e.g. Hansen 1991), iterative conjugate gradient method (CGM) or least squares QR factorization (LSQR) algorithms (Paige & Saunders 1982) and the minimal energy method (Han *et al*. 1995).

Many attempts to solve the ill-posed IHCP are based on the finite element method (FEM; Zienkiewicz & Taylor 1988) and the boundary element method (BEM; Brebbia 1978) techniques. Brief review of applicability of these methods to inverse problems is found in Mackerle & Tanaka (1990). It has been demonstrated that both these methods can be applied to elliptic (steady heat conduction) and parabolic (unsteady heat conduction) equations related to IHCP. For instance, Han *et al*. (1995) and Lesnic *et al*. (1994) considered unsteady heat conduction by using BEM formulation and the minimal energy method. Ingham (1996) investigated two-dimensional problems for both steady and unsteady heat conduction conditions involving measured data by using a BEM formulation and the minimal energy method for solving the system. Cimetiere *et al*. (2002) reported that both FEM and BEM approaches provide satisfactory results in the problem of the determination of harmonic functions by overspecified data on a part of the domain boundary. In some situations, only internal data can be provided or some additional internal data can be obtained and used in order to improve accuracy of solutions. This situation has been studied by Pasquetti & Petit (1995) for the diffusion equation by applying BEM complemented by different minimization techniques.

Besides FEM and BEM, there are other techniques such as the method of fundamental solutions, meshless approaches, the method of control volume, integration free, etc. (see the survey by Babuska *et al*. (2003)). The idea for the method of fundamental solutions (method of functional equations) was suggested by Kupradze & Aleksidze (1964). This method can be referred to as the meshless and integration-free method. Its application to IHCP has been reported by Hon & Wei (2004) for one-dimensional unsteady heat conduction with the use of the Tikhonov regularization and L-curve criterion. Recent developments of this method to solve the IHCP for different problems are found in Yan *et al*. (2008) and Shidfar *et al*. (2009). The method of control volume was used in a three-dimensional IHCP for the estimation of surface heat fluxes from temperature interior points by Huang & Wang (1999). Huang *et al*. (2007) used the steepest descent method to simulate the heat flux of titanium drilling and reported comparisons with real experimental data. Different variants of the Trefftz method are often used in applications (see surveys by Kita & Kamiya 1995; Jirousek & Zielinski 1997) and hybrid types of FEM and Trefftz approaches are discussed (e.g. Qin 2000) for two-dimensional boundary value problems. Kompis (1992) suggested using polynomial approximations within the elements that satisfy governing equations to solve plane boundary value problems. The present study generalizes this approach into the three-dimensional case and for inverse problems. A somewhat similar hybrid approach has been applied by Irša & Galybin (2009*a*) for heat flux reconstruction from temperature measurements inside a two-dimensional region. The method is actually of the Trefftz type since solutions in every element are sought as holomorphic functions whose real parts (presenting temperature) satisfy the Laplace equation. The method employs complex variables and therefore the approach is restricted to two-dimensional problems only.

Inverse potential problems are also found in other areas of applied mathematics; for instance, in fluid flow and elasticity. They are often solved by the methods mentioned earlier. For instance, Zeb *et al*. (2000) applied BEM to an inverse Stokes problem and used the Tikhonov regularization and L-curve criterion for the determination of fluid velocities from boundary data and interior pressure measurements. In elasticity, displacement measurements may complement boundary conditions (Zabaras *et al*. 1989), which leads to a minimization problem of fitting the measurements by the least squares method. Other examples include estimation of residual and contact stresses (Zhang *et al.* 1997) or fracture parameters (e.g. Galybin 2008).

Despite considerable progress in the investigation of inverse problems of different types, there is still a gap in the development of numerical methods capable of dealing with finite data, i.e. with such situations where the problem is strongly underdetermined and, thus, has no unique solution. Such formulations have been considered for the problem of stress-state identification in the Earth’s crust by discrete data on stress orientations. Galybin & Mukhamediev (2004) suggested an approach for two-dimensional cases that has been validated by simulating photoelasticity data. The approach is based on the approximation of complex potentials by polynomials, which restricts higher order approximations for noisy data. One more stable modification of this approach has recently been reported by Irša & Galybin (2009*b*) but it is also limited to a two-dimensional case.

The current study is aimed at inverse three-dimensional problems described by harmonic functions; it is assumed that the data are discrete.

## 2. Approximation problems

The general approximation problem considered further on is formulated as follows. Given data at *N* points of an isotropic three-dimensional domain, *Ω*, reconstruct *K* harmonic functions, *φ*_{k}(**x**),**x**=(*x*_{1},*x*_{2},*x*_{3}), whose known linear combinations provide the best fit for the data.

Let us form a column of unknown harmonic functions **h**(**x**)=(*φ*_{1}(**x**),…, *φ*_{K}(**x**))^{T} and introduce a *J*×*K* matrix **A**(**x**)={*A*_{j,k}(**x**)}, then the product **A**(**x**)**h**(**x**) presents an array **L**(**x**) of *J* linear combinations of *φ*_{k}(**x**), i.e. **L**(**x**)=(*L*_{1}(**x**),…,*L*_{J}(**x**))^{T}, with the components
2.1
where *A*_{j,k}(**x**) are given spatial smooth (double differentiable) functions.

It is assumed that at every datum point, **t**_{n}, at least one combination *L*_{j}(**x**) is known. For convenience, one can associate one combination with a single datum point; thus, the same datum point has to be counted *m* times in the dataset {**t**_{n},*n*=1,…,*N*} if *m* combinations are known at this point. Hence, *n*=*j* and *N*=*J*. Therefore, the following constraints are imposed on *L*_{j}(**x**):
2.2
where *f*_{n} are known.

Let us seek the unknown harmonic functions as certain linear combinations of the form
2.3
Here, the coefficients *c*_{k,m} are unknown constants, while *χ*_{m} are known independent harmonic functions that are further referred to as the basis functions (hereafter, the number of basis functions is always *M*). By stacking the coefficients *c*_{k,m} in a single column of the length *KM* as follows: **C**=(*c*_{1,1},…,*c*_{1,M},*c*_{2,1},…,*c*_{2,M},…,*c*_{K,1},…,*c*_{K,M})^{T}, one can express the column of unknown harmonic functions in the form
2.4
where the matrix **T**(**x**) is defined as follows:
2.5
Substitution of equation (2.4) into equation (2.1) determines linear combinations in the matrix form **L**(**x**)=**A**(**x**)**T**(**x**)**C**, which, being substituted into equation (2.2), leads to the following minimization problem for the determination of the unknown coefficients *c*_{k,m}
2.6
As soon as the minimization problem is solved and vector **C** is determined, one finds the set of sought harmonic function from equation (2.4).

Based on the general formulation, two particular engineering problems are investigated further on in detail. The first one is referred to as the heat reconstruction problem (HRP), in which temperature measurements are used as the data. The second one attempts to reconstruct elastic stress states from displacement measurements; it is referred to as the elastic reconstruction problem (ERP).

### (a) Heat reconstruction problem (HRP )

It is assumed that the temperature *φ* (satisfying Δ*φ*(**x**)=0) is measured at certain data points, which form the dataset {*φ*_{n},*n*=1,…,*N*}. No data on heat flux are imposed for simplicity. Therefore, the problem of reconstruction of a single harmonic function *φ*(**x**) can be radically reduced to the following system of linear algebraic equations:
2.7
Once the coefficients are found, the solution of the problem can be written by the dot product
2.8
where **C**=(*c*_{1},…,*c*_{M})^{T} is the array of coefficients.

### (b) Elastic reconstruction problem (ERP )

The Lame equations in three-dimensional elasticity have the form (no body forces; e.g. Timoshenko & Goodier 1970)
2.9
where *u*_{k} are the components of the displacement vector, *ν* is Poisson’s ratio, operator ∂_{k} stands for the partial differentiation with respect to *x*_{k} and *e* is the dilation specified as follows:
2.10
The dilation is a harmonic function in the absence of body forces.

The general solution for displacements in an isotropic elastic body (often referred to as the Papkovitch–Neuber solution) can be expressed via harmonic functions (see Timoshenko & Goodier (1970) for details). There exist several equivalent forms for the general solution, one of which (found in Leonov (1981)) is the most suitable for illustration of the proposed technique for contact problems. It has the form
2.11
where *φ*_{k} are arbitrary harmonic functions and the harmonic function *Π* is not an independent one—it is connected with the dilation as *e*(**x**)=∂_{3}*Π*(**x**) and, therefore, the following relationship between the four harmonic functions holds:
2.12

The components of the strain and stress tensors can be expressed via three independent harmonic functions *φ*_{k} by the strain-displacement relationships and Hooke’s Law.

Taking into account that 2∂_{k}*e*=Δ(*x*_{k}*e*), one can modify equation (2.11) as follows:
2.13
which is more convenient for the further analysis.

Assume that displacements and strain/stress are monitored at some points within the body and/or its boundary. Then the ERP can be formulated in a way similar to the HRP, but, in construct to the latter, three piecewise harmonic functions should be determined instead of one. This would lead to a system that involves three times more unknowns than that in the case of a single harmonic function reconstruction.

In some practical cases, one can consider two consecutive problems for determination of *e* and *φ*_{k} when the dilation is measured at some points (such techniques exist; for instance, the methods based on overcoring (Sjoberg *et al.* 2003) are able to recover the full strain tensor in a point of rockmass). If this information is known, then the procedure for ERP is as follows.

— The harmonic function

*e*is reconstructed in the considered volume by using the technique described for the HRP.— As soon as the dilation is determined, the data for approximation of harmonic functions

*φ*_{k}are formed from equation (2.13) at the points where displacements were measured.— The harmonic functions

*φ*_{k}are found separately by solving the same system for three different right-hand sides.

This procedure will be employed in an example for the contact problem. Therefore, we further focus on the reconstruction of a single harmonic function from system (2.7).

## 3. Uniqueness and connection with interpolation and boundary value problems

Reconstruction of harmonic functions from limited data is not unique. In two dimensions, this is evident from the fact that harmonic function can be considered as the real part of a holomorphic function, *F*(*z*), of complex variable *z*=*x*_{1}+*ix*_{2}. Indeed, if *φ*(*z*)=Re(*F*(*z*)) is a solution obtained from system (2.7), then it is evident that Re(*φ*(*z*)+*P*_{N}(*z*)*ψ*(*z*)) is also a solution that includes an arbitrary holomorphic function *ψ*(*z*) and the *N*th degree polynomial *P*_{N}(*z*) that have the roots at data points and, therefore, the term Re(*P*_{N}(*z*)*ψ*(*z*)) does not contribute to the left-hand side of equation (2.7). Hence, the reconstruction can be performed with high uncertainty. However, in practice, reconstructed plane stress tensors from the data on principal directions obtained from photoelastic experiments show satisfactory correlation with experimental observations and even allow for the determination of isotropic points of different types in the stress trajectory fields (see Galybin & Mukhamediev (2004) for details). It will be seen later on that in the three-dimensional case the reconstruction gives reasonable results for a modest number of data as well. However, it is also non-unique, which can be demonstrated on the basis of the representation for the two-dimensional case as follows.

Let us consider a real valued function *g*(**x**)=*x*_{3}(*P*_{N}(*x*_{1}+*ix*_{2})*ψ*(*x*_{1}+*ix*_{2})+*P*_{N}(*x*_{1}−*ix*_{2})*ψ*(*x*_{1}−*ix*_{2})), where, similar to the two-dimensional case, *P*_{N}(*x*_{1}+*ix*_{2}) is a polynomial with the roots at the data points and *ψ*(*x*_{1}+*ix*_{2}) is an arbitrary analytical function. By inspection, it is obvious that Δ*g*(**x**)=0, i.e. this harmonic function has zeroes at the data points and, therefore, it can be added to any particular solution of equation (2.7). Other functions with the same property can be obtained from *g*(**x**) by interchanging the variables. Thus, it is evident that the reconstruction in three dimensions can be performed only with uncertainty and, therefore, the problem with finite data is somewhat similar to interpolation.

In some cases, when the dataset is not large, the interpolation can be performed, i.e. the reconstructed function passes exactly through all data points. This, however, does not guarantee uniqueness and the result will be dependent on the basis functions in the vector **h**(**x**). Similar to the one-dimensional case, large oscillations of the interpolant between the data points can be expected if interpolation involved many basis functions. Based on these analogies, one can view the problem considered here as an interpolation problem provided that the dataset is relatively small and the data are mostly within the domain.

The case when the data are sufficiently dense and located on the boundary presents the Dirichlet boundary value problem in a discrete form. For a large dataset, it would be unrealistic to interpolate all boundary data and, therefore, the solution at the data points is always approximate because a correspondent system of linear algebraic equations is overspecified (see next section). This does not necessarily mean that such a solution is less accurate than that obtained by conventional approaches that also have errors due to discretization and interpolation. Although this study is not aimed at solving boundary value problems, we inspect the applicability of the algorithms presented in the next section for direct formulations in heat conduction and elasticity.

Some inverse problems are also covered by the proposed approach. These assume that the boundary conditions are given on a part of the boundary and complemented by measurements at internal points. In some cases, such problems may have unique solutions, if boundary conditions are overspecified on a part of the boundary, i.e. two conditions are imposed on the boundary value of the sough function and its normal derivative (e.g. temperature and heat flux are given simultaneously on a part of the boundary). When these boundary conditions do not overlap, there is no uniqueness, which, however, does not put any restrictions on the reconstruction procedure.

In general, all types of harmonic problems associated with discrete three-dimensional data may be covered by the proposed reconstruction procedure regardless of non-uniqueness of this procedure.

## 4. Numerical approach

Equations (2.3)–(2.7) can be viewed as the application of the Trefftz approach (used for solving boundary value problems) to the approximation problem with finite data. This section presents two different numerical algorithms for solving the problem (2.7) based on either smooth or piecewise harmonic basis functions. In the latter case, certain continuity conditions across the adjacent elements are introduced in order to provide smoothness of the sought harmonic function in average.

### (a) Approximation by smooth harmonic functions

Let us first use the most convenient (at least for programming) method for presentation of the sough function via simple radial basis functions (*χ*_{m}(**x**)=*R*(**x**−**q**_{m}),Δ*R*(**x**−**q**_{m})=0) as follows:
4.1
Here parameters **q**_{m} can be selected as random points lying outside the domain *Ω*; in examples, they are placed on a sphere surrounding the domain.

By substituting the data points into equation (4.1), one obtains a system of linear algebraic equations
4.2
By introducing two columns formed from unknown coefficients **C**=(*c*_{1},…*c*_{M})^{T} and the data **f**=(*φ*_{1},…*φ*_{M})^{T}, one can present equation (4.2) in the symbolic form as
4.3
where the matrix **R** has *N* rows and *M* columns and its components *r*_{n,m} are calculated in accordance with equation (4.1), *r*_{n,m}=∥**x**_{n}−**q**_{m}∥^{−1}, *n*=1,…,*N*, *m*=1,…,*M*. This system is overspecified, *M*<*N*, and its approximate solution is found by the least squares method (Golub & van Loan 1989) as follows:
4.4

The quality of the approximate solution (4.4) is examined by monitoring two quantities. The residual (Res), calculated as the L_{2} norm, *Res*=∥**RC**−**f**∥_{2}, controls the accuracy while the condition number (Cn) controls stability. It is known (e.g. Kita & Kamiya 1995) that, for boundary value problems (both direct and inverse), this approach often leads to ill-conditioned matrices for large systems. This necessitates application of regularization techniques or special iterative algorithms. For overspecified matrices with moderate number of columns, one can readily compute the condition number as the ratio of the major and minor singular values in order to decide whether regularization is necessary or not. In the examples considered, further on the SVD regularization (Hansen 1991) has been performed if Cn>10^{8}. In this case, the approximate solution has the form:
4.5
where **U**(*N*×*M*) and **V** (*M*×*M*) are orthogonal matrices in the SVD decomposition **A**=**UDV**^{T} and **D**(*M*×*M*) is a diagonal matrix formed by the singular values, *d*_{j}, placed in descending order on the main diagonal , and the matrix **D**′ is the regularized inverse of **D** with the rank . Stability of the solution is controlled by the condition number of the regularized system determined as the ratio of *d*_{1}/*d*_{k} (Cn=*d*_{1}/*d*_{M} without regularization).

It should be noted that regularization allows stabilization of the system but cannot reduce the residual to zero. This, apparently, limits the application of this technique to large systems.

### (b) Approximation by piecewise harmonic functions: FEM-type approach

The second approach assumes discretization of the considered domain *Ω* into non-overlapping subdomains *Ω*_{m}, *m*=1…*M*, and approximation of the sought harmonic function by piecewise harmonic polynomials, P(**c**_{m},**x**) (the number of subdomains is equal to the number of the basis functions; therefore, denoted as *M*). Thus, the basis functions in equation (2.7) take the form
4.6
The polynomials P(**c**_{m},**x**) are assumed to be quadratic *P*(**c**_{m},**x**)=*P*_{q}(**c**_{m},**x**) or cubic *P*(**c**_{m},**x**)=*P*_{c}(**c**_{m},**x**) generalized three-dimensional polynomials satisfying the Laplace equation, i.e. of the following form:
4.7
The following arrays can be formed from the unknown coefficients *c*_{m,k}:
4.8
where *s*=8 for the quadratic approximation and *s*=15 for the cubic approximation.

By introducing the matrix functions
4.9
one can represent the polynomials in equation (4.6) in the matrix form
4.10
**Q**(**x**)=**Q**_{q}(**x**) for the quadratic approximation and **Q**(**x**)=(**Q**_{q}(**x**),**Q**_{c}(**x**)) for the cubic approximation.

Taking equation (4.10) into account, one presents the sough function in the form
4.11
Now, the first group of equations can be formed by fitting the data at *N* points **t**_{n}
4.12
The second group of equations obeys continuity of the sough function and its derivatives across the interfaces between adjacent blocks. For this purpose, let us introduce a collocation point **p**_{j} at the centre of every interface and impose the continuity conditions at this point. Owing to the properties of harmonic functions, it can be said that such positions of collocation points provide continuity over the interfaces in average.

The number of interfaces can be calculated readily for the case considered further on; namely, when all subdomains *Ω*_{m} are rectangular parallelepipeds of the same size (further referred to as blocks). Then, the computational domain can also be considered as a rectangular parallelepiped consisting of blocks *Ω*_{m}. It should contain the physical domain where the reconstruction is being performed. The number of collocation points, *N*_{cp}, depends on the partition of *Ω*; and it is calculated as *N*_{cp}=3*n*_{1}*n*_{2}*n*_{3}−*n*_{1}*n*_{2}−*n*_{1}*n*_{3}−*n*_{2}*n*_{3}, where the integers *n*_{j} are assigned to the number of blocks placed along the *x*_{j} axis. It is evident that if all *n*_{j}>1 then *N*_{cp} is greater than the number of all blocks, *M*=*n*_{1}*n*_{2}*n*_{3}. It should be noted that some blocks can have no associated collocation points belonging to the physical domain; in such cases, they may be excluded from consideration and then the total number of blocks (and collocation points) would be less than estimated earlier.

Since the existence of the second derivative is necessary in the harmonic problems, it would be natural to impose continuity of the sough function and all its derivatives up to the second order at all collocation points. This, however, cannot be applied for the quadratic polynomials because, in this case, all coefficients from *c*_{m,4} to *c*_{m,8} would be the same in all elements. Therefore, the continuity of the second-order derivatives is imposed for the cubic approximation only. On the one hand, it is expected that the cubic approximation could lead to better accuracy in the reconstruction, but, on the other hand, this approximation requires more computer recourses than the quadratic one. Indeed, every collocation point generates four equations for the quadratic approximation and nine equations for the cubic approximation (one of the second-order derivatives is not independent owing to the Laplace equation). The total number of equations expressing continuity is 4*N*_{cp} and 9*N*_{cp}, respectively. Bearing in mind that the number of unknowns is 9*M* (quadratic) or 16*M* (cubic), one can notice that the size of the system for cubic approximation is approximately four times larger than that for the quadratic one.

The second group of equations can be presented as follows:
4.13
Here, *α* and *β* designate the blocks with the common interface *j*, the column **O** consists of zeroes and the matrix **T** is
4.14
for the quadratic approximation and
4.15
for the cubic approximation.

Two groups of equations (4.12) and (4.13) form the whole system of linear algebraic equations that can be symbolically presented in the form
4.16
where the matrix **M** is organized such that being applied to the column of unknowns arranged as **C**=(**c**_{1},…**c**_{M})^{T} it results in the right-hand side **b**=(*φ*_{1},…*φ*_{N},0,…,0)^{T}. The dimensions of this matrix are (*N*+4*N*_{cp})×9*M* for the quadratic approximation and (*N*+9*N*_{cp})×16*M* for the cubic application.

The size of system (4.16) rapidly increases with the increase in the number of blocks; for instance, if *n*_{1}=*n*_{2}=*n*_{3}=10 then the number of unknowns is 16 000 and the number of equations is 24 300+*N* (for the cubic approximation). This makes it impossible to apply the least squares method in the form (4.3)–(4.5) and requires application of iterative methods. It has been reported (Paige & Saunders 1982) that the LSQR algorithm is efficient for large overdetermined and ill-conditioned systems. This algorithm is used further on when necessary.

The quality of approximate solution is controlled by the condition number (Cn) and the residual (Res) calculated in a similar way as for system (4.3).

## 5. Numerical experiments

This section presents some results of numerical experiments with data of different types in two potential three-dimensional problems, steady heat conduction and contact problem in elastostatics. The aim of these experiments is to investigate the accuracy in the reconstruction for the different data types and to test the sensitivity to experimental errors. In all experiments, the domain *Ω* is the three-dimensional cube of unit volume with the centre at point (0.5,0.5,0.5) and the edges parallel to the coordinate axes.

The following data types are considered.

— Dense data are given on the whole boundary; this case represents the Dirichlet boundary value problem (Dirichlet BVP).

— Dense data are given on a part of the boundary with several data inside the domain; this case corresponds to an inverse problem, and, further on, it is referred to as inverse BVP.

— Sparse data are distributed randomly inside the domain (including its boundary); this situation can be referred to as an interpolation problem.

We also study the performance of the proposed numerical schemes, i.e. the use of the radial basis functions for approximation and quadratic and cubic approximations in the finite element variant of the Trefftz method. In accordance with the methods used, the results are marked as follows: RBF, FE2 and FE3. In the RBF approach, the parameters **q**_{m} were placed at random points on the sphere of radius with the centre at point (0.5,0.5,0.5) that describes the domain *Ω*.

### (a) Heat reconstruction problem

The following ideal distribution of temperature is assumed: 5.1

This distribution corresponds to the case when the temperature is zero on all sides of the cube except on its top (*x*_{3}.=1), where it reaches a maximum of 100^{°} at point (0.5,0.5,1). Although solution (5.1) is simple, it is evident that a superposition of solutions of this type with different harmonics and interchange of variables covers the wide class of temperature distributions observed in practice.

Equation (5.1) is further used to compile synthetic data. Namely, the dataset is formed as the ideal solution with additional noise at *N* points **t**_{n},
5.2

Here, *η*_{n} are random numbers uniformly distributed on (−1,1) and *ε* is the noise level. The latter is set as *ε*=0.05 for the inverse problem (in case HRP-2 below) and as *ε*=0 for the other cases.

All three types of data mentioned earlier have been generated and used as input into the three numerical schemes. The reconstructed temperature fields are compared with the ideal distribution computed from equation (5.1). Contour plots and profiles along certain directions are used for illustration of the results of reconstructions. They mostly present absolute or relative errors.

Three profiles along the following directions are shown in the graphs.

— Profile 1 is along the diagonal of the top side between vertexes (0,0,1) and (1,1,1) that intersect the point (0.5,0.5,1), where the temperature reaches its maximum.

— Profile 2 is along the diagonal connecting vertexes (0,0,0) and (1,1,1).

— Profile 3 is along the line connecting the centres of the top and bottom sides, i.e. passing through the points (0.5,0.5,0) and (0.5,0.5,1).

Temperature profiles calculated from ideal solution (5.1) are shown in figure 1.

The following cases are studied further on.

*Case HRP-1:* Dirichlet BVP with regularly distributed data on each side, no noise.

— The RBF approach used 100 basis functions. The system was 3456×100, both the condition number and the residual depend on the parameters

**q**_{m}; they were Cn=1.8×10^{7}and Res=21.8, respectively, and no regularization has been performed.— FE2 approach with 25 elements along each edge with a total of 4048 data points employing the LSQR algorithms limited to 10 000 iterations. The number of collocation points was 45 000, the size of the system was 184 048×140 625, Cn was not estimated and Res=15.5 after all iterations.

— FE3 with five elements placed along each edge was solved by direct method (4.4). The size of the system was 3876×2000; the solution quality was Cn=4.7×10

^{4}, Res=25.3.

The results are summarized in figure 2*a*, where the absolute errors are shown for three profiles specified earlier. Although FE2 shows better accuracy in temperature recovery, it should be noted that the other two approaches also have satisfactory accuracy with considerably less computational recourses. Statistical characteristics are shown in table 1.

*Case HRP-2:* Inverse BVP with 930 zero data points on sides *Ω*, no data points on the top side, 10 data points in the cross section *x*_{3}=0.9 and four data points on the line connecting the points (0.5,0.5,0) and (0.5,0.5,1).

— The RBF approach used 100 basis functions. The size of the system was 944×100; the solution quality depends on the parameters

**q**_{m}; for the solutions presented in figure 2*b*, they were Cn=1.6×10^{7}, Res=7.5, respectively; no regularization was performed.— FE2 approach with seven elements along each edge was solved by direct method (4.4). The size of the system was 4472×3087, Cn=2.8×10

^{4}, Res=140, no regularization.— FE3 with five elements placed along each edge was solved by direct method (4.4). The size of the system was 3644×2000, Cn=4.4×10

^{4}, Res=16.1, no regularization.

Some results are summarized in figure 2*b*, where the absolute errors are shown for three profiles specified earlier. Statistical characteristics are shown in table 2. The FE methods require similar computational resources to produce errors of the same order and do not require regularization. The RBF approach is much more sensitive to the choice of parameters **q**_{m} than in the Dirichlet BVP, and therefore the errors vary; they may be of the same order as FE but may be large especially on the top of the domain where no data were specified. The RBF approach is also quite sensitive to the number of basis functions. The condition number often exceeded 10^{8} if more than 120 functions were used to fit the same dataset, which indicates poor stability of the solutions and, therefore, necessitates regularization. In such cases, the solutions were found in accordance with equation (4.5).

In order to model experimental errors, tests with the relative ±5% noise have also been performed. The FE approaches are not very sensitive to such noise because of low condition numbers, while RBF produces less stable results.

*Case HRP-3:* Interpolation problem with eight zero data points in the vertices of the cube and 100 randomly distributed data points inside the no noise area.

— RBF demonstrates reasonable accuracy if

*M*=25–50, i.e. from a quarter to a half of the data, Cn=10^{3}–10^{7}, Res=5–50; errors in temperature reconstruction are small for profiles 1 and 2; on the top of the cube they are more pronounced, reaching up to 20^{°}for particular sets of parameters**q**_{m}; one of the solutions is illustrated in figure 2*c*for*M*=40.— FE2 approach with seven elements along each edge was solved by direct method (4.4) without regularization. The size of the system was 3636×3087, Cn=7.1×10

^{4}and Res=70.1. figure 3 shows ideal and reconstructed temperature distributions on the top of the cube and in the cross section passing vertices (1,0,1), (0,1,1) and (1,1,0).— FE3 with five elements placed along each edge was solved by direct method (4.4) without regularization. The size of the system is 2808×2000, Cn=1.7×10

^{4}and Res=7.9.

Absolute errors are shown in figure 2*c*, and their statistical characteristics are summarized in table 3 for all three approaches. However, the results for RBF are for illustration only owing to sensitivity to parameters **q**_{m} and to the number of basis functions.

FE2 and FE3 produce results of similar quality; however, no regularization is needed for FE3 and therefore the residual is much smaller than that in FE2.

### (b) Isotropic three-dimensional elastostatics

Leonov (1981) presented the solution for a rigid circular punch of unit radius applied to the boundary of an isotropic upper half-space without friction in the form
5.3
This solution is further used to compile synthetic data. Namely, the dataset is formed as the ideal solution with additional noise at *N* points
5.4
where *N*_{e} is the number of points where the dilation has been measured, and *N*_{d} is the number of points where displacements have been monitored, *η*_{m} are random numbers uniformly distributed within the interval (−1,1) and *ε* is the level of noise (±5% in examples below). It is assumed that *N*≤*N*_{e}+*N*_{d} because, in some points, both dilations and displacements may be measured.

When selecting the computational domain, one should note that, although displacements are continuous everywhere, the strain/stress components are discontinuous across the circumference of unit radius representing the punch edge. Therefore, in simulations, the boundary is excluded from the computational domain because the technique assumes the continuity of basis functions. On the other hand, in order to check the effectiveness of reconstruction, the domain should not be selected far from the location where the punch acts and where the displacement field does not vary much. Here, the domain is a cube of unit volume situated in the upper half-space near the punch applied over the circle ; namely, *Ω*={0≤*x*_{1}≤1,0≤*x*_{2}≤1,0.1≤*x*_{3}≤1.1}.

In the tests described later, the reconstruction undergoes two stages; the dilation is reconstructed first and then the harmonic functions *φ*_{k} entering into equation (2.13) are found. For conciseness, the results for reconstructed displacements are reported for the vertical component of displacements only.

*Case ERP-1:* Dirichlet BVP with regularly distributed data on each side, no noise.

— The RBF approach used 100 basis functions. The size of the system was 1176×100, Cn=10

^{6}×10^{7}, Res<0.5 in both problems, no regularization applied.— FE2 approach with 31 elements along each edge with a total of 5766 data points employing the LSQR algorithms. The number of collocation points was 45 000, the system was 351 726×268 119, Res=0.33 (for dilation) and Res=0.43 (for harmonic functions,

*φ*_{3}) after 10 000 iterations. Deviation of approximate solution for vertical displacement from*u*_{3}data was 1.4×10^{−5}(calculated as the L_{2}norm of the difference of data and approximation over the number of data).— FE3 with 3×5×7 elements placed along the coordinate axes respectively was solved by direct method (4.4). The size of the system was 3372×1680 (1176 data), Cn=3.5×10

^{4}, the residuals were Res=0.61 (for dilation reconstruction) and Res=0.20 (for the reconstruction of*φ*_{3}). Deviation of the approximate solution from*u*_{3}data was 1.4×10^{−4}.

The results of calculations are presented in figure 4 for two profiles. Profile A is along the diagonal on the bottom of the cube *Ω* connecting points (0,0,0) with (1,1,0); profile B is along the vertical punch axis from point (0,0,0) to point (0,0,1). Ideal vertical displacements vary monotonically in the range from −0.98 to −0.46, and, therefore, relative errors give proper illustration of accuracy of the results; it is better than 3 per cent along profile A and better than 1 per cent along profile B in all three methods.

*Case ERP-2:* Inverse BVP with 1030 data points regularly distributed over the sides of *Ω* apart from the area located above the punch on the bottom side of *Ω* where no data are specified. The aim of this modelling is the reconstruction of displacements in this area (where the ideal *u*_{3} is close to unity). The lack of internal data is apparently the worst-case scenario and therefore it is demanding to expect perfect reconstruction in this case.

— The RBF approach with 30–150 basis functions has failed to produce reasonable results by both solution methods, (4.4) and (4.5).

— FE2 with eight elements along each side has produced errors up to 24.7% with the mean of 7.4% and standard deviation of 7.9%. The parameters of the solution were Cn=2.010

^{4}(no regularization) 6406×4608, Res=2.15 (for dilation) and Res=0.91 (for harmonic function,*φ*_{3}). Deviation of the approximate solution from the*u*_{3}data was 1.8×10^{−5}.— FE3 with four elements placed along the edges was solved by direct method (4.4). The size of the system is 2326×1024, the condition number was Cn=2.1×10

^{4}and the residuals were Res=0.29 (for dilation reconstruction) and Res=0.11 (for the reconstruction of*φ*_{3}). Deviation of the approximate solution from the*u*_{3}data was 8.3×10^{−5}.

The results of reconstruction are shown in figure 5 for the finite element approaches. The figure presents the contour maps of relative errors in the displacement reconstruction under the punch (these can also be considered as reconstructed maps of vertical displacements because the ideal displacements are close to unity). The maximum error in FE3 is below 15 per cent when compared with that of 24.7 per cent in FE2. In both cases, the maximum vertical displacements have been overestimated.

*Case ERP-3*: Interpolation problem with randomly distributed data, the relative noise of ±5%.

— RBF provides reasonable accuracy if

*M*=25–50, i.e. from a quarter to a half of the data, Cn=10^{3}−8×10^{7}, typical residuals are 0.1–0.2 for dilation and 0.2–0.3 for*φ*_{3}reconstruction; errors strongly depend on the number of basis functions and particular sets of parameters**q**_{m}; typical reconstruction is shown in figure 6 for*M*=40 (maximum error was less than 25%).— FE2 with four elements along each edge was solved by direct method (4.4). The size of the system is 684×576. For the example presented in figure 6, the condition number was Cn=4.8×10

^{3}, and the residual for*φ*_{3}reconstruction was Res=2.16 (no regularization).— FE3 with four elements placed along each edge was solved by direct method (4.4). The size of the system was 1404×1024. For the example presented in figure 6, the condition number was Cn=6.4×10

^{4}, and the residuals were Res=0.1 for dilation and Res=0.19 for*φ*_{3}reconstruction. Deviation from the*u*_{3}data was 10^{−3}.

It should be noted that FE3 showed stable behaviour in all tests; it was quite accurate even if two blocks were placed along every edge.

## 6. Conclusions

Numerical experiments have been conducted by three approaches with three different data types modelling direct, inverse and interpolation problems. All these approaches are based on the Trefftz method; they provide smoothness of the reconstructed three-dimensional harmonic functions, although the finite element approaches FE2 and FE3 satisfy the continuity conditions approximately in collocation points located on the interfaces between adjacent blocks.

All approaches are capable of producing satisfactory results for all direct boundary value problems investigated in this paper. FE2 shows high accuracy when many blocks are used in approximation and when the LSQR algorithm is applied for solving the system.

For inverse boundary value problems, the finite approaches produce results of the same quality; RBF for the heat reconstruction problems was less stable, which limits its application to noisy data and led to large errors in displacement reconstruction in ERP-2.

For interpolation problems, the results produced by RBF are sensitive to the choice of parameters **q**_{m} and to the number of basis functions, and, therefore, this approach is less reliable. It has been found that FE3 usually leads to a system with a lower condition number than that resulting from FE2. Therefore, FE3 is less sensitive to experimental errors than FE2 and, therefore, better approximations can be achieved without regularization.

It should be emphasized that the reconstruction of harmonic functions in three dimensions from finite datasets of measurements is not a unique operation. However, the results of such reconstructions obey all physical relationships imposed on the sought solutions and, therefore, in contrast to conventional interpolation, they can be used for further analysis; for example, for reconstruction of heat fluxes in temperature problems, or strain and stress in elastostatics.

## Acknowledgements

The authors are grateful to the EPSRC for partial support of this work through research grant EP/E032494/1.

## Footnotes

- Received September 7, 2009.
- Accepted January 5, 2010.

- © 2010 The Royal Society