Seek efficient fast manipulation of power series

New to Julia
1

I am investigating changing from Reduce computer algebra to Julia symbolics. My main first need is efficient manipulation of power series (aiming for multiple small variables, and coefficients being complicated algebra and/or vector coefficients). How can it be done?

A simple example. In Reduce I simply do the following
1: on div,revpri;
2: let a^5=>0;
3: (1+a/2)^9;
1 + 9/2a + 9a2 + 21/2*a3 + 63/8a**4
4: ws(1-a/2);
1 + 4a + 27/4a2 + 6*a3 + 21/8*a**4
to indicate how all subsequent algebraic steps are automatically and seamlessly computed to truncations in “a” to error O(a^5).

In Julia, all I can see to do is
julia> using Symbolics
julia> @variables a
julia> taylor((1+a/2)^9,a,0:4)
(1//1) + (9//2)a + (9//1)(a^2) + (21//2)(a^3) + (63//8)(a^4)
julia> taylor(ans*(1-a/2),a,0:4)
(1//1) + (4//1)a + (27//4)(a^2) + (6//1)(a^3) + (21//8)(a^4)
Unfortunately, this code seems slow computationally because the taylor() functions seems slow (correct me if I am wrong). And also, it is slow for human coding because one has to repeatedly explicitly invoke the taylor() function (and also explicitly invoke expand() function for anything that does not need immediate truncation).

Please advise what could be efficient and fast in Julia.

Tony Roberts

2

A good way out is TaylorDiff.jl

You can simply define a as a TaylorScalar{5} (not sure it’s the exact right syntax, see their docs) and it will apply the faa di Bruno formula at compile time to expand f(a) for any Julia function f.

Edit:


a=0.367294
f(x) = (1+x/2)^9


using TaylorDiff
t = TaylorScalar{5}(a, 1.0)
f(t)

This might not work I did not test it but it should

3

Sorry I have no idea what the following results of this are. It seems to have nothing whatsoever to do with the power series expansion of (1+a/2)^9 to errors O(a^5).

4

Divide them by factorial.(0:5) and you’ll have the coefficient you want. Now that I think about it you might want to set a=0 instead or 0.36…

5

Not quite. All I need to do for this simple example is, after ‘using TaylorDiff’,

julia> a=TaylorScalar{5}(0.0,1.0)
TaylorScalar{Float64, 5}(0.0, (1.0, 0.0, 0.0, 0.0, 0.0))
julia> (1+a/2)^9
TaylorScalar{Float64, 5}(1.0, (4.5, 9.0, 10.5, 7.875, 3.9375))
julia> ans*(1-a/2)
TaylorScalar{Float64, 5}(1.0, (4.0, 6.75, 6.0, 2.625, 0.0))

This sets up a to be ‘variable’ a to be special in that it carries all its derivatives up to and including 5th order, but no more so the error here is O(a^6). Then subsequent expressions inherit all those derivatives and their respective computation. Very nice.
Three questions:

  • Although for physically complicated problems I will be happy with floating point coefficients, for many physically basic problems I want the coefficients as rational numbers, with numerator and denominator computed exactly to arbitrarily high precision. How can that be done?
  • Many times the coefficients are complicated algebraic expressions: in basic PDE examples coefficients might be each a finite sum of sin(n*x). Does TaylorDiff handle those?
  • This example indicate human-coding could be simple, but what about computational efficiency? How can I time a chunk of such code?

Thanks, Tony

1 Like
6

  1. Try with a=TaylorScalar{5}(0//1,1//1) ?

  2. don’t know, the easiest way to know would be trying

  3. the expansion happens at compile time so depending on your use case it is either the most efficient way of doing it or one of the worsts : if you want to reuse the same expansion over and over again on different values, then this is perfect. If on the other hand you want to expand only once per different expression, then this might not be the best code performance-wise

I suggests you read a bit about Julia, maybe the getting started page of the manual ? I haven’t seen that was your first post, welcome by the way !