weights with exponential terms via Clenshaw

Author: TSGut

test/test_decompOPs.jl

@@ -108,7 +108,7 @@ end
     end
 end
 
-@testset "Comparison with Lanczos and Classical, with Clenshaw, basics" begin
+@testset "Comparison with Lanczos and Classical, with Clenshaw, polynomial weights" begin
     @testset "w(x) = x^2*(1-x)" begin
         P = Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])
         x = axes(P,1)
@@ -173,4 +173,66 @@ end
         # Comparison with Lanczos
         @test Jchol[1:500,1:500] ≈ Jlanc[1:500,1:500]
     end
+end
+
+@testset "Comparison with Lanczos and Classical, with Clenshaw, exponential weights" begin
+    @testset "w(x) = exp(x)" begin
+        P = Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])
+        x = axes(P,1)
+        J = jacobimatrix(P)
+        Jx = symmjacobim(J)
+        wf(x) = exp(x)
+        # compute Jacobi matrix via cholesky
+        W = P \ (wf.(x) .* P)
+        Jchol = cholesky_jacobimatrix(Symmetric(W))
+        # compute Jacobi matrix via Lanczos
+        Jlanc = jacobimatrix(LanczosPolynomial(@.(wf.(x)),Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])))
+        # Comparison with Lanczos
+        @test Jchol[1:500,1:500] ≈ Jlanc[1:500,1:500]
+    end
+    
+    @testset "w(x) = (1-x)*exp(x)" begin
+        P = Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])
+        x = axes(P,1)
+        J = jacobimatrix(P)
+        Jx = symmjacobim(J)
+        wf(x) = (1-x)*exp(x)
+        # compute Jacobi matrix via cholesky
+        W = P \ (wf.(x) .* P)
+        Jchol = cholesky_jacobimatrix(Symmetric(W))
+        # compute Jacobi matrix via Lanczos
+        Jlanc = jacobimatrix(LanczosPolynomial(@.(wf.(x)),Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])))
+        # Comparison with Lanczos
+        @test Jchol[1:500,1:500] ≈ Jlanc[1:500,1:500]
+    end
+    
+    @testset "w(x) = (1-x^2)*exp(x^2)" begin
+        P = Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])
+        x = axes(P,1)
+        J = jacobimatrix(P)
+        Jx = symmjacobim(J)
+        wf(x) = (1-x^2)*exp(x^2)
+        # compute Jacobi matrix via cholesky
+        W = P \ (wf.(x) .* P)
+        Jchol = cholesky_jacobimatrix(Symmetric(W))
+        # compute Jacobi matrix via Lanczos
+        Jlanc = jacobimatrix(LanczosPolynomial(@.(wf.(x)),Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])))
+        # Comparison with Lanczos
+        @test Jchol[1:500,1:500] ≈ Jlanc[1:500,1:500]
+    end
+    
+    @testset "w(x) = x*(1-x^2)*exp(-x^2)" begin
+        P = Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])
+        x = axes(P,1)
+        J = jacobimatrix(P)
+        Jx = symmjacobim(J)
+        wf(x) = x*(1-x^2)*exp(-x^2)
+        # compute Jacobi matrix via cholesky
+        W = P \ (wf.(x) .* P)
+        Jchol = cholesky_jacobimatrix(Symmetric(W))
+        # compute Jacobi matrix via Lanczos
+        Jlanc = jacobimatrix(LanczosPolynomial(@.(wf.(x)),Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])))
+        # Comparison with Lanczos
+        @test Jchol[1:500,1:500] ≈ Jlanc[1:500,1:500]
+    end
 end