TypedPolynomials.jl icon indicating copy to clipboard operation
TypedPolynomials.jl copied to clipboard

Antidifferentiation

Open bmxam opened this issue 2 years ago • 1 comments

Hi, I was looking for a way to symbolically integrate a MultiVariate polynomial, or at least to find a primitive. I ended up writing the below code to find it. Would you be interested by a PR? If yes, I could add a little bit of documentation, write unit tests etc. Note : as primitive is already a specific julia language word, I opted out for antidifferentiate.

using TypedPolynomials
using MultivariatePolynomials

antidifferentiate(p::Polynomial, v::Variable) = Polynomial(antidifferentiate.(terms(p), v))
antidifferentiate(t::Term, v::Variable) = coefficient(t) * antidifferentiate(monomial(t), v)

function antidifferentiate(m::Monomial{Vars,N}, v::V) where {Vars,N,V<:Variable}
    if TypedPolynomials.inmonomial(v, Vars...)
        return _antidiff(m, exponents(m), v)
    else
        # Insert `v` in the monomial
        return m * v
    end
end

function _antidiff(::Monomial{Vars}, exponents::NTuple{N,Integer}, v::Variable) where {Vars,N}
    vi = TypedPolynomials.find_variable_index(v, Vars)
    new_m = Monomial{Vars,N}(
        ntuple(i -> begin
                if i == vi
                    (exponents[i] == 0) ? 1 : exponents[i] + 1
                else
                    exponents[i]
                end
            end,
            Val{N}()
        )
    )
    return 1 / (exponents[vi] + 1) * new_m
end

@polyvar x y # assigns x (resp. y) to a variable of name x (resp. y)

p = 2x + 3.0x * y^2 + y + 1.5
@show p
@show antidifferentiate(p, x)
@show antidifferentiate(p, y)

p = 2y + 3
@show p
@show antidifferentiate(p, x)
@show antidifferentiate(p, y)

bmxam avatar Dec 04 '22 22:12 bmxam

Yes, I would be interested by a PR. Could you define the part for polynomial and terms in https://github.com/JuliaAlgebra/MultivariatePolynomials.jl like it's done for differentiate ? Also, the function antidifferentiate should be defined by MultivariatePolynomials so that it's part of the API

blegat avatar Dec 05 '22 07:12 blegat