Ask Your Question
1

Undeterministic numerical approximation

asked 2017-06-15 04:29:17 +0100

LaurentClaessens gravatar image

Let the following expression of x :

f=1/152444172305856930250752000000*x^28 + 1/10888869450418352160768000000*x^27 + 1/15511210043330985984000000*x^25 + 1/310224200866619719680000*x^24 + 1/25852016738884976640000*x^23 + 1/51090942171709440000*x^21 + 1/1216451004088320000*x^20 + 1/121645100408832000*x^19 + 1/355687428096000*x^17 + 1/10461394944000*x^16 + 1/1307674368000*x^15 + 1/6227020800*x^13 + 1/239500800*x^12 + 1/39916800*x^11 + 1/362880*x^9 + 1/20160*x^8 + 1/5040*x^7 + 1/120*x^5 + 1/12*x^4 + 1/6*x^3 + x - cos(x) + 2 -exp(x)

I noticed that evaluating that for the value x=numerical_approximation(10) leads to non-deterministic results (from the 12th decimal).

I guess that this is a "cancellation error" : the last term (-exp(x)) and the whole stuff before are both of the order of magnitude 22026.

I'm only interested in the first 3 or 4 decimals.

  • QUESTION 1 How can I "force" sage to provide deterministic results ? Doing numerical_approx(...,prec=30) makes the job, but the "30" is quite heuristic; I'm not sure that this is a good approach.

  • QUESTION 2 (not really sage-related) How should I unit-test a function that has undeterministic bahaviour like that ?

edit retag flag offensive close merge delete

1 Answer

Sort by ยป oldest newest most voted
0

answered 2017-06-15 21:28:25 +0100

dan_fulea gravatar image

updated 2017-06-15 21:31:15 +0100

The following is a possibility to make the result predictible:

def f(x):
    return \
        ( 1/152444172305856930250752000000*x^28
          + 1/10888869450418352160768000000*x^27
          + 1/15511210043330985984000000*x^25
          + 1/310224200866619719680000*x^24
          + 1/25852016738884976640000*x^23
          + 1/51090942171709440000*x^21
          + 1/1216451004088320000*x^20
          + 1/121645100408832000*x^19
          + 1/355687428096000*x^17
          + 1/10461394944000*x^16
          + 1/1307674368000*x^15
          + 1/6227020800*x^13
          + 1/239500800*x^12
          + 1/39916800*x^11
          + 1/362880*x^9
          + 1/20160*x^8
          + 1/5040*x^7
          + 1/120*x^5
          + 1/12*x^4
          + 1/6*x^3
          + x
          - cos(x) + 2 -exp(x) )

R = RealField( 200 )
C = ComplexField( 200 )

f( C(10) )

f( R(10) )
f( R(10) )
f( R(10) )

Now the last four calls of f( C(10) ) and f( R(10) ) lead to...

-0.013417186790728922927885279844822282173830250128225204250732
-0.013417186790728922927885279844822282173830250128225204250732
-0.013417186790728922927885279844822282173830250128225204250732
-0.013417186790728922927885279844822282173830250128225204250732

We can also compare with...

sage: f(x) = exp(x) + cos(x) - 2
sage: p = taylor( f(x), x, 0, 28 )
sage: g(x) = sum( [ c*x^pow for c,pow in p.coefficients() ] ) - f(x) 
sage: g(10)
-cos(10) - e^10 + 185630709328907872984/8427947351530707
sage: g(10.)
-0.0134171867909140
sage: g( RealField(100)(10) )
-0.013417186790728922927885322847
sage: g( RealField(200)(10) )
-0.013417186790728922927885279844822282173830250128225204254272
sage: g( RealField(400)(10) )
-0.0134171867907289229278852798448222821738302501282252042588478039342715326239015599771064232777651085304454608204788898295
sage: gp( g(10) )
-0.01341718679072892292788521428
sage: g(10).n()
-0.0134171867909140
sage: g(10).n( digits=100 )
-0.01341718679072892292788527984482228217383025012822520425884780393427153262390155997710642327776499594
edit flag offensive delete link more

Comments

Thanks for the answer. Since I am only interested in the firsts decimals (these are coordinates of of points on a A4 paper, in centimeters), I didn't though to increase the precision.

I'l have to tests that.

LaurentClaessens gravatar imageLaurentClaessens ( 2017-06-16 16:59:03 +0100 )edit

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower

Stats

Asked: 2017-06-15 04:29:17 +0100

Seen: 746 times

Last updated: Jun 15 '17