PerezHz / TaylorIntegration.jl

ODE integration using Taylor's method, and more, in Julia
Other
126 stars 22 forks source link

Problem using macro with `pow!`, for zero Taylor1 series #99

Closed lbenet closed 4 years ago

lbenet commented 4 years ago

There is a subtle problem when using the macro in functions involving ^r (for r different from 0, 1, 2). I illustrate the problem here with the integration of the Duffing equation, using of the macro (using current master of TaylorSeries):

julia> using TaylorIntegration

julia> @taylorize function duffing!(du, u, params, t)
           du[1] = u[2]
           du[2] = u[1]-0.01*u[1]^3
           nothing
       end

julia> x = taylorinteg(duffing!, [0.0, 1.0], 0.0:0.25:20.0, 20, 1.e-20);

julia> x[2,:]
2-element Array{Float64,1}:
 NaN
 NaN

The problem is related with TaylorSeries.pow! and the fact that it requires one parameter (the order of the first non-zero entry of a polynomial), which by default is set to zero. That default is used by @taylorize, which produces the NaNs in this case, and it does not appear if one computes c = a^3.

julia> a = 0.0*Taylor1(20)   # Taylor polynomial identical to zero
 0.0 + 𝒪(t²¹)

julia> c = Taylor1(20)
 1.0 t + 𝒪(t²¹)

# We compute c=a^3 term by term, using pow!
julia> TaylorSeries.pow!(c, a, 3, 0)

julia> TaylorSeries.pow!(c, a, 3, 1)

julia> c
 NaN t + 𝒪(t²¹)

The problem is solved by https://github.com/JuliaDiff/TaylorSeries.jl/pull/236; we will require a new tag of TaylorSeries.jl and related updates here.

mforets commented 4 years ago

The problem is solved by JuliaDiff/TaylorSeries.jl#236; we will require a new tag of TaylorSeries.jl and related updates here.

Very nice! Strange that this didn't pop up earlier.

Incidentally, is there any technical difference between using u[1]-0.01*u[1]^3 or u[1]-(0.01*(u[1]*u[1])*u[1])? (I mean, performance-wise, by using parentheses to force binary operations).

lbenet commented 4 years ago

Very nice! Strange that this didn't pop up earlier.

We spotted this a week or two ago, precisely playing with the Duffing equation.

Incidentally, is there any technical difference between using u[1]-0.01u[1]^3 or u[1]-(0.01(u[1]u[1])u[1])? (I mean, performance-wise, by using parentheses to force binary operations).

Yes, there is a performance difference. The complexity of the calculation is what matters (for performance). One product of two series has "complexity 1", as well as the power. So exchanging a^3 by (a*a)*a does have a performance penalty. I think there may e precision difference too.

lbenet commented 4 years ago

Fixed by https://github.com/JuliaDiff/TaylorSeries.jl/pull/236 and #100