|
Here we interpolate a function $f(x,y)$ over the two-dimensional unit cube from $(0,0)$ to $(1,1)$, using the multiquadric Radial Basis Function (RBF).
The function to be interpolated is
\begin{equation} f(x,y) = \exp(-9(x-\frac{1}{2})^2 - 9(y-\frac{1}{4})^2), \end{equation}
and the multiquadric RBF $K(r)$ is given by
\begin{equation} K(r) = \sqrt{r^2 + c^2}. \end{equation}
We create a set of basis functions around a set of $N$ particles with positions $\mathbf{x}_i$. In this case the radial coordinate around each point is $r = ||\mathbf{x}_i - \mathbf{x}||$. The function $f$ is approximated using the sum of these basis functions, with the addition of a constant factor $\beta$
\begin{equation} \overline{f}(\mathbf{x}) = \beta + \sum_i \alpha_j K(||\mathbf{x}_i-\mathbf{x}||). \end{equation}
We exactly interpolate the function at $\mathbf{x}_i$ (i.e. set $\overline{f}(\mathbf{x}_i)=f(\mathbf{x}_i)$), leading to
\begin{equation} f(\mathbf{x}_i) = \beta + \sum_j \alpha_j K(||\mathbf{x}_j-\mathbf{x}_i||). \end{equation}
Note that the sum $j$ is over the same particle set.
We also need one additional constraint to find $\beta$
\begin{equation} 0 = \sum_j \alpha_j. \end{equation}
We rewrite the two previous equations as a linear algebra expression
\begin{equation} \mathbf{\Phi} = \begin{bmatrix} \mathbf{G}&\mathbf{P} \\ \mathbf{P}^T & \mathbf{0} \end{bmatrix} \begin{bmatrix} \mathbf{\alpha} \\ \mathbf{\beta} \end{bmatrix} = \mathbf{W} \mathbf{\gamma}, \end{equation}
where $\mathbf{\Phi} = [f(\mathbf{x}_1),f(\mathbf{x}_2),...,f(\mathbf{x}_N),0]$, $\mathbf{G}$ is an $N \times N$ matrix with elements $G_{ij} = K(||\mathbf{x}_j-\mathbf{x}_i||)$ and $\mathbf{P}$ is a $N \times 1$ vector of ones.
The basis function coefficients $\mathbf{\gamma}$ are found by solving this equation. To do this, we use the iterative solvers found in the Eigen package and Aboria's ability to wrap C++ function objects as Eigen matricies.
#include "Aboria.h" #include <random> using namespace Aboria; int main() { auto funct = [](const double x, const double y) { return std::exp(-9 * std::pow(x - 0.5, 2) - 9 * std::pow(y - 0.25, 2)); // return x; }; ABORIA_VARIABLE(alpha, double, "alpha value") ABORIA_VARIABLE(interpolated, double, "interpolated value") ABORIA_VARIABLE(constant2, double, "c2 value") typedef Particles<std::tuple<alpha, constant2, interpolated>, 2, std::vector, SearchMethod> ParticlesType; typedef position_d<2> position; typedef typename ParticlesType::const_reference const_particle_reference; typedef Eigen::Matrix<double, Eigen::Dynamic, 1> vector_type; typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> matrix_type; ParticlesType knots; ParticlesType augment; ParticlesType test; const double c = 0.5; const int N = 1000; const int max_iter = 100; const int restart = 101; typename ParticlesType::value_type p; std::default_random_engine generator; std::uniform_real_distribution<double> distribution(0.0, 1.0); for (int i = 0; i < N; ++i) { get<position>(p) = vdouble2(distribution(generator), distribution(generator)); get<constant2>(p) = std::pow(c, 2); knots.push_back(p); get<position>(p) = vdouble2(distribution(generator), distribution(generator)); get<constant2>(p) = std::pow(c, 2); test.push_back(p); } // knots.init_neighbour_search(min,max,periodic); augment.push_back(p); auto kernel = [](const_particle_reference a, const_particle_reference b) { return std::sqrt((get<position>(b) - get<position>(a)).squaredNorm() + get<constant2>(b)); }; auto one = [](const_particle_reference a, const_particle_reference b) { return 1.0; }; auto G = create_dense_operator(knots, knots, kernel); auto P = create_dense_operator(knots, augment, one); auto Pt = create_dense_operator(augment, knots, one); auto Zero = create_zero_operator(augment, augment); auto W = create_block_operator<2, 2>(G, P, Pt, Zero); auto G_test = create_dense_operator(test, knots, kernel); auto W_test = create_block_operator<2, 2>(G_test, P, Pt, Zero); vector_type phi(N + 1), gamma(N + 1); for (size_t i = 0; i < knots.size(); ++i) { const double x = get<position>(knots[i])[0]; const double y = get<position>(knots[i])[1]; phi[i] = funct(x, y); } phi[knots.size()] = 0; matrix_type W_matrix(N + 1, N + 1); W.assemble(W_matrix); gamma = W_matrix.ldlt().solve(phi); Eigen::GMRES<matrix_type> gmres; gmres.setMaxIterations(max_iter); gmres.set_restart(restart); gmres.compute(W_matrix); gamma = gmres.solve(phi); std::cout << "GMRES: #iterations: " << gmres.iterations() << ", estimated error: " << gmres.error() << std::endl; phi = W * gamma; double rms_error = 0; double scale = 0; for (size_t i = 0; i < knots.size(); ++i) { const double x = get<position>(knots[i])[0]; const double y = get<position>(knots[i])[1]; const double truth = funct(x, y); const double eval_value = phi[i]; rms_error += std::pow(eval_value - truth, 2); scale += std::pow(truth, 2); // TS_ASSERT_DELTA(eval_value,truth,2e-3); } std::cout << "rms_error for global support, at centers = " << std::sqrt(rms_error / scale) << std::endl; TS_ASSERT_LESS_THAN(std::sqrt(rms_error / scale), 1e-6); phi = W_test * gamma; rms_error = 0; scale = 0; for (size_t i = 0; i < test.size(); ++i) { const double x = get<position>(test[i])[0]; const double y = get<position>(test[i])[1]; const double truth = funct(x, y); const double eval_value = phi[i]; rms_error += std::pow(eval_value - truth, 2); scale += std::pow(truth, 2); // TS_ASSERT_DELTA(eval_value,truth,2e-3); } std::cout << "rms_error for global support, away from centers = " << std::sqrt(rms_error / scale) << std::endl; TS_ASSERT_LESS_THAN(std::sqrt(rms_error / scale), 1e-5); }
#include "Aboria.h" using namespace Aboria; int main() { auto funct = [](const double x, const double y) { return std::cos(4*x+4*y); }; auto laplace_funct = [](const double x, const double y) { return -32*std::cos(4*x+4*y); }; ABORIA_VARIABLE(boundary,uint8_t,"is boundary knot") ABORIA_VARIABLE(interpolated,double,"interpolated value") ABORIA_VARIABLE(constant2,double,"c2 value") ABORIA_VARIABLE(alpha,double,"alpha value") typedef Particles<std::tuple<alpha,boundary,constant2,interpolated>,2> ParticlesType; typedef position_d<2> position; typedef typename ParticlesType::const_reference const_particle_reference; typedef Eigen::Map<Eigen::Matrix<double,Eigen::Dynamic,1>> map_type; typedef Eigen::Matrix<double,Eigen::Dynamic,1> vector_type; ParticlesType knots,augment; const double c = 0.5; const int max_iter = 100; const int restart = 100; const int nx = 7; constexpr int N = (nx+1)*(nx+1); const double delta = 1.0/nx; typename ParticlesType::value_type p; for (int i=0; i<=nx; ++i) { for (int j=0; j<=nx; ++j) { get<position>(p) = vdouble2(i*delta,j*delta); if ((i==0)||(i==nx)||(j==0)||(j==nx)) { get<boundary>(p) = true; } else { get<boundary>(p) = false; } get<constant2>(p) = std::pow(c,2); knots.push_back(p); } } augment.push_back(p); auto kernel = []( const_particle_reference a, const_particle_reference b) { const vdouble2 dx = get<position>(b) - get<position>(a); return std::exp(-dx.squaredNorm()/get<constant2>(b)); }; auto laplace_kernel = []( const_particle_reference a, const_particle_reference b) { const vdouble2 dx = get<position>(b) - get<position>(a); return 4.0*(dx.squaredNorm()-get<constant2>(b))* std::exp(-dx.squaredNorm()/get<constant2>(b))/ (get<constant2>(b)*get<constant2>(b)); }; auto G = create_dense_operator(knots,knots, [kernel,laplace_kernel]( const_particle_reference a, const_particle_reference b) { if (get<boundary>(a)) { return kernel(a,b); } else { return laplace_kernel(a,b); } }); auto P = create_dense_operator(knots,augment, []( const_particle_reference a, const_particle_reference b) { if (get<boundary>(a)) { return 1.0; } else { return 0.0; } }); auto Pt = create_dense_operator(augment,knots, []( const_particle_reference a, const_particle_reference b) { if (get<boundary>(b)) { return 1.0; } else { return 0.0; } }); auto Zero = create_zero_operator(augment,augment); auto W = create_block_operator<2,2>(G, P, Pt,Zero); vector_type phi(N+1), gamma(N+1); for (int i=0; i<N; ++i) { const double x = get<position>(knots[i])[0]; const double y = get<position>(knots[i])[1]; if (get<boundary>(knots[i])) { phi[i] = funct(x,y); } else { phi[i] = laplace_funct(x,y); } } phi[N] = 0; std::cout << std::endl; Eigen::GMRES<decltype(W), Eigen::DiagonalPreconditioner<double>> gmres; gmres.set_restart(restart); gmres.setMaxIterations(max_iter); gmres.compute(W); gamma = gmres.solve(phi); std::cout << "GMRES: #iterations: " << gmres.iterations() << ", estimated error: " << gmres.error() << std::endl; phi = W*gamma; for (int i=0; i<N; ++i) { const double x = get<position>(knots[i])[0]; const double y = get<position>(knots[i])[1]; if (get<boundary>(knots[i])) { TS_ASSERT_DELTA(phi[i],funct(x,y),2e-3); } else { TS_ASSERT_DELTA(phi[i],laplace_funct(x,y),2e-3); } } TS_ASSERT_DELTA(phi[N],0,2e-3); map_type alpha_wrap(get<alpha>(knots).data(),N); map_type interp_wrap(get<interpolated>(knots).data(),N); alpha_wrap = gamma.segment<N>(0); vector_type beta = vector_type::Constant(N,gamma[N]); auto K = create_dense_operator(knots,knots,kernel); interp_wrap = K*alpha_wrap + beta; double rms_error = 0; double scale = 0; for (int i=0; i<N; ++i) { const double x = get<position>(knots[i])[0]; const double y = get<position>(knots[i])[1]; const double truth = funct(x,y); const double eval_value = get<interpolated>(knots[i]); rms_error += std::pow(eval_value-truth,2); scale += std::pow(truth,2); TS_ASSERT_DELTA(eval_value,truth,1e-2); } std::cout << "rms_error for global support, at centers = "<<std::sqrt(rms_error/scale)<<std::endl; TS_ASSERT_LESS_THAN(std::sqrt(rms_error/scale),1e-3); }