From bc88a147d73548414d7e22cca3c6133b23b96a6d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 11 Jul 2024 23:24:53 +0200 Subject: [PATCH] Trigonometric basis --- src/MultivariateBases.jl | 3 +- src/chebyshev.jl | 29 ++++++------ src/trigonometric.jl | 96 ++++++++++++++++++++++++++++++++++++++++ 3 files changed, 114 insertions(+), 14 deletions(-) create mode 100644 src/trigonometric.jl diff --git a/src/MultivariateBases.jl b/src/MultivariateBases.jl index 8e35e58..bd23b67 100644 --- a/src/MultivariateBases.jl +++ b/src/MultivariateBases.jl @@ -47,7 +47,7 @@ include("scaled.jl") export AbstractMultipleOrthogonal, ProbabilistsHermite, PhysicistsHermite, Laguerre export AbstractGegenbauer, - Legendre, Chebyshev, ChebyshevFirstKind, ChebyshevSecondKind + Legendre, Chebyshev, ChebyshevFirstKind, ChebyshevSecondKind, Trigonometric export algebra_element, sparse_coefficients, univariate_orthogonal_basis, @@ -62,6 +62,7 @@ include("hermite.jl") include("laguerre.jl") include("legendre.jl") include("chebyshev.jl") +include("trigonometric.jl") include("lagrange.jl") include("quotient.jl") diff --git a/src/chebyshev.jl b/src/chebyshev.jl index 35619c1..46b57cc 100644 --- a/src/chebyshev.jl +++ b/src/chebyshev.jl @@ -24,11 +24,23 @@ const Chebyshev = ChebyshevFirstKind # https://en.wikipedia.org/wiki/Chebyshev_polynomials#Properties # `T_n * T_m = T_{n + m} / 2 + T_{|n - m|} / 2` -function (::Mul{Chebyshev})(a::MP.AbstractMonomial, b::MP.AbstractMonomial) +function univariate_mul!(::Mul{Chebyshev}, terms, var, a, b) + I = eachindex(terms) + for i in I + mono = MP.monomial(terms[i]) * var^(a + b) + terms[i] = MA.mul!!(terms[i], var^abs(a - b)) + terms[i] = MA.operate!!(/, terms[i], 2) + α = MA.copy_if_mutable(MP.coefficient(terms[i])) + push!(terms, MP.term(α, mono)) + end + return +end + +function (mul::Mul{B})(a::MP.AbstractMonomial, b::MP.AbstractMonomial) where {B<:AbstractMonomialIndexed} terms = [MP.term(1 // 1, MP.constant_monomial(a * b))] - vars_a = MP.variables(a) + vars_a = MP.effective_variables(a) var_state_a = iterate(vars_a) - vars_b = MP.variables(b) + vars_b = MP.effective_variables(b) var_state_b = iterate(vars_b) while !isnothing(var_state_a) || !isnothing(var_state_b) if isnothing(var_state_a) || @@ -48,16 +60,7 @@ function (::Mul{Chebyshev})(a::MP.AbstractMonomial, b::MP.AbstractMonomial) else var_a, state_a = var_state_a var_b, state_b = var_state_b - d_a = MP.degree(a, var_a) - d_b = MP.degree(b, var_b) - I = eachindex(terms) - for i in I - mono = MP.monomial(terms[i]) * var_a^(d_a + d_b) - terms[i] = MA.mul!!(terms[i], var_a^abs(d_a - d_b)) - terms[i] = MA.operate!!(/, terms[i], 2) - α = MA.copy_if_mutable(MP.coefficient(terms[i])) - push!(terms, MP.term(α, mono)) - end + univariate_mul!(mul, terms, var_a, MP.degree(a, var_a), MP.degree(b, var_b)) var_state_a = iterate(vars_a, state_a) var_state_b = iterate(vars_b, state_b) end diff --git a/src/trigonometric.jl b/src/trigonometric.jl new file mode 100644 index 0000000..2f5305a --- /dev/null +++ b/src/trigonometric.jl @@ -0,0 +1,96 @@ +""" + struct Trigonometric <: AbstractMonomialIndexed end + +Univariate trigonometric basis is +``` +a0 + a1 cos(ωt) + a2 sin(ωt) + a3 cos(2ωt) + a4 sin(2ωt) +``` +""" +struct Trigonometric <: AbstractMonomialIndexed end + +_cos_id(d) = iszero(d) ? 0 : 2d - 1 +_sin_id(d) = 2d +_is_cos(d) = isodd(d) +_is_sin(d) = isodd(d) +_id(d) = div(d + 1, 2) + +# https://en.wikipedia.org/wiki/Chebyshev_polynomials#Properties +# Using +# sin(a + b) = sin(a) cos(b) + cos(a) sin(b) +# sin(a - b) = sin(a) cos(b) - cos(a) sin(b) +# If a > b +# sin(a) cos(b) = sin(a + b) + sin(a - b) +# sin(a) cos(b) = sin(a + b) + sin(a - b) +function univariate_mul!(::Mul{Trigonometric}, terms, var, a, b) + @assert !iszero(a) + @assert !iszero(b) + I = eachindex(terms) + da = _id(a) + db = _id(b) + for i in I + if _is_cos(a) == _is_cos(b) + # Chebyshev first kind + mono = MP.monomial(terms[i]) * var^(_cos_id(da + db)) + terms[i] = MA.mul!!(terms[i], var^_cos_id(abs(da - db))) + terms[i] = MA.operate!!(/, terms[i], 2) + α = MA.copy_if_mutable(MP.coefficient(terms[i])) + push!(terms, MP.term(α, mono)) + # cos(a + b) = cos(a) cos(b) - sin(a) sin(b) + # cos(a - b) = cos(a) cos(b) + sin(a) sin(b) + # cos(a - b) - cos(a + b) + if _is_sin(a) + terms[end] = MA.operate!!(*, terms[end], -1) + end + else + if _is_cos(a) + da, db = db, da + end + # sin(da) * cos(db) + if da == db + # sin(da) * cos(da) = sin(2da) / 2 + terms[i] = MA.mul!!(terms[i], var^_cos_id(da + db)) + terms[i] = MA.operate!!(/, terms[i], 2) + else + # Using + # sin(a + b) = sin(a) cos(b) + cos(a) sin(b) + # sin(a - b) = sin(a) cos(b) - cos(a) sin(b) + # If a > b + # sin(a) cos(b) = (sin(a + b) + sin(a - b)) / 2 + # If a < b + # sin(a) cos(b) = (sin(b + a) - sin(b - a)) / 2 + mono = MP.monomial(terms[i]) * var^(_sin_id(da + db)) + terms[i] = MA.mul!!(terms[i], var^_sin_id(abs(da - db))) + terms[i] = MA.operate!!(/, terms[i], 2) + α = MA.copy_if_mutable(MP.coefficient(terms[i])) + push!(terms, MP.term(α, mono)) + if da < db + terms[i] = MA.operate!!(*, terms[i], -1) + end + end + end + end + return +end + +function degree_one_univariate_polynomial(::Type{Trigonometric}, variable) + MA.@rewrite(variable + 0) +end + +function recurrence_eval( + ::Type{Trigonometric}, + previous::AbstractVector, + value, + degree, +) + d = _id(degree) + if _is_cos(degree) + # Chebyshev first order + return 2 * value * previous[_cos_id(d - 1) + 1] - previous[_cos_id(d - 2) + 1] + else + return sqrt(1 - previous[degree]^2) + end +end + +function _promote_coef(::Type{T}, ::Type{Trigonometric}) where {T} + return _promote_div(T) +end