Add weighted Hermite Derivative (#125)

* Add weighted Hermite Derivative

* add Quartic example

Author: dlfivefifty

Project.toml

@@ -1,7 +1,7 @@
 name = "ClassicalOrthogonalPolynomials"
 uuid = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd"
 authors = ["Sheehan Olver <[email protected]>"]
-version = "0.7.3"
+version = "0.7.4"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"

examples/quarticev.jl

@@ -0,0 +1,50 @@
+using ClassicalOrthogonalPolynomials
+import ClassicalOrthogonalPolynomials: OrthonormalWeighted
+
+# find the 1m eigenvalue of -u'' + x^4*u
+
+H = Hermite()
+Q = OrthonormalWeighted(H)
+x = axes(H,1)
+X = Q'*(x.*Q)
+D = Derivative(x)
+Δ = Q'*(D^2*Q)
+
+L = -Δ + X^4
+
+n = 1_000_000
+Lₙ = L[1:n,1:n]
+Fₙ = factorize(Symmetric(Lₙ));
+
+v = [1; zeros(n-1)];
+
+# find the smallest eigenvalue with inverse iteration
+for _=1:10
+    w = Fₙ \ v; v = w/norm(w); v
+end
+λ = v'*Lₙ*v
+
+@test λ ≈ eigvals(Symmetric(L[1:100,1:100]))[1]
+
+# we can acellerate using Rayleigh shifts
+
+v = [1; zeros(n-1)];
+
+# find the millionth eigenvalue with inverse iteration
+n = 50_000_000
+E = 218506784;
+Lₙ = (L-E*I)[1:n,1:n]
+
+Fₙ = factorize(Symmetric(Lₙ));
+
+v = randn(n)
+
+
+for _=1:10
+    w = Fₙ \ v; v = w/norm(w); v
+    @show λ = v'*Lₙ*v
+end
+
+
+
+

src/classical/hermite.jl

@@ -51,6 +51,12 @@ massmatrix(::Hermite{T}) where T = Diagonal(sqrt(convert(T,π)) .* convert(T,2)
     H*D
 end
 
+@simplify function *(D::Derivative, W::Weighted{<:Any,<:Hermite})
+    T = promote_type(eltype(D),eltype(W))
+    D = _BandedMatrix(Fill(-one(T), 1, ∞), ℵ₀, 1,-1)
+    W*D
+end
+
 @simplify function *(D::Derivative, Q::OrthonormalWeighted{<:Any,<:Hermite})
     X = jacobimatrix(Q.P)
     Q * Tridiagonal(-X.ev, X.dv, X.ev)

src/normalized.jl

@@ -218,6 +218,11 @@ weight(Q::OrthonormalWeighted) = sqrt.(orthogonalityweight(Q.P))
 
 broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), x::Inclusion, Q::OrthonormalWeighted) = Q * (Q.P \ (x .* Q.P))
 
+@simplify function *(Ac::QuasiAdjoint{<:Any,<:OrthonormalWeighted}, B::OrthonormalWeighted)
+    @assert parent(Ac).P == B.P
+    Eye{promote_type(eltype(Ac),eltype(B))}(∞)
+end
+
 
 
 """

test/test_hermite.jl

@@ -22,9 +22,20 @@ import DomainSets: ℝ
         @test (D*H)[0.1,1:4] ≈ [0,2,8*0.1,24*0.1^2-12]
     end
 
+    @testset "Weighted" begin
+        H = Hermite()
+        W = Weighted(H)
+        x = axes(W,1)
+        D = Derivative(x)
+        @test (D*W)[0.1,1:4] ≈ [-2*0.1, 2-4*0.1^2, 12*0.1 - 8*0.1^3, -4*(3 - 12*0.1^2 + 4*0.1^4)]*exp(-0.1^2)
+    end
+
     @testset "OrthonormalWeighted" begin
         H = Hermite()
         Q = OrthonormalWeighted(H)
+
+        @test Q'Q == Eye(∞)
+
         @testset "evaluation" begin
             x = 0.1
             @test Q[x,1] ≈ exp(-x^2/2)/π^(1/4)