Add basic setup to solve PDEs by collocation (#36)

* add basic setup to solve PDEs by collocation

* format

* remove commented code (not needed anymore)

* fix inundation

* fix docstring of system_matrix

* use WendlandKernel for better stability

* fix typo in docstring

* differential operators and pdes only work with a RadialSymmetricKernel

* remove unbounded type parameter

* remove KernelInterpolation from docs/Project.toml

Project.toml

@@ -4,13 +4,15 @@ authors = ["Joshua Lampert <[email protected]> and contributors"]
 version = "1.0.0-DEV"
 
 [deps]
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 TypedPolynomials = "afbbf031-7a57-5f58-a1b9-b774a0fad08d"
 
 [compat]
+ForwardDiff = "0.10.36"
 LinearAlgebra = "1"
 RecipesBase = "1.1"
 SpecialFunctions = "2"

docs/make.jl

@@ -16,6 +16,7 @@ makedocs(;
                                   edit_link = "main",
                                   assets = String[]),
          pages = ["Home" => "index.md",
+             "Solving PDEs" => "pde.md",
              "Reference" => "ref.md",
              "License" => "license.md"])
 

docs/src/pde.md

@@ -0,0 +1,5 @@
+# Solving PDEs by collocation
+TODO:
+* Explain basic setup
+* Define custom differential operators and PDEs and solve them
+* AD vs analytic derivatives

examples/PDEs/poisson_2d.jl

@@ -0,0 +1,26 @@
+using KernelInterpolation
+using Plots
+
+# right-hand-side of Poisson equation
+f(x) = 5 / 4 * pi^2 * sin(pi * x[1]) * cos(pi * x[2] / 2)
+pde = PoissonEquation(f)
+
+# analytical solution of equation
+u(x) = sin(pi * x[1]) * cos(pi * x[2] / 2)
+
+n = 10
+nodeset_inner = homogeneous_hypercube(n, (0.1, 0.1), (0.9, 0.9); dim = 2)
+n_boundary = 3
+nodeset_boundary = homogeneous_hypercube_boundary(n_boundary; dim = 2)
+values = f.(nodeset_inner)
+# Dirichlet boundary condition (here taken from analytical solution)
+g(x) = u(x)
+values_boundary = g.(nodeset_boundary)
+
+kernel = WendlandKernel{2}(3, shape_parameter = 0.3)
+itp = solve(pde, nodeset_inner, nodeset_boundary, values_boundary, kernel)
+
+many_nodes = homogeneous_hypercube(20; dim = 2)
+
+plot(many_nodes, itp)
+plot!(many_nodes, u)

src/KernelInterpolation.jl

@@ -1,13 +1,16 @@
 module KernelInterpolation
 
-using LinearAlgebra: norm, Symmetric
+using LinearAlgebra: norm, Symmetric, tr
 using RecipesBase: RecipesBase, @recipe, @series
 using SpecialFunctions: besselk, loggamma
 using StaticArrays: StaticArrays, MVector
 using TypedPolynomials: Variable, monomials, degree
+using ForwardDiff: ForwardDiff
 
 include("kernels/kernels.jl")
 include("nodes.jl")
+include("differential_operators.jl")
+include("pdes.jl")
 include("interpolation.jl")
 include("visualization.jl")
 include("util.jl")
@@ -19,12 +22,14 @@ export GaussKernel, MultiquadricKernel, InverseMultiquadricKernel,
        Matern52Kernel, Matern72Kernel, RieszKernel,
        TransformationKernel, ProductKernel, SumKernel
 export phi, Phi, order
+export Laplacian
+export PoissonEquation
 export NodeSet, separation_distance, dim, values_along_dim, random_hypercube,
        random_hypercube_boundary, homogeneous_hypercube, homogeneous_hypercube_boundary,
        random_hypersphere, random_hypersphere_boundary
 export interpolation_kernel, nodeset, coefficients, kernel_coefficients,
        polynomial_coefficients, polynomial_basis, polyvars, system_matrix,
-       interpolate, kernel_inner_product, kernel_norm
+       interpolate, solve, kernel_inner_product, kernel_norm
 export examples_dir, get_examples, default_example, include_example
 
 end

src/differential_operators.jl

@@ -0,0 +1,35 @@
+abstract type AbstractDifferentialOperator end
+
+function (D::AbstractDifferentialOperator)(kernel::RadialSymmetricKernel, x, y)
+    @assert length(x) == length(y)
+    return save_call(D, kernel, x .- y)
+end
+
+# Workaround to avoid evaluating the derivative at zeros to allow automatic differentiation,
+# see https://github.com/JuliaDiff/ForwardDiff.jl/issues/303
+# the same issue appears with Zygote.jl
+function save_call(D::AbstractDifferentialOperator, kernel::RadialSymmetricKernel, x)
+    if all(iszero, x)
+        x[1] = eps(typeof(x[1]))
+    end
+    # x .+= eps(typeof(x[1]))
+    return D(kernel, x)
+end
+
+"""
+    Laplacian()
+
+The Laplacian operator. It can be called with an `RadialSymmetricKernel` and points
+`x` and `y` to evaluate the Laplacian of the `kernel` at `x - y`.
+"""
+struct Laplacian <: AbstractDifferentialOperator
+end
+
+function Base.show(io::IO, ::Laplacian)
+    print(io, "Δ")
+end
+
+function (::Laplacian)(kernel::RadialSymmetricKernel, x)
+    H = ForwardDiff.hessian(x -> Phi(kernel, x), x)
+    return tr(H)
+end

src/interpolation.jl

@@ -8,7 +8,7 @@ Return the kernel from an interpolation object.
 interpolation_kernel(itp::AbstractInterpolation) = itp.kernel
 
 """
-    nodeset(itp)
+	nodeset(itp)
 
 Return the node set from an interpolation object.
 """
@@ -34,7 +34,7 @@ struct Interpolation{Kernel, Dim, RealT, A, Monomials, PolyVars} <:
     kernel::Kernel
     nodeset::NodeSet{Dim, RealT}
     c::Vector{RealT}
-    symmetric_system_matrix::A
+    system_matrix::A
     ps::Monomials
     xx::PolyVars
 end
@@ -110,20 +110,15 @@ order(itp::Interpolation) = maximum(degree.(itp.ps), init = -1) + 1
 @doc raw"""
     system_matrix(itp::Interpolation)
 
-Return the system matrix, i.e., the matrix
+Return the system matrix, i.e., the matrix ``A`` in the linear system
 ```math
-\begin{pmatrix}
-A & P \\
-P^T & 0
-\end{pmatrix},
+    Ac = f,
 ```
-where ``A\in\mathbb{R}^{n\times n}`` is the matrix with entries
-``a_{ij} = K(x_i, x_j)`` for the kernel function `K` and nodes `x_i`
-and ``P\in\mathbb{R}^{n\times q}`` is the matrix with entries
-``p_{ij} = p_j(x_i)``, where ``p_j`` is the ``j``-th multivariate monomial
-of the space of polynomials up to degree ``m``.
+where ``c`` are the coefficients of the kernel interpolant and ``f`` the vector
+of known values. The exact form of ``A`` differs depending on  whether classical interpolation
+or collocation is used.
 """
-system_matrix(itp::Interpolation) = itp.symmetric_system_matrix
+system_matrix(itp::Interpolation) = itp.system_matrix
 
 @doc raw"""
     interpolate(nodeset, values, kernel = GaussKernel{dim(nodeset)}(), m = order(kernel))
@@ -142,7 +137,7 @@ maximum degree of `m - 1`. If `m = 0`, no polynomial is added. The additional co
 are enforced. Returns an [`Interpolation`](@ref) object.
 """
 function interpolate(nodeset::NodeSet{Dim, RealT}, values::Vector{RealT},
-                     kernel = GaussKernel{dim(nodeset)}(),
+                     kernel = GaussKernel{Dim}(),
                      m = order(kernel)) where {Dim, RealT}
     @assert dim(kernel) == Dim
     n = length(nodeset)
@@ -171,6 +166,48 @@ function interpolate(nodeset::NodeSet{Dim, RealT}, values::Vector{RealT},
     return Interpolation(kernel, nodeset, c, symmetric_system_matrix, ps, xx)
 end
 
+"""
+    solve(pde, nodeset_inner, nodeset_boundary, values_boundary, kernel = GaussKernel{dim(nodeset_inner)})
+
+Solve a partial differential equation `pde` with Dirichlet boundary conditions by non-symmetric collocation
+(Kansa method) using the kernel `kernel`. The `nodeset_inner` are the nodes in the domain and `nodeset_boundary`
+are the nodes on the boundary. The `values_boundary` are the values of the boundary condition at the nodes given
+by Dirichlet boundary conditions. Returns an [`Interpolation`](@ref) object.
+"""
+function solve(pde, nodeset_inner::NodeSet{Dim, RealT},
+               nodeset_boundary::NodeSet{Dim, RealT},
+               values_boundary::Vector{RealT},
+               kernel = GaussKernel{Dim}()) where {Dim, RealT}
+    @assert dim(kernel) == Dim
+    n_i = length(nodeset_inner)
+    n_b = length(nodeset_boundary)
+    nodeset = merge(nodeset_inner, nodeset_boundary)
+    n = n_i + n_b
+
+    pde_matrix = Matrix{RealT}(undef, n_i, n)
+    for i in 1:n_i
+        for j in 1:n
+            pde_matrix[i, j] = pde(kernel, nodeset_inner[i], nodeset[j])
+        end
+    end
+    boundary_matrix = Matrix{RealT}(undef, n_b, n)
+    for i in 1:n_b
+        for j in 1:n
+            # Dirichlet boundary condition
+            boundary_matrix[i, j] = kernel(nodeset_boundary[i], nodeset[j])
+        end
+    end
+    system_matrix = [pde_matrix
+                     boundary_matrix]
+    b = [rhs(pde, nodeset_inner); values_boundary]
+    c = system_matrix \ b
+
+    # Do not support additional polynomial basis for now
+    xx = polyvars(Dim)
+    ps = monomials(xx, 0:-1)
+    return Interpolation(kernel, nodeset, c, system_matrix, ps, xx)
+end
+
 # Evaluate interpolant
 function (itp::Interpolation)(x)
     s = 0
@@ -196,6 +233,28 @@ function (itp::Interpolation)(x::RealT) where {RealT <: Real}
     return itp([x])
 end
 
+function (diff_op::AbstractDifferentialOperator)(itp::Interpolation, x)
+    kernel = interpolation_kernel(itp)
+    xs = nodeset(itp)
+    c = kernel_coefficients(itp)
+    s = 0
+    for j in 1:length(c)
+        s += c[j] * diff_op(kernel, x, xs[j])
+    end
+    return s
+end
+
+function (pde::AbstractPDE)(itp::Interpolation, x)
+    kernel = interpolation_kernel(itp)
+    xs = nodeset(itp)
+    c = kernel_coefficients(itp)
+    s = 0
+    for j in 1:length(c)
+        s += c[j] * pde(kernel, x, xs[j])
+    end
+    return s
+end
+
 # TODO: Does this also make sense for conditionally positive definite kernels?
 @doc raw"""
     kernel_inner_product(itp1, itp2)