JuliaApproximation / SingularIntegrals.jl

A Julia package for computing singular integrals
MIT License
6 stars 4 forks source link

logkernel with weighted ChebyshevT #38

Closed kburns closed 6 months ago

kburns commented 8 months ago

I'm trying a simple modification of one of the readme examples, applying a log kernel to a weighted ChebyshevT function. This fails but I can't quite tell what's going on under the hood. Any pointers?

W = Weighted(ChebyshevT())
x = axes(W, 1)
f = expand(W, x -> exp(x) / sqrt(1-x^2))
log.(abs.(0.3 .- x')) * f

returns

ERROR: not implemented
Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:35
  [2] logkernel_layout(::ClassicalOrthogonalPolynomials.WeightedOPLayout{…}, P::Weighted{…}, z::Float64)
    @ SingularIntegrals ~/.julia/packages/SingularIntegrals/LB9Wj/src/logkernel.jl:46
  [3] logkernel(P::Weighted{Float64, ChebyshevT{Float64}}, z::Float64)
    @ SingularIntegrals ~/.julia/packages/SingularIntegrals/LB9Wj/src/logkernel.jl:32
  [4] logkernel_layout(LAY::LazyArrays.ApplyLayout{…}, V::QuasiArrays.ApplyQuasiVector{…}, y::Float64)
    @ SingularIntegrals ~/.julia/packages/SingularIntegrals/LB9Wj/src/logkernel.jl:51
  [5] logkernel_layout(::ContinuumArrays.ExpansionLayout{…}, A::QuasiArrays.ApplyQuasiVector{…}, dims::Float64)
    @ SingularIntegrals ~/.julia/packages/SingularIntegrals/LB9Wj/src/logkernel.jl:54
  [6] logkernel(P::QuasiArrays.ApplyQuasiVector{Float64, typeof(*), Tuple{Weighted{…}, LazyArrays.ApplyArray{…}}}, z::Float64)
    @ SingularIntegrals ~/.julia/packages/SingularIntegrals/LB9Wj/src/logkernel.jl:32
  [7] macro expansion
    @ ~/.julia/packages/SingularIntegrals/LB9Wj/src/logkernel.jl:14 [inlined]
  [8] mul
    @ ~/.julia/packages/ContinuumArrays/amUIf/src/operators.jl:34 [inlined]
  [9] *(A::QuasiArrays.BroadcastQuasiMatrix{…}, B::QuasiArrays.ApplyQuasiVector{…})
    @ QuasiArrays ~/.julia/packages/QuasiArrays/ZCn0F/src/matmul.jl:23
 [10] top-level scope
    @ REPL[108]:1
Some type information was truncated. Use `show(err)` to see complete types.
dlfivefifty commented 8 months ago

Just means it hasn't been implemented yet. I guess I only got around to weighted Chebyshev U... There is a general approach for classical OPs but need to add it.

If you only need to evaluate the log kernel on the interval you can do this:

julia> using ClassicalOrthogonalPolynomials, SingularIntegrals

julia> W = Weighted(ChebyshevT())
Weighted(ChebyshevT())

julia> x = axes(W, 1)
Inclusion(-1.0 .. 1.0 (Chebyshev))

julia> f = expand(W, x -> exp(x) / sqrt(1-x^2))
Weighted(ChebyshevT()) * [1.2660658777520084, 1.1303182079849703, 0.27149533953407656, 0.044336849848663776, 0.0054742404420936785, 0.0005429263119139354, 4.497732295427654e-5, 3.198436462529006e-6, 1.992124804817033e-7, 1.1036771880698548e-8  …  ]

julia> x = axes(W, 1)
Inclusion(-1.0 .. 1.0 (Chebyshev))

julia> @time Lf = log.(abs.(x .- x')) * f;
  0.000020 seconds (7 allocations: 496 bytes)

julia> @time Lf[0.3]
  0.000008 seconds (1 allocation: 16 bytes)
-3.437622681488403
kburns commented 8 months ago

Thanks. How is (log.(abs.(x .- x')) * f)[0.3] different under the hood from (log.(abs.(0.3 .- x')) * f)?

dlfivefifty commented 8 months ago

The former deduces the Chebyshev coefficients by applying a diagonal matrix:

julia> T = ChebyshevT(); W = Weighted(T);

julia> T \ logkernel(W)
ℵ₀×ℵ₀ LinearAlgebra.Diagonal{Float64, LazyArrays.ApplyArray{Float64, 1, typeof(vcat), Tuple{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{Float64, InfiniteArrays.InfUnitRange{Int64}}}}}} with indices OneToInf()×OneToInf():
 -2.17759    ⋅         ⋅        ⋅      …  
   ⋅       -3.14159    ⋅        ⋅         
   ⋅         ⋅       -1.5708    ⋅         
   ⋅         ⋅         ⋅      -1.0472     
   ⋅         ⋅         ⋅        ⋅         
   ⋅         ⋅         ⋅        ⋅      …  
   ⋅         ⋅         ⋅        ⋅         
   ⋅         ⋅         ⋅        ⋅         
   ⋅         ⋅         ⋅        ⋅         
   ⋅         ⋅         ⋅        ⋅         
  ⋮                                    ⋱  

(I don't know if I documented it but log.(abs(x .- x')) * f is a synonym for logkernel(f). )