Checking accuracy

.

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.

edit retag close merge delete

( 2022-04-11 09:22:12 +0200 )edit

Sort by » oldest newest most voted

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 $1O^{-202}$ means that this value has $200$ exact digits, which are all 0. You cane therefore state confidently that this absolute value is inferior to $10^{-200}$, therefore not distingushable of $0$ with 200 digits precision.

Computing an exact value

You seem to aim at computing $\sin\frac{50\pi}{180}=\sin\frac{5\pi}{18}$. Given $\sin{}x$ and $\cos{}x$, Sage' can easily compute $\sin{5x}$ :

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\frac{\pi}{18}$ and $\cos\frac{\pi}{18}$. Remarking that $\sin\frac{\pi}{2}=1\ \cos\frac{\pi}{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 $\mathbb{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])

(-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 $\sin\frac{5\pi}{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,

more

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

( 2022-05-07 13:45:31 +0200 )edit

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.

more