Interpolation for space scattered data by bicubic polynomial natural splines

Bicubic splines interpolate for space scattered data such that the integral of square of partial derivative of two orders to x and to y for the interpolating function is minimal (with natural boundary conditions). The solution is constructed as the sum of a bilinear polynomial and piecewise bicubic polynomials by Hilbert space spline function methods. Its coefficients can be decided by a linear system. The coefficient matrix is so symmetry that the LDLT method can be successed. Results are very simple and can be achieved easily in computer programs.