Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add example evaluating rational function #99

Open
wants to merge 2 commits into
base: master
Choose a base branch
from

Conversation

mforets
Copy link
Contributor

@mforets mforets commented Jan 19, 2021

This documentation PR is a continuation of #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

@mforets mforets mentioned this pull request Jan 20, 2021
@mforets mforets mentioned this pull request Feb 15, 2021
@lbenet
Copy link
Member

lbenet commented Feb 19, 2021

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- 0.6666666666666667+ 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- 1.0+ 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+ 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
Copy link
Member

lbenet commented Feb 19, 2021

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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

2 participants