============================
Symbolic vs numerics in Sage
============================
------------------
Parent and element
------------------
When using Sage, it is important to understand where your objects live. Or more
mathematically in which structure they belong to. The number `1` as an integer
will behave differently as the number `1` as a rational number...
In Sage this *ambient space* is called parent::
sage: a = 4
sage: parent(a)
sage: a = 2/3
sage: parent(a)
sage: a = 2.5
sage: parent(a)
sage: a = pi
sage: parent(pi)
The parent determines the behavior of many operations::
sage: x = SR.var('x')
sage: p = 12*x^2 + 12*x + 24
sage: print parent(p)
sage: print p.factor()
sage: x = polygen(ZZ)
sage: p = 12*x^2 + 12*x + 24
sage: print parent(p)
sage: print p.factor()
sage: x = polygen(QQ)
sage: p = 12*x^2 + 12*x + 24
sage: print parent(p)
sage: print p.factor()
sage: x = polygen(QQbar)
sage: p = 12*x^2 + 12*x + 24
sage: print parent(p)
sage: print p.factor()
Even for real numbers there is a huge variety of them
- exact numbers (integer, rational, algebraic numbers)
- floating point numbers (real or complex)
- intervals and balls (real or complex)
- symbolic (the "big garbage")
Floating point numbers are **inexact** by nature... it is particularly
dangerous in an **unstable** situation. They should be used with a lot of care.
-----------------------------
Exercise 1: stability matters
-----------------------------
Look at the following four cells and try to guess what will be the output
*before* executing them.::
sage: x0 = 1/3
sage: print 4*x0 - 1 == x0
sage: x = x0
sage: for i in range(100):
....: x = 4*x - 1
sage: print x
sage: x0 = RR(1/3)
sage: print 4*x0 - 1 == x0
sage: x = x0
sage: for i in range(100):
....: x = 4*x - 1
sage: print x
sage: #edit
--------------------------
Exercise 2: parent matters
--------------------------
Compare the following two cells. Which computations is the most accurate? the fastest?::
sage: a = sqrt(2) + sqrt(3) + sqrt(5) + sqrt(7) - 6
sage: for i in range(8):
....: a = a^4 - 4*a^3 + 4*a^2
sage: print a.numerical_approx()
sage: a = sqrt(2.) + sqrt(3.) + sqrt(5.) + sqrt(7.) - 6
sage: for i in range(8):
....: a = a^4 - 4*a^3 + 4*a^2
sage: print a.numerical_approx()
sage: #edit
In both cases, add some code to see what is the parent of ``a``.
--------
Coercion
--------
Coercion is what happens when you mix numbers (or more generally objects) of different kinds.
For example the addition of an integer and a rational will be a rational::
sage: parent(5 + 2/3)
But what about:
- the addition of a rational with a polynomial with integer coefficients?
- division of two integers?
- division of two polynomials?
::
sage: # edit
sage: # edit
sage: # edit
----------------------------
Interval and ball arithmetic
----------------------------
Interval and balls are two kinds of floating point arithmetic that take care of
error propagation. These types correspond to pairs of floating point numbers
that represent respectively:
- the left and right endpoints of an interval
- the center and the radius of a ball
The parents are ``RIF`` (for *real interval field*) and ``RBF`` (for *real ball field*)::
sage: RIF(sqrt(2))
sage: RBF(sqrt(2))
The arithmetic is done in such way that for any operation `f(X)` applied to an
interval or a ball `X` *contains* all of `\{f(x): x \in X\}`::
sage: a = RIF(1/2,1) # the interval [1/2,1]
sage: print a.endpoints()
sage: b = a^2
sage: print b.endpoints()
sage: a = RBF(RIF(1/2,1)) # the same seen as a ball
sage: print a.endpoints()
sage: b = a^2
sage: print b.endpoints()
You can prove theorems using balls or intervals (which is much harder using
floating point)!
You can apply most of standard functions on intervals and balls::
sage: a = RIF(1/3)
sage: print a.cos()
sage: print a.tan()
sage: # etc
**Exercise:**
Find an example of a singleton interval `X = [x,x]` and a function `f` where the
`f(X)` is not a singleton::
sage: # edit
sage: # edit
sage: # edit
---------
Precision
---------
Most floating point parents admit a precision attribute::
sage: R = RealField(64)
sage: R(1/3)
sage: R = RealField(256)
sage: R(1/3)
sage: R = RealIntervalField(64)
sage: R(1/3)
sage: R = RealIntervalField(256)
sage: R(1/3)
------------------------------
Playing with algebraic numbers
------------------------------
Since sage-7.1 it is possible to work with embedded number fields in the set of
real numbers::
sage: R. = PolynomialRing(ZZ)
sage: K. = NumberField(x^3 - 2, embedding=AA(2)**(1/3))
sage: 1 < cbrt2 < 2
sage: continued_fraction(cbrt2)
These numbers are exact by nature, but comparing them can lead to expensive
computations. Moreover the complexity of addition/multiplication increase
(linearly) with the degree of the number field.
---------------------------------
Integrating Abelian differentials
---------------------------------
In this section, you are invited to play with Abelian differential in the
complex plane. The aim is to:
- compute holonomy of paths $$ \int_\gamma f(x) dx = \int_0^1 f(\gamma(t)) \gamma'(t) dt.$$
- integrate $f(x) dx$ in other words find numerica solutions to the one
parameter family of ODEs $$f \left(\gamma(t) \right)\ \gamma'(t) = e^{i \theta}.$$
Two basic functions in Sage to perform this task are:
- `numerical_integral`
- `ode_solver`
::
sage: class DifferentialCaller:
....: def __init__(self, f):
....: vars = f.variables()
....: assert len(vars) == 1
....: var = vars[0]
....:
....: x = SR.var('x')
....: y = SR.var('y')
....: theta = SR.var('theta')
....: I = SR('I')
....: g = exp(I * theta) / f(x + I*y)
....: g1 = g.real()
....: self.g1 = fast_callable(g1, vars=[x,y,theta], domain=float)
....: g2 = g.imag()
....: self.g2 = fast_callable(g2, vars=[x,y,theta], domain=float)
....: dF = exp(I * theta) * f.derivative(var)(x + I*y) / (f(x+I*y) * f(x+I*y))
....: dF = dF.full_simplify()
....: print dF
....: self.dg1x = fast_callable(dF.real(), vars=[x,y,theta], domain=float)
....: self.dg1y = fast_callable((I * dF).real(), vars=[x,y,theta], domain=float)
....: self.dg2x = fast_callable(dF.imag(), vars=[x,y,theta], domain=float)
....: self.dg2y = fast_callable((I * dF).imag(), vars=[x,y,theta], domain=float)
....: def function(self, t, z, params):
....: return (self.g1(z[0], z[1], params[0]), self.g2(z[0], z[1], params[0]))
....:
....: def jacobian(self, t, z, params):
....: return ((self.dg1x(z[0], z[1], params[0]), self.dg2x(z[0],z[1],params[0])),
....: (self.dg1y(z[0], z[1], params[0]), self.dg2y(z[0],z[1],params[0])))
sage: f(z) = exp(1/z)
sage: D = DifferentialCaller(f)
-e^(I*theta - 1/(x + I*y))/(x^2 + 2*I*x*y - y^2)
sage: T = ode_solver()
sage: T.function = D.function
sage: #T.jacobian = D.jacobian
sage: G = Graphics()
sage: Z0 = [(-1,0.3)]
sage: theta_step = 0.5
sage: for z0 in Z0:
....: theta = 0.0
....: while theta < 2*RR.pi()+theta_step/2:
....: T.ode_solve(y_0=z0, t_span=[0,1], params=[theta], num_points=100)
....: G += line2d([z for t,z in T.solution], color='blue')
....: theta += theta_step
sage: G.show()