 # Calculation of formal solution of dynamical systems

``````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

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.

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

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.

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

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