Fixed bjgs

Author: dlfivefifty

examples/ultrasphericalspectralmethod.jl

@@ -0,0 +1,14 @@
+using ContinuumArrays, FillArrays, InfiniteArrays, Plots
+
+T = Chebyshev()
+C = Ultraspherical(2)
+D = Derivative(axes(T,1))
+A = (C\(D^2*T))+100(C\T)
+
+n = 100
+c = [T[[-1,1],1:n]; A[1:n-2,1:n]] \ [1;2;zeros(n-2)]
+u = T*Vcat(c,Zeros(∞))
+
+xx = range(-1,1;length=1000)
+plot(xx,u[xx])
+

src/ContinuumArrays.jl

@@ -14,7 +14,9 @@ import QuasiArrays: cardinality, checkindex, QuasiAdjoint, QuasiTranspose, Inclu
                     ApplyQuasiArray, ApplyQuasiMatrix, LazyQuasiArrayApplyStyle, AbstractQuasiArrayApplyStyle,
                     LazyQuasiArray, LazyQuasiVector, LazyQuasiMatrix, LazyLayout, LazyQuasiArrayStyle
 
-export Spline, LinearSpline, HeavisideSpline, DiracDelta, Derivative, JacobiWeight, Jacobi, Legendre, Chebyshev, Ultraspherical,
+export Spline, LinearSpline, HeavisideSpline, DiracDelta, Derivative, 
+            Jacobi, Legendre, Chebyshev, Ultraspherical,
+            JacobiWeight, ChebyshevWeight, UltrasphericalWeight,
             fullmaterialize
 
 ####

src/bases/jacobi.jl

@@ -51,7 +51,7 @@ function jacobimatrix(J::Jacobi)
     b,a = J.b,J.a
     n = 0:∞
     B = @. 2*(n+1)*(n+a+b+1) / ((2n+a+b+1)*(2n+a+b+2))
-    A = (a^2-b^2) ./ ((2n.+a.+b).*(2n.+a.+b.+2))
+    A = Vcat((a-b) / (a+b+2), (a^2-b^2) ./ ((2n.+a.+b.+2).*(2n.+a.+b.+4)))
     C = @. 2*(n+a)*(n+b) / ((2n+a+b)*(2n+a+b+1))
 
     _BandedMatrix(Vcat(C',A',B'), ∞, 1,1)

src/bases/orthogonalpolynomials.jl

@@ -13,7 +13,7 @@ end
 
 function forwardrecurrence!(v::AbstractVector{T}, b::AbstractVector, a::AbstractVector, c::AbstractVector, x) where T
     isempty(v) && return v
-    v[1] = one(T) # assume OPs are normalized to one for now
+    v[1] = one(x) # assume OPs are normalized to one for now
     length(v) == 1 && return v
     v[2] = (x-a[1])/c[1]
     @inbounds for n=3:length(v)
@@ -24,7 +24,7 @@ end
 
 function forwardrecurrence!(v::AbstractVector{T}, b::AbstractVector, ::Zeros{<:Any,1}, c::AbstractVector, x) where T
     isempty(v) && return v
-    v[1] = one(T) # assume OPs are normalized to one for now
+    v[1] = one(x) # assume OPs are normalized to one for now
     length(v) == 1 && return v
     v[2] = x/c[1]
     @inbounds for n=3:length(v)
@@ -37,7 +37,7 @@ end
 function forwardrecurrence!(v::AbstractVector{T}, b::AbstractVector, ::Zeros{<:Any,1}, c::Vcat{<:Any,1,<:Tuple{<:Number,<:AbstractVector}}, x) where T
     isempty(v) && return v
     c0,c∞ = c.args
-    v[1] = one(T) # assume OPs are normalized to one for now
+    v[1] = one(x) # assume OPs are normalized to one for now
     length(v) == 1 && return v
     v[2] = x/c0
     @inbounds for n=3:length(v)
@@ -53,7 +53,7 @@ function forwardrecurrence!(v::AbstractVector{T}, b_v::AbstractFill, ::Zeros{<:A
     c∞ = getindex_value(c∞_v) 
     mbc  = -b/c∞
     xc = x/c∞
-    v[1] = one(T) # assume OPs are normalized to one for now
+    v[1] = one(x) # assume OPs are normalized to one for now
     length(v) == 1 && return v
     v[2] = x/c0
     @inbounds for n=3:length(v)

src/bases/splines.jl

@@ -31,9 +31,10 @@ function getindex(B::HeavisideSpline{T}, x::Number, k::Int) where T
     p = B.points
     n = length(p)
 
-    x < p[k] && return zero(T)
-    k < n && x > p[k+1] && return zero(T)
-    return one(T)
+    p[k] < x < p[k+1] && return one(T)
+    p[k] == x && return one(T)/2
+    p[k+1] == x && return one(T)/2
+    return zero(T)
 end
 
 ## Sub-bases

src/bases/ultraspherical.jl

@@ -9,6 +9,8 @@ struct UltrasphericalWeight{T,Λ} <: AbstractJacobiWeight{T}
     λ::Λ
 end
 
+UltrasphericalWeight(λ) = UltrasphericalWeight{typeof(λ),typeof(λ)}(λ)
+
 function getindex(w::UltrasphericalWeight, x::Number)
     x ∈ axes(w,1) || throw(BoundsError())
     (1-x^2)^(w.λ-one(w.λ)/2)
@@ -34,10 +36,8 @@ Ultraspherical(λ::Λ) where Λ = Ultraspherical{Float64,Λ}(λ)
 
 jacobimatrix(C::Chebyshev{T}) where T = _BandedMatrix(Vcat(Fill(one(T)/2,1,∞), Zeros(1,∞), Hcat(one(T), Fill(one(T)/2,1,∞))), ∞, 1, 1)
 
-function jacobimatrix(P::Ultraspherical)
-    T = promote_type(eltype(P), eltype(X), eltype(U))
+function jacobimatrix(P::Ultraspherical{T}) where T
     λ = P.λ
-    λ == U.λ || throw(ArgumentError())
     _BandedMatrix(Vcat((((2λ-1):∞) ./ (2 .*((zero(T):∞) .+ λ)))',
                         Zeros{T}(1,∞),
                         ((one(T):∞) ./ (2 .*((zero(T):∞) .+ λ)))'), ∞, 1, 1)
@@ -50,6 +50,7 @@ end
 
 # Ultraspherical(1)\(D*Chebyshev())
 @simplify function \(J::Ultraspherical, *(D::Derivative{<:Any,<:ChebyshevInterval}, S::Chebyshev))
+    T = promote_type(eltype(J),eltype(D),eltype(S))
     (J.λ == 1) || throw(ArgumentError())
     _BandedMatrix((zero(eltype(M)):∞)', ∞, -1,1)
 end
@@ -59,6 +60,19 @@ end
     ApplyQuasiMatrix(*, Ultraspherical{eltype(S)}(1), A)
 end
 
+# Ultraspherical(1/2)\(D*Legendre())
+@simplify function \(J::Ultraspherical, *(D::Derivative{<:Any,<:ChebyshevInterval}, S::Legendre))
+    T = promote_type(eltype(J),eltype(D),eltype(S))
+    (J.λ == 3/2) || throw(ArgumentError())
+    _BandedMatrix(Ones{T}(1,∞), ∞, -1,1)
+end
+
+@simplify function *(D::Derivative{<:Any,<:ChebyshevInterval}, S::Legendre)
+    A = apply(\,Ultraspherical{eltype(S)}(3/2),applied(*,D,S))
+    ApplyQuasiMatrix(*, Ultraspherical{eltype(S)}(3/2), A)
+end
+
+
 # Ultraspherical(λ+1)\(D*Ultraspherical(λ))
 @simplify function \(J::Ultraspherical, *(D::Derivative{<:Any,<:ChebyshevInterval}, S::Ultraspherical))
     (J.λ == S.λ+1) || throw(ArgumentError())
@@ -76,21 +90,28 @@ end
 ##########
 
 
-@simplify function \(U::Ultraspherical, C::Chebyshev)
-    (U.λ == 1) || throw(ArgumentError())
-    T = promote_type(eltype(U), eltype(C))
-    BandedMatrix(0 => Vcat([one(T)],Fill(one(T)/2,∞)), 2 => Vcat([-one(T)/2],Fill(-one(T)/2,∞)))
+@simplify function \(U::Ultraspherical{<:Any,<:Integer}, C::Chebyshev)
+    if U.λ == 1
+        T = promote_type(eltype(U), eltype(C))
+        BandedMatrix(0 => Vcat([one(T)],Fill(one(T)/2,∞)), 2 => Vcat([-one(T)/2],Fill(-one(T)/2,∞)))
+    elseif U.λ > 0
+        (U\Ultraspherical(1)) * (Ultraspherical(1)\C)
+    else
+        error("Not implemented")
+    end
 end
 
-@simplify function \(C2::Ultraspherical, C1::Ultraspherical)
+@simplify function \(C2::Ultraspherical{<:Any,<:Integer}, C1::Ultraspherical{<:Any,<:Integer})
     λ = C1.λ
     T = promote_type(eltype(C2), eltype(C1))
     if C2.λ == λ+1 
         _BandedMatrix( Vcat(-(λ ./ (1:∞ .+ λ))', Zeros(1,∞), (λ ./ (1:∞ .+ λ))'), ∞, 0, 2)
     elseif C2.λ == λ
         Eye{T}(∞)
+    elseif C2.λ > λ
+        (C2 \ Ultraspherical(λ+1)) * (Ultraspherical(λ+1)\C1)
     else
-        throw(ArgumentError())
+        error("Not implemented")
     end
 end
 

src/operators.jl

@@ -151,6 +151,7 @@ Derivative(axis) = Derivative{Float64}(axis)
 
 axes(D::Derivative) = (D.axis, D.axis)
 ==(a::Derivative, b::Derivative) = a.axis == b.axis
+copy(D::Derivative) = Derivative(copy(D.axis))
 
 
 function copy(M::QMul2{<:Derivative,<:SubQuasiArray})

test/test_orthogonalpolynomials.jl

@@ -1,4 +1,7 @@
-import ContinuumArrays: jacobimatrix
+using ContinuumArrays, FillArrays, LazyArrays, BandedMatrices
+import ContinuumArrays: jacobimatrix, SimplifyStyle, ∞
+import LazyArrays: ApplyStyle, colsupport
+
 
 @testset "Ultraspherical" begin
     T = Chebyshev()
@@ -19,7 +22,7 @@ import ContinuumArrays: jacobimatrix
     @test colsupport(D₀,1) == 1:0
 
     D₁ = C\(D*U)
-    @test D₁ isa BandedMatrix
+    @test_broken D₁ isa BandedMatrix
     @test apply(*,D₁,D₀.args...)[1:10,1:10] == diagm(2 => 4:2:18)
     @test (D₁*D₀)[1:10,1:10] == diagm(2 => 4:2:18)
 
@@ -32,6 +35,36 @@ import ContinuumArrays: jacobimatrix
     @test S₁ == diagm(0 => 1 ./ (1:10), 2=> -(1 ./ (3:10)))
 end
 
+@testset "Legendre" begin
+    @test jacobimatrix(Jacobi(0.,0.))[1,1] == 0.0
+    @test jacobimatrix(Jacobi(0.,0.))[1:10,1:10] == jacobimatrix(Legendre())[1:10,1:10] == jacobimatrix(Ultraspherical(1/2))[1:10,1:10]
+    @test Jacobi(0.,0.)[0.1,1:10] ≈ Legendre()[0.1,1:10] ≈ Ultraspherical(1/2)[0.1,1:10]
+
+    P = Legendre()
+    D = Derivative(axes(P,1))
+    @test Ultraspherical(3/2)\(D*P) isa BandedMatrix{Float64,<:Ones}
+end
+
+
+
+@testset "Jacobi" begin
+    b,a = 0.1,0.2
+    S = Jacobi(b,a)
+    x = 0.1
+    @test S[x,1] === 1.0
+    X = jacobimatrix(S)
+    @test X[1,1] ≈ (a^2-b^2)/((a+b)*(a+b+2))
+    @test X[2,1] ≈ 2/(a+b+2)
+    @test S[x,2] ≈ 0.065
+    @test S[x,10] ≈ 0.22071099583604945
+
+    w = JacobiWeight(b,a)
+    @test w[x] ≈ (1+x)^b * (1-x)^a
+    wS = w.*S
+    @test wS[0.1,1] ≈ w[0.1]
+    @test wS[0.1,1:2] ≈ w[0.1] .* S[0.1,1:2]
+end
+
 @testset "Jacobi integer" begin
     S = Jacobi(true,true)
     W = Diagonal(JacobiWeight(true,true))
@@ -103,18 +136,6 @@ end
     @test W[0.1,0.2] ≈ 0.0
 end
 
-@testset "Jacobi" begin
-    b,a = 0.1,0.2
-    S = Jacobi(b,a)
-    x = 0.1
-    @test S[x,1] === 1.0
-    X = jacobimatrix(S)
-    @test X[1,1] ≈ (a^2-b^2)/((a+b)*(a+b+2))
-    @test X[2,1] ≈ 2/(a+b+2)
-    @test S[x,2] ≈ 0.065
-    @test S[x,10] ≈ 0.22071099583604945
-end
-
 @testset "P-FEM" begin
     S = Jacobi(true,true)
     W = Diagonal(JacobiWeight(true,true))