Lecture 24: Symbolic-Numeric Computation¶
With interval arithmetic we compute floating-point approximations and estimates for the errors. The manipulation of symbolic expressions is used to formulate equations to solve, which we illustrate via the application of the Lagrange multiplier technique to a constrained optimization problem.
There is one important type of arithmetic we have not covered yet, and that is interval arithmetic.
R = RealIntervalField() R
R is a
Real Interval Field with 53 bits of precision.
By default, the
without argument (also abbreviated as RIF)
works with double precision floating-point arithmetic.
We can coerce or convert into this field.
api = R(pi) print(api, 'has type', type(api)) print(api.str(style='brackets'))
and we respectively see
3.141592653589794? (observe the question mark at the end)
and the interval
[3.1415926535897931 .. 3.1415926535897936]
Interval arithmetic provides operations to add, substract, multiply, and divide intervals. The width of the interval gives an upper bound on the roundoff error. With interval arithmetic we compute at the same time a floating-point approximation and an estimate for the error on this approximation.
To motivate the use of interval arithmetic, we take an example from the paper of Stefano Taschini in the proceedings of SciPy 2008. The problem is to evaluate
at \((77617, 33096)\).
To calculate the exact value, we convert the floating-point coefficients to rational numbers.
F(x,y) = (QQ(333.75) - x**2)*y**6 + x**2*(11*x**2*y**2 - 121*y**4 - 2) + QQ(5.5)*y**8 + x/(2*y) print F (A, B) = (77617, 33096) exact = F(A, B) print('the exact value :', exact, RDF(exact))
and the exact value is
or approximated in double precision
With increasing multiprecision, we may gain more confidence in the results.
f53 = F(A.n(53), B.n(53)) print(f53) for k in range(54, 105, 10): print(F(A.n(k), B.n(k)))
and we see
3.62381792172646e20 -8.18209863729137e20 3.64375473385373696e17 1.46852011835392000000e14 6.72051363840000000000000e10 9.62723840000000000000000000e7 -196610.0000000000000000000000000
Let us continue.
for k in range(100, 130, 5): print(F(A.n(k), B.n(k)))
and then we get the following:
3.9976940000000000000000000000e6 65534.00000000000000000000000000 -2.0000000000000000000000000000000 -2.000000000000000000000000000000000 -2.0000000000000000000000000000000000 -0.7500000000000000000000000000000000000
The sequence of repeated
-2.0 values is especially troublesome
as one would reach the wrong conclusion that the precision is enough
once the sequence has stabilized.
This trial and error process is not really encouraging…
R = RealIntervalField(130) z = F(R(A), R(B)) print(R(z).str(style='brackets'))
and we get
[-0.82812500000000000000000000000000000000000 .. -0.82031250000000000000000000000000000000000]
In a numeric factorization of a polynomial in one variable,
we compute the complex roots of the polynomial.
In a symbolic factorization, we extend the field of rational
numbers till the polynomial can be written as a product of linear factors.
A symbolic-numeric factorization uses the complex interval field,
that is: the roots of the polynomial are computed with complex interval
arithmetic in a specific precision. In SageMath, a symbolic-numeric
factorization is computed using the
QQbar set of numbers.
P.<x> = PolynomialRing(QQbar)
The above statement defines
P as the univariate polynomial ring
x over the algebraic field. Any expression in
x will be
interpreted as belonging to
P, as is the case with
q = x^3 + x + 1 factor(q)
We see the three factors of
q with the roots written in the
interval notation and the
I imaginary unit.
To see the roots of
QQbar we can do
and we same the same intervals as in the output of
If we want a more accurate symbolic-numeric factorization, for example with 200 bits of precision, then we first define a new complex interval field before computing the roots in that new defined field, as done by the code below.
C = ComplexIntervalField(prec=200) q.roots(C)
Grabbing the output of
q.roots(C) defines then the factors
in a symbolic-numeric factorization. Computing the product shows
the backward error of this computation, that is: the polynomial
q for which
q.roots(C) are the exact roots.
The method of Lagrange multipliers allows to solve a constrained optimization problem.
var('x,y,z') target = x^2 + 2*x*y*z - z^2 constraint = x^2 + y^2 + z^2 - 1 == 0 print('optimize', target, 'subject to', constraint)
and we get the printed message
optimize 2*x*y*z + x^2 - z^2 subject to x^2 + y^2 + z^2 - 1 == 0.
To apply the method of Lagrange multipliers, we need to compute the gradients.
For this, we convert the expressions into functions.
Observe that the constraint was given as an equation,
so we must select the left hand side of the equation
when we convert the expression to a function.
If we want that the target evaluates at 2, then we see there is no solution, see Fig. 30 made with the commands below.
tplot = implicit_plot3d(target-2, (x,-2,2), (y,-2,2), (z,-2,2), color='red') cplot = implicit_plot3d(constraint.lhs(), (x,-2,2), (y,-2,2), (z,-2,2)) (tplot+cplot).show()
If we want solutions, then the contraint and target must touch, which is the case when the value of the target equals one. Let us make the plot, see Fig. 31.
tplot = implicit_plot3d(target-1, (x,-2,2), (y,-2,2), (z,-2,2), color='red') cplot = implicit_plot3d(constraint.lhs(), (x,-2,2), (y,-2,2), (z,-2,2)) (tplot+cplot).show()
ft(x,y,z) = target fc(x,y,z) = constraint.lhs() print('target as function :', ft) print('constraint as function :', fc) gt = ft.diff() gc = fc.diff() print('gradient of target :', gt) print('gradient of constraint :', gc)
and the results are
target as function : (x, y, z) |--> 2*x*y*z + x^2 - z^2 constraint as function : (x, y, z) |--> x^2 + y^2 + z^2 - 1 gradient of target : (x, y, z) |--> (2*y*z + 2*x, 2*x*z, 2*x*y - 2*z) gradient of constraint : (x, y, z) |--> (2*x, 2*y, 2*z)
Now we set up the system.
w = var('w') eqs = [gt(x,y,z)[k] - w*gc(x,y,z)[k] == 0 for k in range(3)] eqs.append(constraint) eqs
The polynomial system we obtain is thus
[-2*w*x + 2*y*z + 2*x == 0, -2*w*y + 2*x*z == 0, 2*x*y - 2*w*z - 2*z == 0, x^2 + y^2 + z^2 - 1 == 0]
and then we just do solve…
sols = solve(eqs,[x,y,z,w]) print(sols) print('found', len(sols), 'solutions')
to see that we found the candidate optimal solutions.
x = RIF(-1, +1); print 1/x. Interpret the result.
Compare the evaluation of \(p = x^2 - 4x\) and \(q = x (x-4)\) at the interval [1, 4]. Which form of the expression, \(p\) or \(q\) leads to the most accurate result? Justify your answer.
Find those points on the surface \(x^2 - x y + y^2 - z^2 = 1\) closest to the origin. Set up the system with Lagrange multipliers and solve.
Consider the surface defined by \(g(x,y) = x^2 + 3 x y + 2 y^2 + 7 x = 0\). Find those points \((x,y)\) that satisfy \(g(x,y) = 0\) and that are closest to the origin.
Define this problem as a polynomial system.
Bring this system in triangular form. How many solutions can this system have? Justify your answer.
Find all real solutions.
Consider the function \(f(x) = 1 - 10 x + 0.01 e^x\). Use interval arithmetic to show that \(f(x) = 0\) has a solution in the interval \([9, 10]\).
Compute a symbolic-numeric factorization of \(p = x^4 + 3 x^3 + 2 x^2 + x + 9\) with 20 bits of precision. Use the brackets style to display the roots. What is backward error?