Loading [MathJax]/jax/output/HTML-CSS/jax.js

First time here? Check out the FAQ!

Ask Your Question
1

Checking accuracy

asked 2 years ago

jarkky gravatar image

updated 2 years ago

Max Alekseyev gravatar image

.

sage: a=-1/2
sage: b=sqrt(3)/2
sage: s=(3/8*(a^3+sqrt(-a^4*b^2-2*a^2*b^4-b^6)+a*b^2)^(1/3)-(-5184*a^2-5184*b^2)/(13824*(a^3+sqrt(-a^4*b^2-2*a^2*b^4-b^6+a*b^2)^(1/3))+a/4)^(1/3)
sage: N(s,digits=200)       
    0.76604444311897803520239265055541667393583245708039524585404528464215538885687472352822927668054849344996248839987004553419542483163334907118418648890238602810825947820108992619243612738937829040787703
sage: N(sin(50/180*pi),digits=200)    
    0.76604444311897803520239265055541667393583245708039524585404528464215538885687472352822927668054849344996248839987004553419542483163334907118418648890238602810825947820108992619243612738937829040787703
sage: N(sin(50/180*pi)-s,digits=200)  
    -8.1651261039391273250903677461338074230805364330263310047708623423822773669493747438272352951232512935097984755386384444942102172477771762845811026642442173532418045576803620147104496959316258647247981e-202

The last calculation shows strange result. The difference should be exactly zero... (???) The text formatting does not look working well here.

Preview: (hide)

Comments

your expression for s can't be evaluated... Please check !

Emmanuel Charpentier gravatar imageEmmanuel Charpentier ( 2 years ago )

2 Answers

Sort by » oldest newest most voted
1

answered 2 years ago

Emmanuel Charpentier gravatar image

Furthering tmonteil's answer

N will try to refine its answer until the `digits' digits no longer change after refinement, and prints the last refined value. The fact that it stops and returns an absolute value of the order of 1O202 means that this value has 200 exact digits, which are all 0. You cane therefore state confidently that this absolute value is inferior to 10200, therefore not distingushable of 0 with 200 digits precision.

Computing an exact value

You seem to aim at computing sin50π180=sin5π18. Given sinx and cosx, `Sage' can easily compute sin5x :

sage: sin(5*x)==sin(5*x).trig_expand()
sin(5*x) == 5*cos(x)^4*sin(x) - 10*cos(x)^2*sin(x)^3 + sin(x)^5

Your problem is now therefore reduced to computing sinπ18 and cosπ18. Remarking that sinπ2=1 cosπ2=0, we can compute them by noting that these quantities are roots of the polynomial equations :

Sys=[u(9*x).trig_expand()==u(pi/2) for u in (sin, cos)]
Sys
[9*cos(x)^8*sin(x) - 84*cos(x)^6*sin(x)^3 + 126*cos(x)^4*sin(x)^5 - 36*cos(x)^2*sin(x)^7 + sin(x)^9 == 1,
  cos(x)^9 - 36*cos(x)^7*sin(x)^2 + 126*cos(x)^5*sin(x)^4 - 84*cos(x)^3*sin(x)^6 + 9*cos(x)*sin(x)^8 == 0]

Here, we could try to solve this polynomial system in SR by typing solve([u.subs({sin(x):S, cos(x):C}) for u in Sys],[S,C]) ; this would give us a whole lot of solutions in C, expressed using floating point numbers (i. e. inexact solutions).

A better solution is to get exact roots by moving to the ring of polynomials over the (real) algebraic number field. Sage is able to compute "exact" roots of such polynomials :

R1.<s,c>=AA[]
PSys=[R1((u.lhs()-u.rhs()).subs({sin(x):s, cos(x):c})) for u in Sys]
PSys

[s^9 - 36*s^7*c^2 + 126*s^5*c^4 - 84*s^3*c^6 + 9*s*c^8 - 1,
  9*s^8*c - 84*s^6*c^3 + 126*s^4*c^5 - 36*s^2*c^7 + c^9]

The roots of this polynomial system are given by :

PSol=R1.ideal(PSys).variety()

Note that there are 9 such real solutions. To select the right one, we can use a tad of elementary geometric intuition and :

  • restrict ourselves to the first quadrant (i. e. s' andc' both positive) ;

  • note that in this quadrant, pi/18 has the smallest sine and the largest cosine.

It turns out that both have a radical expression :

Q1=[u for u in PSol if u[s]>0 and u[c]>0]
S1=min([u[s] for u in Q1])
C1=max([u[c] for u in Q1])
(S1.radical_expression(), C1.radical_expression())

(-1/2*(I*sqrt(3) + 1)*(1/16*I*sqrt(3) - 1/16)^(1/3) - 1/8*(-I*sqrt(3) + 1)/(1/16*I*sqrt(3) - 1/16)^(1/3),
 1/4*(4*(1/128*I*sqrt(3) + 1/128)^(1/3) + 1)/(1/128*I*sqrt(3) + 1/128)^(1/6))

Computing sin5π18 is now trivial :

S50=sin(5*x).trig_expand().subs({sin(x):S1, cos(x):C1})
S50

0.7660444431189781?

Note that this is not a mere numerical approximation, but a representation of an algorithm allowing to compute numerical approximations to arbitrary precision, thus representing a unique algebraic real number.

This result can be exactly verified:

bool(S50==sin(50*pi/180))

True

HTH,

Preview: (hide)
link

Comments

Is there any method to get rid of the "irreduclible cuberoots"?

jarkky gravatar imagejarkky ( 2 years ago )
1

answered 2 years ago

tmonteil gravatar image

I can not reproduce your computation (the definition of s seems uncomplete), but you should notice that -8e-202 is extremly small (in particular, the first 200 digits of s and N(sin(50/180*pi)-s,digits=200) are equal).

If you want a result with an explicit bound on the error, you should probably use RealBallField or RealIntervalField.

Preview: (hide)
link

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: 2 years ago

Seen: 282 times

Last updated: Apr 12 '22

Related questions