Add first moment inside interval (#8)

* add moment between -1 and 1

* expand on power law tests

src/power.jl

@@ -20,7 +20,11 @@ end
 
 function powerlawmoment(::Val{1}, α, λ, z::Real)
     T = promote_type(typeof(α), typeof(λ), typeof(z))
-    -sign(z)α*λ*beta(one(T)/2, λ+one(T)/2)*abs(z)^(α-1)*_₂F₁((1-α)/2, 1-α/2, 2+λ, 1/z^2)/(1+λ)
+    if -1 ≤ z ≤ 1
+        -2^(2-α)*λ*sqrt(convert(T,π))*gamma(α+2)*gamma(λ+one(T)/2)/(gamma(α/2)*gamma(α/2+λ+1))*((2*z^2*(α+λ)+1)*_₂F₁(-α/2,-λ-α/2,one(T)/2,z^2)+(z^2-1)*_₂F₁(-α/2,-λ-α/2,-one(T)/2,z^2))/(α*(α+1)*(α+2*λ+1)*z)+2*λ*z*gamma((α+1)/2)*gamma(λ+one(T)/2)/(gamma(λ+α/2+1))*_₂F₁(-α/2,-λ-α/2,one(T)/2,z^2)
+    else
+        -sign(z)α*λ*beta(one(T)/2, λ+one(T)/2)*abs(z)^(α-1)*_₂F₁((1-α)/2, 1-α/2, 2+λ, 1/z^2)/(1+λ)
+    end
 end
 
 

test/test_power.jl

@@ -1,37 +1,38 @@
 using SingularIntegrals, ClassicalOrthogonalPolynomials, Test
 using SingularIntegrals: PowerKernelPoint, powerlawmoment, powerlawrecurrence, RecurrenceArray
 
-
-# using HypergeometricFunctions
-
-
-# α,λ = 0.2,0.4
-# T = Float64
-
-# x = 10.0
-
-
 @testset "Weights" begin
-    @testset "on interval" begin
-        z = 0.2
-        x = axes(ChebyshevT(),1)
-        α = 0.1
-        L = abs.(z .- x') .^ α
-        λ = 1
-        L0 = powerlawmoment(Val(0), α, λ, z)
-        # L1 = powerlawmoment(Val(1), α, λ, z)
+    z = 10.0
+    x = axes(ChebyshevT(),1)
+    α = 0.1
+    L = abs.(z .- x') .^ α
+    # from Mathmatica
+    λ = 1
+    L0 = powerlawmoment(Val(0), α, λ, z)
+    L1 = powerlawmoment(Val(1), α, λ, z)
+    @test L* UltrasphericalWeight(λ) ≈ L0  ≈ 1.9772924292721128
+    @test L1 ≈ -0.009901716900034385
+    A, B, C = powerlawrecurrence(α, λ)
+    @test (A[2]z + B[2])*L1-C[2]L0 ≈ -0.00022324029766696007
+    r = RecurrenceArray(z, (A,B,C), [L0,L1])
+    @test r[5] ≈ -2.5742591209035326E-7 
+    @test r[1:10] ≈ (L * Weighted(Ultraspherical(λ)))[1,1:10]
 
-        @test L* UltrasphericalWeight(λ) ≈ L0  ≈ 1.407056672015012 # from Mathmatica
-    end
-    @testset "off interval" begin
+    @test L* LegendreWeight() ≈ 2.517472100701719
+    @test (L * Legendre())[5] ≈ -1.3328397976790363E-7 
+end
+
+@testset "Compare with Mathematica results" begin
+    @testset "z >> 1" begin
         z = 10.0
         x = axes(ChebyshevT(),1)
         α = 0.1
         L = abs.(z .- x') .^ α
+        # from Mathematica
         λ = 1
         L0 = powerlawmoment(Val(0), α, λ, z)
         L1 = powerlawmoment(Val(1), α, λ, z)
-        @test L* UltrasphericalWeight(λ) ≈ L0  ≈ 1.9772924292721128 # from Mathmatica
+        @test L* UltrasphericalWeight(λ) ≈ L0  ≈ 1.9772924292721128
         @test L1 ≈ -0.009901716900034385
         A, B, C = powerlawrecurrence(α, λ)
         @test (A[2]z + B[2])*L1-C[2]L0 ≈ -0.00022324029766696007
@@ -42,4 +43,49 @@ using SingularIntegrals: PowerKernelPoint, powerlawmoment, powerlawrecurrence, R
         @test L* LegendreWeight() ≈ 2.517472100701719
         @test (L * Legendre())[5] ≈ -1.3328397976790363E-7 
     end
-end
+    @testset "z > 1" begin
+        z = 1.1
+        x = axes(ChebyshevT(),1)
+        α = 0.1
+        L = abs.(z .- x') .^ α
+        # from Mathematica
+        λ = 1
+        L0 = powerlawmoment(Val(0), α, λ, z)
+        L1 = powerlawmoment(Val(1), α, λ, z)
+        @test L* UltrasphericalWeight(λ) ≈ L0  ≈ 1.56655191643602910
+        @test L1 ≈ -0.0850992609853987
+        A, B, C = powerlawrecurrence(α, λ)
+        r = RecurrenceArray(z, (A,B,C), [L0,L1])
+        @test r[5] ≈ -0.00381340800899034
+    end
+    @testset "z < -1" begin
+        z = -1.43
+        x = axes(ChebyshevT(),1)
+        α = 0.1
+        L = abs.(z .- x') .^ α
+        # from Mathematica
+        λ = 1
+        L0 = powerlawmoment(Val(0), α, λ, z)
+        L1 = powerlawmoment(Val(1), α, λ, z)
+        @test L* UltrasphericalWeight(λ) ≈ L0  ≈ 1.617769472203235
+        @test L1 ≈ 0.06180254575531147
+        A, B, C = powerlawrecurrence(α, λ)
+        r = RecurrenceArray(z, (A,B,C), [L0,L1])
+        @test r[5] ≈ -0.000807806617344142
+    end
+    @testset "1 < z < 1" begin
+        z = 0.7398
+        x = axes(ChebyshevT(),1)
+        α = 0.1
+        L = abs.(z .- x') .^ α
+        # from Mathematica
+        λ = 1
+        L0 = powerlawmoment(Val(0), α, λ, z)
+        L1 = powerlawmoment(Val(1), α, λ, z)
+        @test L* UltrasphericalWeight(λ) ≈ L0  ≈ 1.480874346601822
+        @test L1 ≈ -0.1328257513137366
+        A, B, C = powerlawrecurrence(α, λ)
+        r = RecurrenceArray(z, (A,B,C), [L0,L1])
+        @test r[5] ≈ 0.02537001925265781
+    end
+end