# Calculation of formal solution of dynamical systems

**URL:** <https://discourse.julialang.org/t/calculation-of-formal-solution-of-dynamical-systems/11046>\
**Category:** General Usage\
**Tags:** question\
**Created:** [May 21, 2018, 2:40pm UTC](https://discourse.julialang.org/t/calculation-of-formal-solution-of-dynamical-systems/11046 "2018-05-21T14:40:17Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![gmouts](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gmouts/32/9896_2.png) [@gmouts](https://discourse.julialang.org/u/gmouts)\
**Post date:** [May 21, 2018, 2:40pm UTC](https://discourse.julialang.org/t/calculation-of-formal-solution-of-dynamical-systems/11046/1 "2018-05-21T14:40:17Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![foobar\_lv2](https://avatars.discourse-cdn.com/v4/letter/f/ee59a6/32.png) [@foobar\_lv2](https://discourse.julialang.org/u/foobar_lv2)\
**Post date:** [May 21, 2018, 2:57pm UTC](https://discourse.julialang.org/t/calculation-of-formal-solution-of-dynamical-systems/11046/2 "2018-05-21T14:57:09Z")

</div>

Have you looked at [nemo.jl](http://nemocas.github.io/Nemo.jl/latest/index.html)? I’d hope that it has all the goodies you need built in.

---

<div class="post-metadata">

**Author:** ![jlapeyre](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jlapeyre/32/4514_2.png) [@jlapeyre](https://discourse.julialang.org/u/jlapeyre)\
**Post date:** [May 21, 2018, 6:35pm UTC](https://discourse.julialang.org/t/calculation-of-formal-solution-of-dynamical-systems/11046/3 "2018-05-21T18:35:23Z")

</div>

Can you post an example of the Mathematica code ?

---

<div class="post-metadata">

**Author:** ![gmouts](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gmouts/32/9896_2.png) [@gmouts](https://discourse.julialang.org/u/gmouts)\
**Post date:** [May 22, 2018, 1:03am UTC](https://discourse.julialang.org/t/calculation-of-formal-solution-of-dynamical-systems/11046/4 "2018-05-22T01:03:22Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![gmouts](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gmouts/32/9896_2.png) [@gmouts](https://discourse.julialang.org/u/gmouts)\
**Post date:** [May 22, 2018, 1:03am UTC](https://discourse.julialang.org/t/calculation-of-formal-solution-of-dynamical-systems/11046/5 "2018-05-22T01:03:24Z")

</div>

I hadn’t before, I did now 🙂

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)

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [May 22, 2018, 6:38am UTC](https://discourse.julialang.org/t/calculation-of-formal-solution-of-dynamical-systems/11046/6 "2018-05-22T06:38:13Z")

</div>

This sounds like a job for [`TaylorSeries.jl`](https://github.com/JuliaDiff/TaylorSeries.jl).  
However, I’m not completely clear on what you would like to do.

cc @lbenet

---

<div class="post-metadata">

**Author:** ![jlapeyre](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jlapeyre/32/4514_2.png) [@jlapeyre](https://discourse.julialang.org/u/jlapeyre)\
**Post date:** [May 22, 2018, 10:21am UTC](https://discourse.julialang.org/t/calculation-of-formal-solution-of-dynamical-systems/11046/7 "2018-05-22T10:21:50Z")

</div>

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](https://github.com/JuliaAlgebra), as well as a few other polynomial-related packages.

@dpsanders suggested [TaylorSeries.jl](https://github.com/JuliaDiff/TaylorSeries.jl), which looks like a more likely candidate than the others.

Following is your Mathematica code translated into [Symata.jl](https://github.com/jlapeyre/Symata.jl). But, I don’t have access to Mathematica at the moment, so I can’t test it.

_getcoeff.sj_

```julia
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])))))

```

```julia
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](https://github.com/JuliaPy/SymPy.jl). Or maybe [SymEngine.jl](https://github.com/symengine/SymEngine.jl), but it is not well-developed compared to `SymPy`.

---

<div class="post-metadata">

**Author:** ![gmouts](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gmouts/32/9896_2.png) [@gmouts](https://discourse.julialang.org/u/gmouts)\
**Post date:** [May 22, 2018, 11:24am UTC](https://discourse.julialang.org/t/calculation-of-formal-solution-of-dynamical-systems/11046/8 "2018-05-22T11:24:15Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![gmouts](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gmouts/32/9896_2.png) [@gmouts](https://discourse.julialang.org/u/gmouts)\
**Post date:** [May 22, 2018, 11:24am UTC](https://discourse.julialang.org/t/calculation-of-formal-solution-of-dynamical-systems/11046/9 "2018-05-22T11:24:27Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![gmouts](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gmouts/32/9896_2.png) [@gmouts](https://discourse.julialang.org/u/gmouts)\
**Post date:** [May 22, 2018, 12:11pm UTC](https://discourse.julialang.org/t/calculation-of-formal-solution-of-dynamical-systems/11046/10 "2018-05-22T12:11:01Z")

</div>

@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?

---

<div class="post-metadata">

**Author:** ![kristoffer.carlsson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kristoffer.carlsson/32/22_2.png) [@kristoffer.carlsson](https://discourse.julialang.org/u/kristoffer.carlsson)\
**Post date:** [May 22, 2018, 12:31pm UTC](https://discourse.julialang.org/t/calculation-of-formal-solution-of-dynamical-systems/11046/11 "2018-05-22T12:31:47Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![gmouts](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gmouts/32/9896_2.png) [@gmouts](https://discourse.julialang.org/u/gmouts)\
**Post date:** [May 22, 2018, 1:22pm UTC](https://discourse.julialang.org/t/calculation-of-formal-solution-of-dynamical-systems/11046/12 "2018-05-22T13:22:15Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![jlapeyre](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jlapeyre/32/4514_2.png) [@jlapeyre](https://discourse.julialang.org/u/jlapeyre)\
**Post date:** [May 22, 2018, 1:26pm UTC](https://discourse.julialang.org/t/calculation-of-formal-solution-of-dynamical-systems/11046/13 "2018-05-22T13:26:35Z")

</div>

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:

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

```

---

<div class="post-metadata">

**Author:** ![foobar\_lv2](https://avatars.discourse-cdn.com/v4/letter/f/ee59a6/32.png) [@foobar\_lv2](https://discourse.julialang.org/u/foobar_lv2)\
**Post date:** [May 22, 2018, 2:59pm UTC](https://discourse.julialang.org/t/calculation-of-formal-solution-of-dynamical-systems/11046/14 "2018-05-22T14:59:22Z")

</div>

> [@gmouts](#):
>
> …and it is mindbogglingly fast! I expanded 1/(1−t)81/(1-t)^8 up to order 400 and it takes just less than a second.

Are you kidding me?

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

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 22, 2018, 3:10pm UTC](https://discourse.julialang.org/t/calculation-of-formal-solution-of-dynamical-systems/11046/15 "2018-05-22T15:10:33Z")

</div>

No one has mentioned Memoize.jl?

[https://github.com/simonster/Memoize.jl](https://github.com/simonster/Memoize.jl)

---

<div class="post-metadata">

**Author:** ![jlapeyre](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jlapeyre/32/4514_2.png) [@jlapeyre](https://discourse.julialang.org/u/jlapeyre)\
**Post date:** [May 22, 2018, 3:29pm UTC](https://discourse.julialang.org/t/calculation-of-formal-solution-of-dynamical-systems/11046/16 "2018-05-22T15:29:09Z")

</div>

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

```

---

<div class="post-metadata">

**Author:** ![jlapeyre](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jlapeyre/32/4514_2.png) [@jlapeyre](https://discourse.julialang.org/u/jlapeyre)\
**Post date:** [May 22, 2018, 3:32pm UTC](https://discourse.julialang.org/t/calculation-of-formal-solution-of-dynamical-systems/11046/17 "2018-05-22T15:32:27Z")

</div>

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 🙂

---

<div class="post-metadata">

**Author:** ![foobar\_lv2](https://avatars.discourse-cdn.com/v4/letter/f/ee59a6/32.png) [@foobar\_lv2](https://discourse.julialang.org/u/foobar_lv2)\
**Post date:** [May 22, 2018, 3:47pm UTC](https://discourse.julialang.org/t/calculation-of-formal-solution-of-dynamical-systems/11046/18 "2018-05-22T15:47:49Z")

</div>

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).

---

<div class="post-metadata">

**Author:** ![lbenet](https://avatars.discourse-cdn.com/v4/letter/l/35a633/32.png) [@lbenet](https://discourse.julialang.org/u/lbenet)\
**Post date:** [May 22, 2018, 6:45pm UTC](https://discourse.julialang.org/t/calculation-of-formal-solution-of-dynamical-systems/11046/19 "2018-05-22T18:45:36Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![gmouts](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gmouts/32/9896_2.png) [@gmouts](https://discourse.julialang.org/u/gmouts)\
**Post date:** [May 22, 2018, 7:13pm UTC](https://discourse.julialang.org/t/calculation-of-formal-solution-of-dynamical-systems/11046/20 "2018-05-22T19:13:22Z")

</div>

@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.

[Next page](https://discourse.julialang.org/t/calculation-of-formal-solution-of-dynamical-systems/11046.md?page=2)