Support Clenshaw with Weighted (#262)

* Support Clenshaw with Weighted

* non symmetric example

Author: dlfivefifty

Project.toml

@@ -1,6 +1,6 @@
 name = "ClassicalOrthogonalPolynomials"
 uuid = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd"
-version = "0.15.12"
+version = "0.15.13"
 authors = ["Sheehan Olver <[email protected]>"]
 
 [deps]

examples/nonlinearode.jl

@@ -0,0 +1,27 @@
+####
+# solve
+#
+# -u'' + u^2 = 1 + x
+# u(±1) = 0
+# 
+# using integrated Legendre w/ Newton iteration.
+###
+
+using ClassicalOrthogonalPolynomials, Plots
+
+P = Legendre(); C = Ultraspherical(3/2); W = Weighted(C)
+x = axes(P,1)
+Δ = diff(W)'diff(W)
+M = W'W
+b = W'*(1 .+ x)
+F = u -> Δ*u + M * (W\((W*u) .^ 2)) - b
+J = u -> Δ + 2*M * (W \ ((P * (P\(W*u))) .* W))
+
+
+u = zeros(∞)
+while norm(F(u)) ≥ 1E-12
+    u -= J(u) \ F(u)
+end
+
+
+plot(W*u)

src/clenshaw.jl

@@ -107,12 +107,18 @@ Base.@propagate_inbounds getindex(f::Mul{<:WeightedOPLayout,<:AbstractPaddedLayo
     weight(f.A)[x] * (unweighted(f.A) * f.B)[x, j...]
 
 
-# TODO: generalise this to be trait based
 function layout_broadcasted(::Tuple{ExpansionLayout{<:AbstractOPLayout},AbstractOPLayout}, ::typeof(*), a, P)
     axes(a,1) == axes(P,1) || throw(DimensionMismatch())
     P * Clenshaw(a, P)
 end
 
+function layout_broadcasted(::Tuple{ExpansionLayout{<:AbstractOPLayout},WeightedOPLayout}, ::typeof(*), a, W)
+    P = W.P
+    axes(a,1) == axes(P,1) || throw(DimensionMismatch())
+    W * Clenshaw(a, P)
+end
+
+
 # TODO: layout_broadcasted
 function broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), a::ApplyQuasiVector{<:Any,typeof(*),<:Tuple{SubQuasiArray{<:Any,2,<:OrthogonalPolynomial,<:Tuple{AbstractAffineQuasiVector,Slice}},Any}}, V::SubQuasiArray{<:Any,2,<:OrthogonalPolynomial,<:Tuple{AbstractAffineQuasiVector,Any}})
     axes(a,1) == axes(V,1) || throw(DimensionMismatch())

test/test_chebyshev.jl

@@ -595,6 +595,13 @@ import BandedMatrices: isbanded
         op = 2 * f .* T
         @test op[0.1,1:5] ≈ 2f[0.1] * T[0.1,1:5]
     end
+
+    @testset "Clenshaw for Weighted" begin
+        T = ChebyshevT(); W = Weighted(T)
+        a = expand(T, exp)
+        @test W \ (a .* W) isa Clenshaw
+        @test (a .* W)[0.1,1:5] ≈ a[0.1] * W[0.1,1:5]
+    end
 end
 
 struct QuadraticMap{T} <: Map{T} end