Calculation of formal solution of dynamical systems

General Usage
21 
julia> using Nemo
julia> using BenchmarkTools

julia> T, z = PowerSeriesRing(QQ, 400, "z")
julia> f = sum([z^i*1//rand(1:10^5) for i=0:399]);
julia> gg=x->x^8;
julia> @btime gg($f);
  53.771 ms (2838 allocations: 13.35 MiB)

That’s faster than your mathematica timings by a factor of ~200. The construction code for f is pretty slow the way I wrote it, the right way is

julia> A=[1//rand(1:10^5) for n=0:399];
julia> f=T(A, 400, 400, 0);
julia> @btime gg($f);
  59.129 ms (2838 allocations: 14.00 MiB)
1 Like
22

I came up with the following idea:

mutable struct InductivePowerSeries
    coeffs::Array{Rational{BigInt},1}
    knownUpTo::Int64
    name::String
end


InductivePowerSeries(N::Int,s::String) = InductivePowerSeries(Array{Rational{BigInt},1}(N),0,s)

function getCoeff(s::InductivePowerSeries,n::Int)::Union{Rational{BigInt},Symbol}
    if n <= s.knownUpTo
        return s.coeffs[n]
    else
        return parse(s.name)
    end
end

function setNextCoeff(s::InductivePowerSeries,n::Int,value::Rational{BigInt})
    if n == s.knownUpTo + 1
        s.coeffs[n] = value
        s.knownUpTo += 1
        return
    else
        error("Error: the power series is known up to order $(s.knownUpTo) so the term at position $n cannot be set")
    end
end

setNextCoeff(s::InductivePowerSeries,n::Int,value::Union{Int,Rational{Int},BigInt}) = setNextCoeff(s,n,Rational{BigInt}(value))

Is there any obvious way that this code can be improved? (Before you complain, n is part of the arguments of setNextCoeff for sanity checking.)

This is halfway there as it returns the value when it exists and a symbol when it doesn’t. This means that I will eventually get something like 3+5*:x. How can I solve this?

Or better to my actual problem, if I have two such equations, say 3*:x+:y-1 and :x-2*:y+1, how can I tell julia to solve this system? It seems like something that AbstractAlgebra.jl can do easily, but I read a bit the documentation (admittedly not carefully enough) and I cannot figure this out.

23

Ok, I’m not sure whether anyone is actually interested in this kind of problems, but I thought I should post my solution just in case. The code is rather hacky but works better than what I expected.

I ended up using SymEngine.jl because I don’t know how else to produce the equations.

It calculates the first 100 coefficients in about 0.14s and the first 40 in less than 10s! That may not sound too impressive, but you need to keep in mind that the coefficients grow as factorial so the numerators and denominators are ridiculously large numbers. Anyway, here’s the code:

using SymEngine
using Memoize

mutable struct InductivePowerSeries
    coeffs::Array{Rational{BigInt},1}
    startingOrder::Int64
    knownUpTo::Int64
    symbol::SymEngine.Basic
    dummy::SymEngine.Basic
end

function InductivePowerSeries(N::Int,s::String)
    eval(parse(s*" = symbols(\""*s*"\")"))
    eval(parse(s*"0 = symbols(\""*s*"0\")"))
    return InductivePowerSeries(Array{Rational{BigInt},1}(N),1,0,eval(s),eval(s*"0"))
end

function getCoeff(s::InductivePowerSeries,n::Int)::Union{Rational{BigInt},SymEngine.Basic}
    if n < s.startingOrder
        return Rational{BigInt}(0)
    end
    if n <= s.knownUpTo
        return s.coeffs[n]
    end
    if n == s.knownUpTo + 1
        return s.symbol
    end
    if n == s.knownUpTo + 2
        return s.dummy
    end
    error("Error: the power series is known up to order $(s.knownUpTo) so the term at position $n cannot be returned.")
end

function setNextCoeff!(s::InductivePowerSeries,n::Int,value::Rational{BigInt})
    if n == s.knownUpTo + 1
        s.coeffs[n] = value
        s.knownUpTo += 1
        return
    else
        error("Error: the power series is known up to order $(s.knownUpTo) so the term at position $n cannot be set.")
    end
end

setNextCoeff!(s::InductivePowerSeries,n::Int,value::Union{Int,Rational{Int},BigInt}) = setNextCoeff!(s,n,Rational{BigInt}(value))


function getCoeffOfPower(s::InductivePowerSeries,p::Int,n::Int)::Union{Rational{BigInt},SymEngine.Basic}
    if n < p * s.startingOrder
        return Rational{BigInt}(0)
    end
    if n < p * s.startingOrder + s.knownUpTo
        return getCoeffOfPowerNumber(s,p,n)
    end
    if p == 1
        return getCoeff(s,n)
    end
    if n <= p + s.knownUpTo + 1
        sum = 0*s.symbol
        for i = s.startingOrder*(p-1):n-s.startingOrder
            sum += getCoeffOfPower(s,p-1,i)*getCoeff(s,n - i)
        end
        return sum
    end
    error("Error: the power series is known up to order $(s.knownUpTo) and the order $n of power $p was asked.")
end


@memoize function getCoeffOfPowerNumber(s::InductivePowerSeries,p::Int,n::Int)::Rational{BigInt}
    if p == 1
        return getCoeff(s,n)
    elseif p == 2
        sum = Rational{BigInt}(0)
        for i = s.startingOrder:n-s.startingOrder
            sum += getCoeff(s,i)*getCoeff(s,n-i)
        end
        return sum
    else
        sum = Rational{BigInt}(0)
        for i = s.startingOrder*(p-1):n-s.startingOrder
            sum += getCoeffOfPowerNumber(s,p-1,i)*getCoeff(s,n-i)
        end
        return sum
    end
end


function translationCoeff(p::Int,n::Int)::BigInt
    if (p + n) % 2 == 0
        return binomial(BigInt(n-1),p-1)
    else
        return -binomial(BigInt(n-1),p-1)
    end
end


function getTranslatedCoeffNumber(s::InductivePowerSeries,n::Int)::Rational{BigInt}
    sum = BigInt(0)
    for i = 1:n
        sum += translationCoeff(i,n)*getCoeff(s,i)
    end
    return sum
end

function getTranslatedCoeffButTheLastNumber(s::InductivePowerSeries,n::Int)::Rational{BigInt}
    sum = BigInt(0)
    for i = 1:n-1
        sum += translationCoeff(i,n)*getCoeff(s,i)
    end
    return sum
end

function getTranslatedCoeffButTheLast2Number(s::InductivePowerSeries,n::Int)::Rational{BigInt}
    sum = BigInt(0)
    for i = 1:n-2
        sum += translationCoeff(i,n)*getCoeff(s,i)
    end
    return sum
end

function getTranslatedCoeff(s::InductivePowerSeries,n::Int)::Union{Rational{BigInt},SymEngine.Basic}
    if n > s.knownUpTo + 2
        error("Error: the power series is known up to order $(s.knownUpTo) but the translation coefficient asked for is at order $n")
    end
    if n <= s.knownUpTo
        return getTranslatedCoeffNumber(s,n)
    end
    if n == s.knownUpTo + 1
        return s.symbol + getTranslatedCoeffButTheLastNumber(s,n)
    end
    if n == s.knownUpTo + 2
        return s.dummy-(n-1)*s.symbol + getTranslatedCoeffButTheLast2Number(s,n)
    end
end

function calculateCoeffs(Z::InductivePowerSeries,MaxDeg::Int)
    for n = Z.knownUpTo+2:MaxDeg + 1
        p = expand(getTranslatedCoeff(Z,n) - getCoeffOfPower(Z,1,n) - getCoeffOfPower(Z,2,n) - getCoeffOfPower(Z,3,n))
        value = -Rational{BigInt}(p(0))//Rational{BigInt}((p-p(0))(1))
        setNextCoeff!(Z,n-1,value)
    end
end

In case someone really cares to understand the code, I can try to explain it, but I would bet otherwise. Thank you all for the discussion, it was actually very useful.

2 Likes
24

Good to hear it. And the time does look good for a tool that’s meant to support general-purpose symbolic languages.

EDIT: I think there is a lot of interest in general for this kind of thing. If you were to write a few lines and posted it, it might turn out to be useful sometime.

25

Well, in that case I will do it properly. I will write the code to solve my actual problem and then I will upload it on GitHub (hopefully in a week or so)

1 Like
26

You can also use Reduce to define GetCoeff like this

julia> using Reduce; Algebra.operator(:z)

julia> function GetCoeff(p,n)
           if p == 1
               :(z(n))
           elseif p == 2
               Algebra.sum(:(z(i)*z($n-i)),:i,0,n)
           else
               Algebra.:+([Algebra.:*(GetCoeff(p-1,i),:(z($n-$i))) for i ∈ 0:n]...)
           end
       end
GetCoeff (generic function with 1 method)

julia> GetCoeff(3,5)
:(3 * (z(5) * z(0) ^ 2 + 2 * z(4) * z(1) * z(0) + 2 * z(3) * z(2) * z(0) + z(3) * z(1) ^ 2 + z(2) ^ 2 * z(1)))
1 Like