Need an efficient way to truncate power series on-the-fly

[ProfAJRoberts]

My core research computes asymptotic series with complicated coefficients, such as the simple example u=asin(x)+a^3/12sin(3*x)+O(a^4), and sometimes with numerical matrix coefficients (up to thousands of elements), sometimes mixed numerics-algebra coefficients, and sometimes constructed to high-order (20th--50th). The first hurdle is to efficiently truncate algebra done with such asymptotic power series. I find that using the Symbolics.taylor() function is far too slow. Is there any tweak to Symbolics to truncate asymptotics efficiently?

For the very simplest example, the following code in Julia.Symbolics takes about one second to find the series solution of x(a) to 20th order, whereas the corresponding code in Reduce-algebra takes just 2 ms to find it to 31st order.

@variables a
oErr=20       # >20 fails on factorial(21)
x=1           # set initial approx to root
for iter=1:99 # iterative refinement
    global x,Res
    Res=taylor(x^2+a*x-1,a,0:oErr)  # current residual of quadratic
    x=expand(x-Res/2)               # correct approx using residual
    if isequal(Res,0); println("Success: iter=",iter); break; end
end      

[ufechner7]

Can you put the code example in triple backticks?

[ChrisRackauckas]

It might be good to revisit this after the SymbolicUtils v4 overhaul.

[hersle]

The current implementation (added by me in #1298) is very barebones and far from optimal. I think the best approach for improvements would be make Symbolics composable with a special-purpose Taylor series package like TaylorSeries.jl. See e.g. JuliaDiff/TaylorSeries.jl#242.

[ProfAJRoberts]

I have no idea what the following means "The current implementation (added by me in #1298) is very barebones and far from optimal."

TaylorSeries works and is comparable to Reduce for the very simplest example: for the following btime gives 10ms versus Reduce 2ms: ''' using TaylorSeries
set_taylor1_varname("a")
a=Taylor1(Rational{BigInt},32) # truncate power series
#a=set_variables(Rational{BigInt},"a",order=9) # fails in the arithmetic??
x=1 # set initial approx to root
for iter=1:99 # iterative refinement
global x
Res=x^2+a*x-1 # current residual of quadratic
x=x-Res/2 # correct approx using residual
if Res==0; println("Success: iter=",iter); break; end
end
'''
Although why set_variables() fails is completely mysterious.

So the next challenge is to manipulate expressions like u=asin(x)+a^3/12sin(3*x)+O(a^4) How is this to be done?