From 560c6d830a8dfa47676671b858024e9c552bf843 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Tue, 23 Jul 2024 07:26:21 +0200 Subject: [PATCH] more Fourier operator functions (#280) * more methods for multiplication and division of Fourier operators and numbers * add tests * operations with UniformScaling --- Project.toml | 2 +- src/fourier_operators.jl | 88 ++++++++++++++++++++++++++++++++++ test/fourier_operators_test.jl | 55 ++++++++++++++------- 3 files changed, 128 insertions(+), 17 deletions(-) diff --git a/Project.toml b/Project.toml index 894d4335..d9720300 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "SummationByPartsOperators" uuid = "9f78cca6-572e-554e-b819-917d2f1cf240" author = ["Hendrik Ranocha"] -version = "0.5.64" +version = "0.5.65" [deps] ArgCheck = "dce04be8-c92d-5529-be00-80e4d2c0e197" diff --git a/src/fourier_operators.jl b/src/fourier_operators.jl index dbef6028..28c78d1f 100644 --- a/src/fourier_operators.jl +++ b/src/fourier_operators.jl @@ -530,6 +530,23 @@ function Base.:+(rat1::Union{FourierDerivativeOperator,FourierPolynomialDerivati FourierRationalDerivativeOperator(rat1) + rat2 end +function Base.:+(rat::FourierRationalDerivativeOperator, scaling::UniformScaling) + rat + scaling.λ * rat.D1^0 +end + +function Base.:+(scaling::UniformScaling, rat::FourierRationalDerivativeOperator) + rat + scaling +end + +function Base.:-(rat::FourierRationalDerivativeOperator) + @unpack num_coef = rat + for idx in eachindex(num_coef) + num_coef = Base.setindex(num_coef, -num_coef[idx], idx) + end + + FourierRationalDerivativeOperator(rat.D1, num_coef, rat.den_coef) +end + function Base.:-(rat1::FourierRationalDerivativeOperator, rat2::FourierRationalDerivativeOperator) T = eltype(rat1) @argcheck T == eltype(rat2) ArgumentError @@ -551,6 +568,14 @@ function Base.:-(rat1::Union{FourierDerivativeOperator,FourierPolynomialDerivati FourierRationalDerivativeOperator(rat1) - rat2 end +function Base.:-(rat::FourierRationalDerivativeOperator, scaling::UniformScaling) + rat - scaling.λ * rat.D1^0 +end + +function Base.:-(scaling::UniformScaling, rat::FourierRationalDerivativeOperator) + scaling + (-rat) +end + function Base.:*(rat1::FourierRationalDerivativeOperator, rat2::FourierRationalDerivativeOperator) T = eltype(rat1) @argcheck T == eltype(rat2) ArgumentError @@ -572,6 +597,27 @@ function Base.:*(rat1::Union{FourierDerivativeOperator,FourierPolynomialDerivati FourierRationalDerivativeOperator(rat1) * rat2 end +function Base.:*(factor::Union{Real,Integer}, rat::FourierRationalDerivativeOperator) + @unpack num_coef = rat + for idx in eachindex(num_coef) + num_coef = Base.setindex(num_coef, factor * num_coef[idx], idx) + end + + FourierRationalDerivativeOperator(rat.D1, num_coef, rat.den_coef) +end + +function Base.:*(poly::FourierRationalDerivativeOperator, factor::Union{Real,Integer}) + factor * poly +end + +function Base.:*(rat::FourierRationalDerivativeOperator, scaling::UniformScaling) + scaling.λ * rat +end + +function Base.:*(scaling::UniformScaling, rat::FourierRationalDerivativeOperator) + rat * scaling +end + function Base.:/(rat1::FourierRationalDerivativeOperator, rat2::FourierRationalDerivativeOperator) T = eltype(rat1) @argcheck T == eltype(rat2) ArgumentError @@ -593,6 +639,48 @@ function Base.:/(rat1::Union{FourierDerivativeOperator,FourierPolynomialDerivati FourierRationalDerivativeOperator(rat1) / rat2 end +function Base.:/(factor::Union{Real,Integer}, rat::FourierRationalDerivativeOperator) + @unpack den_coef = rat + for idx in eachindex(den_coef) + den_coef = Base.setindex(den_coef, factor * den_coef[idx], idx) + end + + FourierRationalDerivativeOperator(rat.D1, den_coef, rat.num_coef) +end + +function Base.:/(rat::FourierRationalDerivativeOperator, factor::Union{Real,Integer}) + @unpack num_coef = rat + for idx in eachindex(num_coef) + num_coef = Base.setindex(num_coef, num_coef[idx] / factor, idx) + end + + FourierRationalDerivativeOperator(rat.D1, num_coef, rat.den_coef) +end + +function Base.:\(factor::Union{Real,Integer}, rat::FourierRationalDerivativeOperator) + rat / factor +end + +function Base.:\(rat::FourierRationalDerivativeOperator, factor::Union{Real,Integer}) + inv(rat) * factor +end + +function Base.:/(rat::FourierRationalDerivativeOperator, scaling::UniformScaling) + rat / scaling.λ +end + +function Base.:/(scaling::UniformScaling, rat::FourierRationalDerivativeOperator) + scaling.λ / rat +end + +function Base.:\(rat::FourierRationalDerivativeOperator, scaling::UniformScaling) + rat \ scaling.λ +end + +function Base.:\(scaling::UniformScaling, rat::FourierRationalDerivativeOperator) + scaling.λ \ rat +end + function mul!(dest::AbstractVector{T}, rat::FourierRationalDerivativeOperator, u::AbstractVector{T}) where {T} @unpack D1, num_coef, num_coef_end, den_coef, den_coef_end = rat diff --git a/test/fourier_operators_test.jl b/test/fourier_operators_test.jl index 7c5c9aec..697392ba 100644 --- a/test/fourier_operators_test.jl +++ b/test/fourier_operators_test.jl @@ -54,13 +54,13 @@ for T in (Float32, Float64) end -# Check Fourier polynomial operators +# Check Fourier polynomial/rational operators for T in (Float32, Float64) xmin = -one(T) xmax = one(T) for N in (8, 9) - D = fourier_derivative_operator(xmin, xmax, N) - x = grid(D) + D = @inferred fourier_derivative_operator(xmin, xmax, N) + x = @inferred grid(D) u = @. sinpi(x) - cospi(x)^2 + exp(sinpi(x)) println(devnull, D) @@ -72,37 +72,60 @@ for T in (Float32, Float64) @test !issymmetric(I + D) @test !issymmetric(D - I) - @test SummationByPartsOperators.xmin(D^2) ≈ xmin - @test SummationByPartsOperators.xmax(D^2) ≈ xmax + @test @inferred(SummationByPartsOperators.xmin(D^2)) ≈ xmin + @test @inferred(SummationByPartsOperators.xmax(D^2)) ≈ xmax - poly = (I + 2D + 5*D^2) * (2I * D - D^3 * 5I) * (D*2 - D^2 * 5) + poly = @inferred (I + 2D + 5*D^2) * (2I * D - D^3 * 5I) * (D*2 - D^2 * 5) @test poly.coef == (0.0, 0.0, 4.0, -2.0, -10.0, -45.0, 0.0, 125.0) println(devnull, poly) - @test (I + one(T)/2*D) * u ≈ (u + D*u ./ 2) + @test @inferred(I + one(T)/2*D) * u ≈ (u + D*u ./ 2) v = (I - D^2) * u @test inv(I - D^2) * v ≈ u - @test SummationByPartsOperators.xmin(inv(I - D^2)) ≈ xmin - @test SummationByPartsOperators.xmax(inv(I - D^2)) ≈ xmax + @test @inferred(SummationByPartsOperators.xmin(inv(I - D^2))) ≈ xmin + @test @inferred(SummationByPartsOperators.xmax(inv(I - D^2))) ≈ xmax - v = (I - D^2) \ u + v = @inferred(I - D^2) \ u @test D * v ≈ (D / (I - D^2)) * u - rat = (I - D^2) / (I + D^4) + rat = @inferred((I - D^2) / (I + D^4)) println(devnull, rat) v = rat * u @test (I - D^2) \ (v + D^4 * v) ≈ u - rat1 = (I - D^2) / (I + D^4) - rat2 = (I + D^4) / (I - D^2) - rat3 = (I - D^4) / (I - D^2) - @test (rat2 + rat3) * u ≈ 2 * ((I - D^2) \ u) - @test (rat1 * rat2) * u ≈ u + rat1 = @inferred((I - D^2) / (I + D^4)) + rat2 = @inferred((I + D^4) / (I - D^2)) + rat3 = @inferred((I - D^4) / (I - D^2)) + @test @inferred(rat2 + rat3) * u ≈ 2 * ((I - D^2) \ u) + @test @inferred(rat1 * rat2) * u ≈ u @test integrate(u, D) ≈ sum(mass_matrix(D) * u) @test integrate(u->u^2, u, D) ≈ dot(u, mass_matrix(D), u) + + # combine rational operators and scalars + @test @inferred(2 * rat1) * u ≈ rat1 * (2 * u) + @test @inferred(rat1 * 2) * u ≈ rat1 * (2 * u) + @test @inferred(2 / rat1) * u ≈ inv(rat1) * (2 * u) + @test @inferred(rat1 / 2) * u ≈ rat1 * (u / 2) + @test @inferred(2 \ rat1) * u ≈ rat1 * (u / 2) + @test @inferred(rat1 \ 2) * u ≈ inv(rat1) * (2 * u) + + # combine rational operators and uniform scaling via *, /, \ + @test @inferred((2 * I) * rat1) == 2 * rat1 + @test @inferred(rat1 * (2 * I)) == 2 * rat1 + @test @inferred((2 * I) / rat1) == 2 / rat1 + @test @inferred(rat1 / (2 * I)) == rat1 / 2 + @test @inferred((2 * I) \ rat1) == rat1 / 2 + @test @inferred(rat1 \ (2 * I)) == 2 / rat1 + + # combine rational operators and uniform scaling via +, - + @test @inferred(I + rat1) == D^0 + rat1 + @test @inferred(rat1 + I) == D^0 + rat1 + @test @inferred(I - rat1) == D^0 - rat1 + @test @inferred(rat1 - I) == rat1 - D^0 + @test @inferred(-rat1) == @inferred(-1 * rat1) end end