Calculation of formal solution of dynamical systems

General Usage
1

I want to calculate the formal solution of a discreet dynamical system (obviously not the whole infinite series, but let’s say all terms up to order N).

For simplicity let’s get the example of f(z) = z + z^2 + z^3. I want to find the formal solution of the equations z(t+1) = f(z(t)). I know that there is a formal solution of the form z(t) = -1/t + O(1/t^3) with rational coefficients, so I want to get the coefficients of the series.

I have done this in the past using mathematica. The way I did that is that I defined a function, GetCoeff, that takes an integer n, a power series z and a monomial x^n, and returns the n-th coefficient of z(t)^n. This was done recursively, i.e. it calls itself with lower degree monomials and I used memoization (which improved the time by a factor of 10 or so).

I would like to do the same with Julia, but I am not sure exactly how. Is there a module that I can use for power series? I would like to avoid reinventing the wheel. Can I define a function like GetCoeff?

2

Have you looked at nemo.jl? I’d hope that it has all the goodies you need built in.

3

Can you post an example of the Mathematica code ?

4

The function that I care about is the following:

GetCoeff[p_, n_] := GetCoeff[p, n] = Which[
   p == 1, z[n],
   p == 2, Sum[z[i] z[n - i], {i, 0, n}],
   True, Sum[GetCoeff[p - 1, i]  z[n - i], {i, 0, n}]
   ]

This assumes that the coefficients are stored using the variable z, with z[n] being the n-th coefficient. This is convenient in mathematica because everything is treated like a symbol.

So in practice I avoid using the power series as a power series, because this is too slow. So I am not sure what type I should use to store my data. Is it a fixed length array better maybe?

Edit: I should also say that because mathematica treats everything as symbols, asking for a value that is not there yet returns z[n] that can be used as an unknown in an equation. This is very convenient in this situation.

5

I hadn’t before, I did now :slight_smile:

I am afraid that it is going to be slow as I don’t really want to substitute a power series in a polynomial. This is too slow. I try to explain this in the post above (a reply to jlapeyre)

6

This sounds like a job for TaylorSeries.jl.
However, I’m not completely clear on what you would like to do.

cc @lbenet

7

It’s still not clear to me what you are looking for. You say

GetCoeff, that takes an integer n, a power series z and a monomial x^n,

But, the function you posted GetCoeff takes two arguments, not three; and they both appear to be integers, rather than a power series or monomial. You refer once to x; in x^n. You refer to a fixed polynomial z + z^2 + z^3. The function GetCoeff returns a polynomial in z[i] for i \in [0,1,...n]. Do these refer to the same z ? Perhaps i represents integer-valued t ? It would help if you explain things more precisely.

In general, these algebra tools involve a trade-off between flexibility-expressiveness on one hand and efficiency on the other. As you noted, Mathematica is more expressive than efficient. (Although a few decades and a few hundred million dollars can go a long way to closing this gap.) I’m not sure, but it looks like Nemo.jl and AbstractAlgebra.jl can’t do what you want easily; they want numeric coefficients. But, maybe asking the developers is worthwhile. There is also JuliaAlgebra, as well as a few other polynomial-related packages.

@dpsanders suggested TaylorSeries.jl, which looks like a more likely candidate than the others.

Following is your Mathematica code translated into Symata.jl. But, I don’t have access to Mathematica at the moment, so I can’t test it.

getcoeff.sj

GetCoeff(p_, n_) := GetCoeff(p, n) = Expand(If(
   p == 1, z(n),
   If(p == 2, Sum(z(i) * z(n - i), [i, 0, n]),
    Apply(Plus, Table(GetCoeff(p - 1, i) * z(n - i), [i, 0, n])))))

symata 12> Get("getcoeff.sj")

symata 12> GetCoeff(3,5)
Out(12) = 7z(0)*z(2)*z(3) + 7z(0)*z(1)*z(4) + 7*z(0)^2*z(5)

symata 13> ? GetCoeff
GetCoeff(2,0)=z(0)^2
GetCoeff(2,1)=(2z(0)*z(1))
GetCoeff(2,2)=(3z(0)*z(2))
GetCoeff(2,3)=(4z(0)*z(3))
GetCoeff(2,4)=(5z(0)*z(4))
GetCoeff(2,5)=(6z(0)*z(5))
GetCoeff(3,5)=(7z(0)*z(2)*z(3) + 7z(0)*z(1)*z(4) + 7*z(0)^2*z(5))
GetCoeff(p_,n_):=GetCoeff(p,n)=CompoundExpression(nothing,Expand(If(p == 1,z(n),If(p == 2,Sum((z(i)*z(n - i)),[i,0,n]),Apply(Plus,Table((GetCoeff((p - 1),i)*z(n - i)),[i,0,n]))))))

You can probably gain some efficiency by using Symata from the Julia side so that you can control the evaluation. Or, you can work directly in SymPy.jl. Or maybe SymEngine.jl, but it is not well-developed compared to SymPy.

8

Sorry for the confusion. I thought it would be clearer if I omitted the power series as an argument since in this case we have only one to worry about. The function I provided takes the exponent of the monomial p and n denotes that we care about the n-th term of the series.

So, what I am actually doing is the following:

(Still using f(z) = z + z^2 + z^3 as an example)

I have a function that gets f and returns GetCoeff[1,n] + GetCoeff[2,n] + GetCoeff[3,n] (this function is a completely different problem, but I will deal with it later)

And a function TransCoeff that gives the n-th coefficient of the series z(t+1) (this is a pretty straightforward thing to write)

Let’s say that we know all the coefficients up to order 9.

Then I evaluate the expression

TransCoeff[10] - (GetCoeff[1,10] + GetCoeff[2,10] + GetCoeff[3,10])

This evaluates to a*z[10] + b with a and b rational numbers. Then by solving the linear equation, I get the coefficient z[10] (this equation is linear because the series does not have a constant term, but this doesn’t matter in the present discussion)

The advantage that this has, is that in order to evaluate GetCoeff[3,10] we need the values of GetCoeff[2,i] with 1 \le i \le 9 . But these computations have been done previously, so memoization saves the day here. This becomes more prominent as the degree of f increases.

I want to avoid expanding polynomials as this is ends up being very redundant. My problem is actually the following. I could possibly use an array to store coefficients, but I need to “access” a coefficient before it is set. Mathematica has not problems with that, but how do I do it efficiently in julia?

@jlapeyre I would like to use native julia code, because I want to make it as fast as possible.

@dpsanders I haven’t seen this package before, I am looking at it now. Does the package avoid the redundancy I talked about above?

1 Like
9

I am afraid that my question has actually changed as I am trying to explain myself. I am sorry for this.

10

@dpsanders I played a bit with TaylorSeries.jl with BigInt rationals and it is mindbogglingly fast! I expanded 1/(1-t)^8 up to order 400 and it takes just less than a second. How is this even possible?

Regardless, my problem remains. I need to “access” coefficients that are not there yet and set them inductively? Is there any way to do that using TaylorSeries.jl?

11

So in short, you need something that would help you conveniently memoize function results? Or is that not a correct summary?

12

No, that’s not it. I know how to memoize a function. My problem is how to copy mathematica’s behaviour. For example:

Let’s say that we know z[1] = 1 and z[2] = -1, but we have not compute any other coefficient yet. Then mathematica evaluates:
GetCoeff[1,1] to 1
GetCoeff[1,2] to -1
GetCoeff[1,3] to z[3]

This is important because I can then use z[3] as an unknown in a linear equation. This is what I don’t know how to do in julia.

13

That is impressively fast. Another question is “why is Mathematica so slow ?” One reason is that it has to allow that any of your symbols might behave like this:

symata 1> n = 0
Out(1) = 0

symata 2> a :=  If(n < 2, (n += 1, 3), 0)

symata 2> x^a
Out(2) = x^3

symata 3> x^a
Out(3) = x^3

symata 4> x^a
Out(4) = 1
14

Are you kidding me?

julia> using Nemo
julia> using BenchmarkTools

julia> T, z = PowerSeriesRing(QQ, 400, "z")
julia> @btime divexact(1, (1-z)^8);
  28.043 ��s (90 allocations: 16.48 KiB)
julia> @btime divexact(1, (1-z+z^2 -z^17)^8);
  305.701 ��s (1108 allocations: 47.40 KiB)

julia> T, z = PowerSeriesRing(QQ, 4000, "z");
julia> @btime divexact(1, (1-z)^8);
  1.592 ms (9105 allocations: 324.85 KiB)
julia> @btime divexact(1, (1-z+z^2-z^17)^8);
  10.811 ms (12366 allocations: 2.38 MiB)

julia> T, z = PowerSeriesRing(QQ, 20, "z");
julia> divexact(1, 1-z+z^2)
1+z-z^3-z^4+z^6+z^7-z^9-z^10+z^12+z^13-z^15-z^16+z^18+z^19+O(z^20)
1 Like
15

No one has mentioned Memoize.jl?

https://github.com/simonster/Memoize.jl

2 Likes
16

I can’t get the binary dependencies for Nemo.jl to build in my v0.6.2 setup. So I have to settle for pure Julia, which is slower:

julia> using AbstractAlgebra;

julia> T, z = PowerSeriesRing(QQ, 4000, "z");

julia> using BenchmarkTools;

julia> @btime  divexact(1,(1-z)^8);
  103.203 ms (2043063 allocations: 50.94 MiB)
1 Like
17

I think that’s what @kristoffer.carlsson was getting at. But, caching results is only part of the problem. This does look like a fun problem, though. Wish I had fewer chores right now :slight_smile:

18

Wow, that’s still kinda OK. I mean, a factor of ~100 for nemo instead of abstractalgebra is nothing to sniff at, but abstractalgebra.jl still appears to beat the TaylorSeries.jl variant by a factor of ~100 (when comparing to OPs posted time; and it appears that whatever mathematica code OP used before was even slower).

19

Hi,

If I understand correctly the description of your problem, you are looking to build some sort of Taylor expansion of an invariant for the dynamical system z(t+1) = f(z). Is this interpretation correct? Then, shouldn’t you only get the periodic orbits z_0 = 0, z_0=-1, so no t dependence then?

Regarding the comparison of TaylorSeries.jl with AbstractAlgebra.jl, indeed, TaylorSeries.jl is slower. Yet, there are tricks you can exploit in TaylorSeries.jl to have some sort of unknown variable, which is what eventually leads to the linear equation you solve. Essentially, you deal with a polynomial in t and a, where a is a variable representing the unknown coefficient you are after, which you propagate under the algebra. I guess all this you can do it in AbstractAlgebra.jl as well.

20

@lbenet

Hi,

you got the idea. However, the fixed point of the system is parabolic. This means that the linearization is trivial and does not define the dynamics in any neighbourhood of the fixed point. The goal is to measure things about the separatrices of the fixed point. It is an easy observation that this system admits a formal solution of the form I wrote in the beginning. This solution however is actually formal, i.e. the radius of convergence is 0. But we can still use a truncation of it to get information about the separatrices, hence I care about the coefficients.

What you suggest is a good idea. I didn’t think of that before. But trying to explain what I need, I got an idea that I will try first.

In a more general note I thought I should try and see how fast mathematica actually is in this case. So this is what I got:

The command Timing[Series[1/(1 - t)^8, {t, 0, 400}]] evaluates in 0.02s, which is really good and much better than what I’ve experienced. So I assume that there is some kind of optimization going on, because the Taylor series of (1-t)^{-1} is very simple.

I fix n \in\{ 1,2,3,4,5\} and define

f[t_] := Evaluate[Sum[t^i/(RandomInteger[10^n] + 1), {i, 0, 400}]]

Then I evaluate

Timing[Series[f[t]^8, {t, 0, 400}]]

and I get on average 0.46s when n=1, 1.7s when n=2, 4.5s when n=3, 8.7s when n=4 and 13.4s when n=5. So my code was extremely slow when I tried to do that with the actual coefficients.

I would like to do the same test julia, but honestly I have no clue how.