{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "In lecture 25 of MCS 320, we consider symbolic-numeric computing" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# 1. Interval Arithmetic" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "What is the right precision to obtain accurate results?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To motivate the use of interval arithmetic, we take an example\n", "from the paper of Stefano Taschini in the proceedings of SciPy 2008.\n", "We convert the floating-point coefficients to rational numbers." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(x, y) |--> 11/2*y^8 - 1/4*(4*x^2 - 1335)*y^6 + (11*x^2*y^2 - 121*y^4 - 2)*x^2 + 1/2*x/y" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "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)\n", "F" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "-54767/66192" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "(A, B) = (77617, 33096)\n", "exact = F(A, B)\n", "exact" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "-0.8273960599" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "exact.n(digits=10)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[-6.662869420e25,\n", " 1.0742565174e25,\n", " -1.37766841984e23,\n", " 1.334889447551e22,\n", " 3.9041564432186e21,\n", " -8.18209863729137e20,\n", " -6.553155272242627e18,\n", " -6.5531536229751849e18,\n", " 3.64375473385373696e17,\n", " 4.087503195734016000e15,\n", " 4.0875015851212800000e15,\n", " 1.46852011835392000000e14,\n", " 6.114519285760000000000e12,\n", " 6.1145192857600000000000e12,\n", " 6.72051363840000000000000e10,\n", " -1.514340352000000000000000e9,\n", " -1.5143403520000000000000000e9,\n", " 9.62723840000000000000000000e7,\n", " -3.794534400000000000000000000e7,\n", " 3.9976940000000000000000000000e6,\n", " -196610.000000000000000000000000,\n", " -2.000000000000000000000000000000,\n", " -2.0000000000000000000000000000000,\n", " -2.00000000000000000000000000000000,\n", " -2.000000000000000000000000000000000,\n", " -2.0000000000000000000000000000000000,\n", " -1.00000000000000000000000000000000000,\n", " -0.8125000000000000000000000000000000000,\n", " -0.82812500000000000000000000000000000000,\n", " -0.827148437500000000000000000000000000000]" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "[F(A.n(digits=d), B.n(digits=d)) for d in range(10, 40)]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We see that we eventually get to an accurate result. However, observe the plateaus which could give the impression that, as the result no longer changes, we have reached an accurate result." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Real Interval Field with 53 bits of precision" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "RIF" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.?e22" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "y = F(RIF(A), RIF(B))\n", "y" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To see better the range of numbers, we ask to see the number in the ``brackets`` style." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "'[-3.1793933514545646e21 .. 3.9041567246936242e21]'" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "RIF(y).str(style='brackets')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Observe the magnitude of the numbers. The range of the interval is huge!" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Real Interval Field with 100 bits of precision" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "R100 = RealIntervalField(prec=100)\n", "R100" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "'[-4.3909140000000000000000000000000e6 .. 1.2386302000000000000000000000000e7]'" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "y100 = F(R100(A), R100(B))\n", "R100(y100).str(style='brackets')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "After doubling the precision to 100 bits, observe that the magnitude dropped from $10^{21}$ down to $10^7$." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "'[-0.82739605994682136814116587127405770020482922433302519493736327 .. -0.82739605994682136814115925382915727598343025306348863523453474]'" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "R200 = RealIntervalField(prec=200)\n", "y200 = F(R200(A), R200(B))\n", "R200(y200).str(style='brackets')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Working with 200 bits gives an accurate result. The side of the interval gives a bound on the error." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "-0.82739605994682136814116?" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "y200" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "6.6174449004242213989712695365597028285264968872070312500000e-24" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ry200 = R200(y200)\n", "ry200.upper() - ry200.lower()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The error is defined as the difference between the upper and lower limit of the interval and we see that we have about 24 correct decimal places." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# 2. Symbolic-Numeric Factorization" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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 with algebraic number till the polynomial factors as a product of linear forms." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A symbolic-numeric factorization works over the ``QQbar`` number field, which uses the complex interval field." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Univariate Polynomial Ring in x over Algebraic Field" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "P. = PolynomialRing(QQbar)\n", "P" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(x - 0.3411639019140097? - 1.161541399997252?*I) * (x - 0.3411639019140097? + 1.161541399997252?*I) * (x + 0.6823278038280193?)" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "q = x^3 + x + 1\n", "f = factor(q)\n", "f" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[(-0.6823278038280193?, 1),\n", " (0.3411639019140097? - 1.161541399997252?*I, 1),\n", " (0.3411639019140097? + 1.161541399997252?*I, 1)]" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "q.roots(QQbar)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If we want a more accurate symbolic-numeric factorization, then we can compute the roots at a higher precision, for example with 200 bits." ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Complex Interval Field with 200 bits of precision" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "C = ComplexIntervalField(prec=200)\n", "C" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[(-0.682327803828019327369483739711048256891188581897998577803728?, 1),\n", " (0.34116390191400966368474186985552412844559429094899928890187? - 1.16154139999725193608791768724717407484314725802151429063617?*I,\n", " 1),\n", " (0.34116390191400966368474186985552412844559429094899928890187? + 1.16154139999725193608791768724717407484314725802151429063617?*I,\n", " 1)]" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "q.roots(C)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# 3. Constrained Optimization" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To end the ``calculus`` chapter and to transition into the solving chapter, consider a problem of constrained optimization." ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "optimize 2*x*y*z + x^2 - z^2 subject to x^2 + y^2 + z^2 - 1 == 0\n" ] } ], "source": [ "x, y, z = var('x, y, z')\n", "target = x^2 + 2*x*y*z - z^2\n", "constraint = x^2 + y^2 + z^2 - 1 == 0\n", "print('optimize', target, 'subject to', constraint)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us make plot, for the value of the target equal to two." ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n" ], "text/plain": [ "Graphics3d Object" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "tplot = implicit_plot3d(target-2, (x,-2,2), (y,-2,2), (z,-2,2), color='red')\n", "cplot = implicit_plot3d(constraint.lhs(), (x,-2,2), (y,-2,2), (z,-2,2))\n", "show(tplot+cplot)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If the value of the target equals one, then we see that the red sheets touch the blue ball." ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n" ], "text/plain": [ "Graphics3d Object" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "tplot = implicit_plot3d(target-1, (x,-2,2), (y,-2,2), (z,-2,2), color='red')\n", "cplot = implicit_plot3d(constraint.lhs(), (x,-2,2), (y,-2,2), (z,-2,2))\n", "show(tplot+cplot)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To apply the method of Lagrange multipliers we need the gradients." ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "target as function : (x, y, z) |--> 2*x*y*z + x^2 - z^2\n", "constraint as function : (x, y, z) |--> x^2 + y^2 + z^2 - 1\n" ] } ], "source": [ "ft(x,y,z) = target\n", "fc(x,y,z) = constraint.lhs()\n", "print('target as function :', ft)\n", "print('constraint as function :', fc)" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(x, y, z) |--> (2*y*z + 2*x, 2*x*z, 2*x*y - 2*z)" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "gt = ft.diff()\n", "gt" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(x, y, z) |--> (2*x, 2*y, 2*z)" ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "gc = fc.diff()\n", "gc" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we set up the system. Because ``lambda`` is reserved in Python, using ``lambda`` as the name of a variable is wrong." ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[-2*w*x + 2*y*z + 2*x == 0,\n", " -2*w*y + 2*x*z == 0,\n", " 2*x*y - 2*w*z - 2*z == 0,\n", " x^2 + y^2 + z^2 - 1 == 0]" ] }, "execution_count": 24, "metadata": {}, "output_type": "execute_result" } ], "source": [ "w = var('w')\n", "eqs = [gt(x,y,z)[k] - w*gc(x,y,z)[k] == 0 for k in range(3)]\n", "eqs.append(constraint)\n", "eqs" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Observe that the constraint was appended to the list of equations. We then just solve." ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[[x == 1, y == 0, z == 0, w == 1], [x == -1, y == 0, z == 0, w == 1], [x == 0, y == 0, z == -1, w == -1], [x == 0, y == 0, z == 1, w == -1], [x == 0, y == -1, z == 0, w == 0], [x == 0, y == 1, z == 0, w == 0]]" ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sols = solve(eqs,[x,y,z,w])\n", "sols" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To evaluate the solutions at the target, we must convert the coordinates of the solutions into tuples and drop the value for w." ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[[x == 1, y == 0, z == 0, w == 1],\n", " [x == -1, y == 0, z == 0, w == 1],\n", " [x == 0, y == 0, z == -1, w == -1],\n", " [x == 0, y == 0, z == 1, w == -1],\n", " [x == 0, y == -1, z == 0, w == 0],\n", " [x == 0, y == 1, z == 0, w == 0]]" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ "[sol for sol in sols]" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[[1, 0, 0, 1],\n", " [-1, 0, 0, 1],\n", " [0, 0, -1, -1],\n", " [0, 0, 1, -1],\n", " [0, -1, 0, 0],\n", " [0, 1, 0, 0]]" ] }, "execution_count": 27, "metadata": {}, "output_type": "execute_result" } ], "source": [ "L = [[e.rhs() for e in sol] for sol in sols]\n", "L" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[(1, 0, 0), (-1, 0, 0), (0, 0, -1), (0, 0, 1), (0, -1, 0), (0, 1, 0)]" ] }, "execution_count": 28, "metadata": {}, "output_type": "execute_result" } ], "source": [ "T = [tuple(e[0:3]) for e in L]\n", "T" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[1, 1, -1, -1, 0, 0]" ] }, "execution_count": 29, "metadata": {}, "output_type": "execute_result" } ], "source": [ "[ft(x,y,z) for (x,y,z) in T]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This problem has two maxima, valued at 1, and two minima, with value -1." ] } ], "metadata": { "kernelspec": { "display_name": "SageMath 10.3", "language": "sage", "name": "sagemath" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10" } }, "nbformat": 4, "nbformat_minor": 4 }