Add normalized hermite eigevalue equation

Project.toml

@@ -1,7 +1,7 @@
 name = "ClassicalOrthogonalPolynomials"
 uuid = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd"
 authors = ["Sheehan Olver <[email protected]>"]
-version = "0.2.2"
+version = "0.2.3"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -37,7 +37,7 @@ FillArrays = "0.11"
 InfiniteArrays = "0.10"
 InfiniteLinearAlgebra = "0.5.2"
 IntervalSets = "0.5"
-LazyArrays = "0.20.6"
+LazyArrays = "0.20.8"
 LazyBandedMatrices = "0.5"
 QuasiArrays = "0.4.6"
 SpecialFunctions = "0.10, 1"

src/hermite.jl

@@ -7,9 +7,11 @@ function getindex(w::HermiteWeight, x::Number)
     exp(-x^2)
 end
 
+sum(::HermiteWeight{T}) where T = sqrt(convert(T, π))
 
 struct Hermite{T} <: OrthogonalPolynomial{T} end
 Hermite() = Hermite{Float64}()
+orthogonalityweight(::Hermite{T}) where T = HermiteWeight{T}()
 
 ==(::Hermite, ::Hermite) = true
 axes(::Hermite{T}) where T = (Inclusion(ℝ), oneto(∞))
@@ -33,4 +35,9 @@ end
     T = promote_type(eltype(D),eltype(H))
     D = _BandedMatrix((zero(T):2:∞)', ℵ₀, -1,1)
     H*D
+end
+
+@simplify function *(D::Derivative, Q::OrthonormalWeighted{<:Any,<:Hermite})
+    X = jacobimatrix(Q.P)
+    Q * Tridiagonal(-X.ev, X.dv, X.ev)
 end

src/normalized.jl

@@ -47,7 +47,7 @@ end
 normalizationconstant(P) = NormalizationConstant(P)
 Normalized(P::AbstractQuasiMatrix{T}) where T = Normalized(P, normalizationconstant(P))
 Normalized(Q::Normalized) = Q
-
+normalized(P) = Normalized(P)
 
 struct NormalizedBasisLayout{LAY<:AbstractBasisLayout} <: AbstractBasisLayout end
 
@@ -166,3 +166,25 @@ end
 Base.array_summary(io::IO, C::NormalizationConstant{T}, inds) where T = print(io, "NormalizationConstant{$T}")
 show(io::IO, Q::Normalized) = print(io, "Normalized($(Q.P))")
 show(io::IO, ::MIME"text/plain", Q::Normalized) = show(io, Q)
+
+
+
+struct OrthonormalWeighted{T, PP<:AbstractQuasiMatrix{T}} <: Basis{T}
+    P::Normalized{T, PP, NormalizationConstant{T, PP}}
+end
+
+function OrthonormalWeighted(P)
+    Q = normalized(P)
+    OrthonormalWeighted{eltype(Q),typeof(P)}(Q)
+end
+
+axes(Q::OrthonormalWeighted) = axes(Q.P)
+copy(Q::OrthonormalWeighted) = Q
+
+==(A::OrthonormalWeighted, B::OrthonormalWeighted) = A.P == B.P
+
+function getindex(Q::OrthonormalWeighted, x::Union{Number,AbstractVector}, jr::Union{Number,AbstractVector})
+    w = orthogonalityweight(Q.P)
+    sqrt.(w[x]) .* Q.P[x,jr]
+end
+broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), x::Inclusion, Q::OrthonormalWeighted) = Q * (Q.P \ (x .* Q.P))

test/test_hermite.jl

@@ -1,23 +1,47 @@
 using ClassicalOrthogonalPolynomials, ContinuumArrays, FillArrays, Test
-import ClassicalOrthogonalPolynomials: jacobimatrix, oneto
+import ClassicalOrthogonalPolynomials: jacobimatrix, oneto, OrthonormalWeighted
 import DomainSets: ℝ
 
-@testset "Hermite" begin
-    H = Hermite()
-    w = HermiteWeight()
-
-    @test axes(H) == (Inclusion(ℝ), oneto(∞))
-    x = axes(H,1)
-    @test H[0.1,1:4] ≈ [1,2*0.1,4*0.1^2-2,8*0.1^3-12*0.1]
-    @test w[0.1] ≈ exp(-0.1^2)
+@testset "Hermite" begin    
+    @testset "Basics" begin
+        H = Hermite()
+        w = HermiteWeight()
+        @test axes(H) == (Inclusion(ℝ), oneto(∞))
+        x = axes(H,1)
+        @test H[0.1,1:4] ≈ [1,2*0.1,4*0.1^2-2,8*0.1^3-12*0.1]
+        @test w[0.1] ≈ exp(-0.1^2)
 
-    X = jacobimatrix(H)
-    @test 0.1 * H[0.1,1:10]' ≈ H[0.1,1:11]'*X[1:11,1:10]
+        X = jacobimatrix(H)
+        @test 0.1 * H[0.1,1:10]' ≈ H[0.1,1:11]'*X[1:11,1:10]
 
-    @test (H'*(w .* H))[1,1] ≈ sqrt(π)
-    @test (H'*(w .* H))[2,2] ≈ 2sqrt(π)
-    @test (H'*(w .* H))[3,3] ≈ 8sqrt(π)
+        @test (H'*(w .* H))[1,1] ≈ sqrt(π)
+        @test (H'*(w .* H))[2,2] ≈ 2sqrt(π)
+        @test (H'*(w .* H))[3,3] ≈ 8sqrt(π)
 
-    D = Derivative(x)
-    @test (D*H)[0.1,1:4] ≈ [0,2,8*0.1,24*0.1^2-12]
+        D = Derivative(x)
+        @test (D*H)[0.1,1:4] ≈ [0,2,8*0.1,24*0.1^2-12]
+    end
+
+    @testset "OrthonormalWeighted" begin
+        H = Hermite()
+        Q = OrthonormalWeighted(H)
+        @testset "evaluation" begin
+            x = 0.1
+            @test Q[x,1] ≈ exp(-x^2/2)/π^(1/4)
+            @test Q[x,2] ≈ 2x*exp(-x^2/2)/(sqrt(2)π^(1/4))
+            # Trap test of L^2 orthonormality
+            x = range(-20,20; length=1_000)
+            @test Q[x,1:5]'Q[x,1:5] * step(x) ≈ I
+        end
+
+        @testset "Differentiation" begin
+            x = axes(Q,1)
+            D = Derivative(x)
+            D¹ = Q \ (D * Q)
+            @test D¹[1:10,1:10] ≈ -D¹[1:10,1:10]'
+            D² = Q \ (D^2 * Q)
+            X = Q \ (x .* Q)
+            @test (D² - X^2)[1:10,1:10] ≈ -Diagonal(1:2:19)
+        end
+    end
 end