tpapp / spectralkit.jl Goto Github PK

Building blocks of spectral methods for Julia.

License: Other

Julia 100.00%

chebyshev spectral-methods julia julia-language package approximation

spectralkit.jl's Introduction

SpectralKit.jl

Building blocks of spectral methods for Julia. Currently includes Chebyshev polynomials on univariate and Smolyak (multivariate) grids, with domain transformations to semi-infinite and infinite domains.

Introduction

Mostly useful for algorithms along the lines of

Boyd, John P. Chebyshev and Fourier spectral methods. 2001.

The aim is to provide simple, well-tested, robust, and fast building blocks for spectral algorithms, which can be easily combined into algorithms.

At the moment, the package API is experimental and subject to change.

Help

Asking for help in issues is fine, you can also ping me as @Tamas_Papp on the Discourse forum

Pretty pictures

Some examples generated this library. Circles mark values at the limit, shifted horizontally when this is needed to avoid overlap. Infinite limits shown at finite values, so of course they don't match (this is a visual check of continuity, naturally it is unit tested).

Chebyshev polynomials and their derivatives

Chebyshev rational functions on [0,∞)

Up close, you can see the oscillation.

Let's zoom out a bit to see convergence to 0 at ∞.

Derivatives die out faster.

Chebyshev rational functions on (-∞,∞)

Up close, you can see the oscillation.

Let's zoom out a bit to see convergence at -∞ and ∞.

Derivatives die out slower than for the [0,∞) transformation.

A Smolyak grid

With B = 3.

Bibliography

Boyd, J. P. (2001). Chebyshev and fourier spectral methods. Courier Corporation.
Xu, K. (2016). The chebyshev points of the first kind. Applied Numerical Mathematics, 102, 17–30.

spectralkit.jl's People

Contributors

Stargazers

Watchers

spectralkit.jl's Issues

linear_combination should chech length

Currently it doesn't, which can lead to insidious bugs.

julia> using SpectralKit

julia> basis = Chebyshev(InteriorGrid(), 5)
Chebyshev polynomials (1st kind), InteriorGrid(), dimension: 5

julia> linear_combination(basis, zeros(dimension(basis) + 1), 0.0) # NOTE + 1
0.0

TagBot trigger issue

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

make dimension part of the type

This would allow allocation-free code (using SVector) for a lot of calculations. Of course optionally, as users can still go for Vectors etc.

Caveat: output of collocation_matrix may be too large. Opt-in to static matrices?

domain checking

Evaluation of basis functions and linear combinations should support various forms of domain checking, eg

no checking (as currently done), this will just error because of acos for Chebyshev
strict checking (terminate with an informative DomainError
"soft" extrapolation: when within a small distance of boundaries, just use the boundary value (mostly for dealing with floating point error)

Cf #7, a refactor should implement both.

construct notation for partial derivatives independently of a value (∂(Val(N), (), (1,)))
wrap a value for evaluation with partial derivatives (∂(∂spec, x))
a shorthand that combines the two, (∂(x, partials...))

Methods are

julia> methods(∂)
# 5 methods for generic function "∂" from SpectralKit:
 [1] ∂(x::Tuple, partials...)
     @ ~/.julia/packages/SpectralKit/gEVDA/src/derivatives.jl:240
 [2] ∂(x::StaticArraysCore.SVector{N}, partials...) where N
     @ ~/.julia/packages/SpectralKit/gEVDA/src/derivatives.jl:235
 [3] ∂(∂specification::SpectralKit.∂Specification{M}, x::Union{Tuple, AbstractVector}) where M
     @ ~/.julia/packages/SpectralKit/gEVDA/src/derivatives.jl:224
 [4] ∂(N::Integer, partials...)
     @ ~/.julia/packages/SpectralKit/gEVDA/src/derivatives.jl:179
 [5] ∂(::Val{N}, partials...) where N
     @ ~/.julia/packages/SpectralKit/gEVDA/src/derivatives.jl:171

This is too complex and inconsistent. Yet it kind of made sense in the first iteration.

Proposal:

have a single ∂(spec...) operator that is dimension-independent, to the extend that it encodes the minimum dimension it needs to work with. Otherwise the entries in the lookup are implicitly extented with zeros. This dispenses with the Val as the first argument.
make this callable, wrapping a vector-like object (SVector, Tuple, etc).

Eg a user could do

∂spec = ∂((), (1,), (2,))
linear_combination(basis, coefficients, ∂spec(x))

and have it work for all x such that length(x) >= 2.

implement / for derivatives

Was missed in #35, necessary for infinite and semi-infinite mappings.

check that grid spec is valid

Currently something like

julia> Chebyshev(AbstractMatrix, 10)
Chebyshev polynomials (1st kind), AbstractMatrix, dimension: 10

will fail later on instead of immediately.

Didn't check for Smolyak, do that too.