update triangle.jl to remove ApproxFun dependency

Author: MikaelSlevinsky

examples/sphere.jl

@@ -8,14 +8,12 @@
 # To verify, we sample the function on a 5×9 tensor product grid
 # at equispaced points-in-angle defined by:
 #
-#   θₙ = (n+1/2)π/N, for 0 ≤ n < N,
-#
-# and
+#   θₙ = (n+1/2)π/N, for 0 ≤ n < N, and
 #
 #   φₘ = 2π m/M, for 0 ≤ m < M;
 #
 # we convert the function samples to Fourier coefficients using
-# `FastTransforms.plan_analysis`; and finally, we transform
+# `plan_sph_analysis`; and finally, we transform
 # the Fourier coefficients to spherical harmonic coefficients using
 # `plan_sph2fourier`.
 #
@@ -33,7 +31,7 @@ function threshold!(A::AbstractArray, ϵ)
     A
 end
 
-using FastTransforms
+using FastTransforms, LinearAlgebra
 
 # The colatitudinal grid (mod π):
 N = 5

examples/triangle.jl

@@ -1,47 +1,89 @@
 #############
-# This demonstrates the Triangle harmonic transform and inverse transform,
-# explaining precisely the normalization and points
+# In this example, we sample a bivariate function:
 #
-# Note we use the duffy map
-# x == (s+1)/2
-# y== (t+1)/2*(1-(s+1)/2)
+#   f(x,y) = 1/(1+x^2+y^2),
+#
+# on the reference triangle with vertices (0,0), (0,1), and (1,0) and analyze it
+# in a Proriol series. Then, we find Proriol series for each component of its
+# gradient by term-by-term differentiation of our expansion, and we compare them
+# with the true Proriol series by sampling an exact expression for the gradient.
+#
+# We analyze f(x,y) on an N×M mapped tensor product grid defined by:
+#
+#   x = (1+u)/2, and y = (1-u)*(1+v)/4, where:
+#
+#   uₙ = cos[(n+1/2)π/N], for 0 ≤ n < N, and
+#
+#   vₘ = cos[(m+1/2)π/M], for 0 ≤ m < M;
+#
+# we convert the function samples to mapped Chebyshev² coefficients using
+# `plan_tri_analysis`; and finally, we transform the mapped Chebyshev²
+# coefficients to Proriol coefficients using `plan_tri2cheb`.
+#
+# For the storage pattern of the arrays, please consult the documentation.
 #############
 
+using FastTransforms, LinearAlgebra
 
-using ApproxFun, FastTransforms
+f = (x,y) -> 1/(1+x^2+y^2)
+fx = (x,y) -> -2x/(1+x^2+y^2)^2
+fy = (x,y) -> -2y/(1+x^2+y^2)^2
 
-jacobinorm(n,a,b) = if n ≠ 0
-        sqrt((2n+a+b+1))*exp((lgamma(n+a+b+1)+lgamma(n+1)-log(2)*(a+b+1)-lgamma(n+a+1)-lgamma(n+b+1))/2)
-    else
-        sqrt(exp(lgamma(a+b+2)-log(2)*(a+b+1)-lgamma(a+1)-lgamma(b+1)))
-    end
-njacobip(n,a,b,x) = jacobinorm(n,a,b) * jacobip(n,a,b,x)
+N = 10
+M = N
+
+α, β, γ = 0, 0, 0
 
-P = (ℓ,m,x,y) -> (2*(1-x))^m*njacobip(ℓ-m,2m,0,2x-1)*njacobip(m,-0.5,-0.5,2y/(1-x)-1)
+P = plan_tri2cheb(Float64, N, α, β, γ)
+Px = plan_tri2cheb(Float64, N, α+1, β, γ+1)
+Py = plan_tri2cheb(Float64, N, α, β+1, γ+1)
+PA = plan_tri_analysis(Float64, N, M)
 
+u = [sinpi((N-2n-1)/(2N)) for n in 0:N-1]
+v = [sinpi((M-2m-1)/(2M)) for m in 0:M-1]
 
+# Instead of using the u, v grid, we use one with more accuracy near the origin.
+x = [sinpi((2N-2n-1)/(4N))^2 for n in 0:N-1]
+w = [sinpi((2M-2m-1)/(4M))^2 for m in 0:M-1]
 
-p_T = chebyshevpoints(40)
-f = (x,y) -> exp(x + cos(y))
-    f̃ = (s,t) -> f((s+1)/2, (t+1)/2*(1-(s+1)/2))
+(1 .+ u)./2 ≈ x
+(1 .- u).*(1 .+ v')/4 ≈ reverse(x).*w'
 
-    F = f̃.(p_T, p_T')
-        for j = 1:size(F,2)
-            F[:,j] = chebyshevtransform(F[:,j])
-        end
-        for k = 1:size(F,1)
-            F[k,:] = chebyshevtransform(F[k,:])
-        end
+# On the mapped tensor product grid, our function samples are:
+F = [f(x[n+1], x[N-n]*w[m+1]) for n in 0:N-1, m in 0:M-1]
 
-    F̌ = cheb2tri(F, 0.0, -0.5, -0.5)
+# Its Proriol-(α,β,γ) coefficients are:
+U = P\(PA*F)
 
+# Similarly, our function's gradient samples are:
+Fx = [fx(x[n+1], x[N-n]*w[m+1]) for n in 0:N-1, m in 0:M-1]
+Fy = [fy(x[n+1], x[N-n]*w[m+1]) for n in 0:N-1, m in 0:M-1]
 
-f̃ = function(x,y)
-        ret = 0.0
-        for j=1:size(F,2), k=1:size(F,1)-j+1
-            ret += F̌[k,j] * P(k+j-2,j-1,x,y)
-        end
-        ret
+# For the partial derivative with respect to x, Olver et al.
+# derive simple expressions for the representation of this component
+# using a Proriol-(α+1,β,γ+1) series. For the partial derivative with respect
+# to y, the analogous formulae result in a Proriol-(α,β+1,γ+1) series.
+# These expressions are adapted from https://arxiv.org/abs/1902.04863.
+Gx = zeros(Float64, N, M)
+for m = 0:M-2
+    for n = 0:N-2
+        cf1 = m == 0 ? sqrt((n+1)*(n+2m+α+β+γ+3)/(2m+β+γ+2)*(m+γ+1)*8) : sqrt((n+1)*(n+2m+α+β+γ+3)/(2m+β+γ+1)*(m+β+γ+1)/(2m+β+γ+2)*(m+γ+1)*8)
+        cf2 = sqrt((n+α+1)*(m+1)/(2m+β+γ+2)*(m+β+1)/(2m+β+γ+3)*(n+2m+β+γ+3)*8)
+        Gx[n+1, m+1] = cf1*U[n+2, m+1] + cf2*U[n+1, m+2]
     end
+end
+Ux = Px\(PA*Fx)
+
+Gy = zeros(Float64, N, M)
+for m = 0:M-2
+    for n = 0:N-2
+        Gy[n+1, m+1] = 4*sqrt((m+1)*(m+β+γ+2))*U[n+1, m+2]
+    end
+end
+Uy = Py\(PA*Fy)
+
+# The 2-norm relative error in differentiating the Proriol series
+# for f(x,y) term-by-term and its sampled gradient is:
+hypot(norm(Ux-Gx), norm(Uy-Gy))/hypot(norm(Ux), norm(Uy))
 
-f̃(0.1,0.2) ≈ f(0.1,0.2)
+# This error can be improved upon by increasing N and M.