JuliaIntervals / TaylorModels.jl

Rigorous function approximation using Taylor models in Julia
Other
63 stars 14 forks source link

Add example evaluating rational function #99

Open mforets opened 3 years ago

mforets commented 3 years ago

This documentation PR is a continuation of https://github.com/JuliaIntervals/TaylorModels.jl/pull/98 and can be reviewed independently. Here I considered the evaluation of a rational function from this discourse question and merged ideas suggested by @stevengj and @dpsanders.


In the last paragraph, num(tm) and denom(tm) are Taylor models for the numerator and denominator of the rational function g. Evaluating their division fails (see stacktrace below). I think that's due to the absence of a zero order term in the denominator. Do you have any suggestion to approximate g with a single Taylor model in this case?

julia> num(tm) / denom(tm)

ArgumentError: Division does not define a Taylor1 polynomial;
order k=0 => coeff[0]=∅.

Stacktrace:
 [1] divfactorization at /home/mforets/.julia/packages/TaylorSeries/4PYzx/src/arithmetic.jl:462 [inlined]
 [2] /(::Taylor1{Interval{Float64}}, ::Taylor1{Interval{Float64}}) at /home/mforets/.julia/packages/TaylorSeries/4PYzx/src/arithmetic.jl:421
 [3] inv at ./number.jl:199 [inlined]
 [4] _rpa(::Type{TaylorModel1}, ::typeof(inv), ::Interval{Float64}, ::Interval{Float64}, ::Int64) at /home/mforets/.julia/dev/TaylorModels/src/rpa_functions.jl:26
 [5] rpa(::Function, ::TaylorModel1{Interval{Float64},Float64}) at /home/mforets/.julia/dev/TaylorModels/src/rpa_functions.jl:105
 [6] inv at /home/mforets/.julia/dev/TaylorModels/src/rpa_functions.jl:284 [inlined]
 [7] basediv at /home/mforets/.julia/dev/TaylorModels/src/arithmetic.jl:48 [inlined]
 [8] /(::TaylorModel1{Interval{Float64},Float64}, ::TaylorModel1{Interval{Float64},Float64}) at /home/mforets/.julia/dev/TaylorModels/src/arithmetic.jl:171
 [9] top-level scope at In[48]:1
 [10] include_string(::Function, ::Module, ::String, ::String) at ./loading.jl:1091
lbenet commented 3 years ago

Sorry for responding with such a loooong delay...

You are right, the problem of using TaylorModel1 is related to the fact that there is no zero-order term in the denominator (and numerator). Notice that this problem does not arise for Taylor1. The trick is related to the fact that division is overloaded for Taylor1 in such a way that it factors out the common powers (in the independent variable) in the numerator and denominator, which in this case are both of order 2, and then the answer is a constant (1 in this example). I should look up the thesis of Mioara Joldes, but if I recall correctly, this can't be done with TaylorModel1s, simply because the remainder is absolute.

Yet, this is actually possible with RTaylorModel1s, which is the type of example motivating their introduction. The difference with TaylorModel1 is that the remainder includes an appropriate power of the independent variable, and in that sense is relative. This subtlety allows to factorize out the common powers, like in this case.

The example is really a nice application of the relative Taylor models:

julia> using TaylorModels

julia> nn(x) = 2*(1 - x*exp(-x)-exp(-x))
nn (generic function with 1 method)

julia> dd(x) = (1-exp(-x))^2
dd (generic function with 1 method)

julia> f(x) = nn(x) / dd(x)
f (generic function with 1 method)

julia> rtm = RTaylorModel1(10, 0.0, -0.125 .. 0.125)
 1.0 t + [0, 0] t¹¹

julia> nn(rtm)
 1.0 t² - 0.6666666666666667 t³ + 0.25 t⁴ - 0.06666666666666667 t⁵ + 0.013888888888888888 t⁶ - 0.002380952380952381 t⁷ + 0.00034722222222222224 t⁸ - 4.409171075837743e-5 t⁹ + 4.96031746031746e-6 t¹⁰ + [-9.13422e-06, 4.21723e-06] t¹¹

julia> dd(rtm)
 1.0 t² - 1.0 t³ + 0.5833333333333333 t⁴ - 0.24999999999999997 t⁵ + 0.08611111111111111 t⁶ - 0.024999999999999998 t⁷ + 0.006299603174603175 t⁸ - 0.0014054232804232803 t⁹ + 0.00028163580246913576 t¹⁰ + [-5.32954e-05, -4.91783e-05] t¹¹

julia> frtm = f(rtm)
 1.0 + 0.33333333333333326 t - 0.01111111111111103 t³ + 6.938893903907228e-18 t⁴ + 0.0003968253968253798 t⁵ - 1.2197274440461925e-17 t⁶ - 1.3227513227517302e-5 t⁷ - 6.979551485375435e-19 t⁸ + [851.817, 1741.59] t⁹

While the remainder of the last expression seems pretty large, note that its contribution is actually not so large

julia> remainder(frtm) * (domain(rtm))^9
[-1.29759e-05, 1.29759e-05]
lbenet commented 3 years ago

I'm not sure if the truncated exponential trick that you illustrate in this PR corresponds to a proper rigorous bound. Yet, from the RTaylorModel1 I obtained above, you can convert it into a TaylorModel1 essentially using the bound I obtained above.