# L-7 MCS 507 Wed 6 Sep 2023 : intvalsolved.py

"""
The example is taken from the paper of Stefano Taschini:
"Interval Arithmetic: Python Implementation and Applications"
Proceedings of the 7th Python in Science Conference (SciPy 2008).

We show how to obtain a correct numerical approximation for this example
using the interval arithmetic in the mpmath of SymPy.
The width of the interval gives us an upper bound for the error.
"""

import sympy as sp
X, Y = sp.var('x,y')
G = (sp.Rational(33375)/100 - X**2)*Y**6 \
  + X**2*(11*X**2*Y**2 - 121*Y**4 - 2) \
  + sp.Rational(55)/10*Y**8 \
  + sp.Rational(1)*X/(2*Y)
(A, B) = (77617, 33096)
print('evaluating', G, 'at', (A, B))
E = sp.Subs(G, (X, Y), (A, B)).doit()
E15 = E.evalf(15)
print('exact value :', E, '~', E15)
# now we use interval arithmetic
from mpmath import iv
print('using 35 decimal places ...')
iv.dps = 35
IV_F = lambda x, y: (iv.mpf('333.75') \
- x**2)*y**6 \
+ x**2*(iv.mpf('11')*x**2*y**2 \
- iv.mpf('121')*y**4 - iv.mpf('2')) \
+ iv.mpf('5.5')*y**8 + x/(iv.mpf('2')*y)
IV_A = iv.mpf(str(A))
IV_B = iv.mpf(str(B))
IV_Z = IV_F(IV_A, IV_B)
print(IV_Z)
print('using 36 decimal places ...')
iv.dps = 36
IV_F = lambda x, y: (iv.mpf('333.75') \
- x**2)*y**6 \
+ x**2*(iv.mpf('11')*x**2*y**2 \
- iv.mpf('121')*y**4 - iv.mpf('2')) \
+ iv.mpf('5.5')*y**8 + x/(iv.mpf('2')*y)
IV_A = iv.mpf(str(A))
IV_B = iv.mpf(str(B))
IV_Z = IV_F(IV_A, IV_B)
print(IV_Z)