{ "cells": [ { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "## Math 157: Intro to Mathematical Software\n", "## UC San Diego, winter 2018\n" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "## January 24, 2018: SageMath (part 2 of 3)" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "Administrivia:\n", "\n", "- Homework 2 is due Friday, January 26 at 8pm. If you do not have it, contact course staff immediately.\n", "- We are continuing to monitor the course waitlist. If you are on it and still want to join the course, please continue to keep up with lectures and homework, and watch your email for further instructions." ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "# Finding Roots: symbolic\n", "You can use the solve command to solve for zeroes of a function." ] }, { "cell_type": "code", "execution_count": 67, "metadata": { "collapsed": false }, "outputs": [ { "ename": "SyntaxError", "evalue": "can't assign to operator (, line 1)", "output_type": "error", "traceback": [ "\u001b[0;36m File \u001b[0;32m\"\"\u001b[0;36m, line \u001b[0;32m1\u001b[0m\n\u001b[0;31m x**Integer(2)+Integer(3) = Integer(5) # This is an error; a Python assignment cannot have a complex expression on its left-hand side.\u001b[0m\n\u001b[0;31mSyntaxError\u001b[0m\u001b[0;31m:\u001b[0m can't assign to operator\n" ] } ], "source": [ "x^2+3 = 5 # This is an error; a Python assignment cannot have a complex expression on its left-hand side." ] }, { "cell_type": "code", "execution_count": 0, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "x^2 + 3 == 5" ] }, "execution_count": 1, "metadata": { }, "output_type": "execute_result" } ], "source": [ "x^2 + 3 == 5 ## This doesn't make sense in Python, but in Sage it is a valid object!" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "" ] }, "execution_count": 2, "metadata": { }, "output_type": "execute_result" } ], "source": [ "eqn = x^2 + 3 == 5\n", "show(eqn)" ] }, { "cell_type": "code", "execution_count": 48, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "f(x) = sin(x)+exp(x)" ] }, { "cell_type": "code", "execution_count": 49, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "" ] }, "execution_count": 49, "metadata": { }, "output_type": "execute_result" } ], "source": [ "show(f)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "x^2 == 2" ] }, "execution_count": 3, "metadata": { }, "output_type": "execute_result" } ], "source": [ "eqn.add_to_both_sides(-3)" ] }, { "cell_type": "code", "execution_count": 52, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "685f5ac125f8ac38119bdafc97232e090779b78b" } } ], "source": [ "v = [(0,0), (1,1), (1,2), (2,0), (3,1)]\n", "line(v, thickness=3, color='blue', zorder=-1) + points(v, pointsize=200, color='grey', zorder=1)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "pi == pi" ] }, "execution_count": 4, "metadata": { }, "output_type": "execute_result" } ], "source": [ "pi == pi" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "3.14159265358979" ] }, "execution_count": 6, "metadata": { }, "output_type": "execute_result" } ], "source": [ "N(pi)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 7, "metadata": { }, "output_type": "execute_result" } ], "source": [ "bool(pi == pi)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "False" ] }, "execution_count": 7, "metadata": { }, "output_type": "execute_result" } ], "source": [ "bool(eqn) # Should return False, because this is not an identity." ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "[x == -1/2*sqrt(((2*I*sqrt(15) + 2)^(2/3) + 4)/(2*I*sqrt(15) + 2)^(1/3)) - 1/2*sqrt(-(2*I*sqrt(15) + 2)^(1/3) - 4/(2*I*sqrt(15) + 2)^(1/3) + 4/sqrt(((2*I*sqrt(15) + 2)^(2/3) + 4)/(2*I*sqrt(15) + 2)^(1/3))), x == -1/2*sqrt(((2*I*sqrt(15) + 2)^(2/3) + 4)/(2*I*sqrt(15) + 2)^(1/3)) + 1/2*sqrt(-(2*I*sqrt(15) + 2)^(1/3) - 4/(2*I*sqrt(15) + 2)^(1/3) + 4/sqrt(((2*I*sqrt(15) + 2)^(2/3) + 4)/(2*I*sqrt(15) + 2)^(1/3))), x == 1/2*sqrt(((2*I*sqrt(15) + 2)^(2/3) + 4)/(2*I*sqrt(15) + 2)^(1/3)) - 1/2*sqrt(-(2*I*sqrt(15) + 2)^(1/3) - 4/(2*I*sqrt(15) + 2)^(1/3) - 4/sqrt(((2*I*sqrt(15) + 2)^(2/3) + 4)/(2*I*sqrt(15) + 2)^(1/3))), x == 1/2*sqrt(((2*I*sqrt(15) + 2)^(2/3) + 4)/(2*I*sqrt(15) + 2)^(1/3)) + 1/2*sqrt(-(2*I*sqrt(15) + 2)^(1/3) - 4/(2*I*sqrt(15) + 2)^(1/3) - 4/sqrt(((2*I*sqrt(15) + 2)^(2/3) + 4)/(2*I*sqrt(15) + 2)^(1/3)))]" ] }, "execution_count": 8, "metadata": { }, "output_type": "execute_result" } ], "source": [ "solve(x^4 + 2*x + 3 == 0, x) # There is an analogue of the quadratic formula for polynomials of degree 3 and 4..." ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "" ] }, "execution_count": 9, "metadata": { }, "output_type": "execute_result" } ], "source": [ "show(solve(x^4 + 2*x + 3 == 0, x)[0]) " ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "[x == -1/2*(1/18*sqrt(23)*sqrt(3) - 25/54)^(1/3)*(I*sqrt(3) + 1) - 1/18*(-I*sqrt(3) + 1)/(1/18*sqrt(23)*sqrt(3) - 25/54)^(1/3) + 1/3, x == -1/2*(1/18*sqrt(23)*sqrt(3) - 25/54)^(1/3)*(-I*sqrt(3) + 1) - 1/18*(I*sqrt(3) + 1)/(1/18*sqrt(23)*sqrt(3) - 25/54)^(1/3) + 1/3, x == (1/18*sqrt(23)*sqrt(3) - 25/54)^(1/3) + 1/9/(1/18*sqrt(23)*sqrt(3) - 25/54)^(1/3) + 1/3, x == -1/2*I*sqrt(3) - 1/2, x == 1/2*I*sqrt(3) - 1/2]" ] }, "execution_count": 10, "metadata": { }, "output_type": "execute_result" } ], "source": [ "solve(x^5+x+1 == 0, x) # ... but unless you get lucky..." ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "[0 == x^5 + 3*x + 1]" ] }, "execution_count": 11, "metadata": { }, "output_type": "execute_result" } ], "source": [ "solve(x^5+3*x+1 == 0, x) # ... not in degree 5. I'll cover this in Math 100C next quarter!" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "[x == 4]" ] }, "execution_count": 12, "metadata": { }, "output_type": "execute_result" } ], "source": [ "solve(sqrt(x) == 2, x)" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "[x == 0]" ] }, "execution_count": 13, "metadata": { }, "output_type": "execute_result" } ], "source": [ "solve(sin(x) == 0, x) # This doesn't give *all* of the solutions, but you can infer the rest." ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "[x == log(1/2*I*5^(1/3)*sqrt(3) - 1/2*5^(1/3)), x == log(-1/2*I*5^(1/3)*sqrt(3) - 1/2*5^(1/3)), x == 1/3*log(5)]" ] }, "execution_count": 14, "metadata": { }, "output_type": "execute_result" } ], "source": [ "solve(e^(3*x) == 5, x)" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "[{x: log(1/2*I*5^(1/3)*sqrt(3) - 1/2*5^(1/3))},\n", " {x: log(-1/2*I*5^(1/3)*sqrt(3) - 1/2*5^(1/3))},\n", " {x: 1/3*log(5)}]" ] }, "execution_count": 15, "metadata": { }, "output_type": "execute_result" } ], "source": [ "v = solve(e^(3*x) == 5, x, solution_dict=True); v" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "[[x == 1.35720887245841, y == -1.35720887245841], [x == (-0.6786044041487247 + 1.175377306225595*I), y == (0.6786044041487285 - 1.175377306225602*I)], [x == (-0.6786044041487247 - 1.175377306225595*I), y == (0.6786044041487285 + 1.175377306225602*I)]]" ] }, "execution_count": 16, "metadata": { }, "output_type": "execute_result" } ], "source": [ "var('x, y')\n", "solve([x^2 == y^2, x^3 == y^3 + 5], [x,y])" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "[{y: -1.35720887245841, x: 1.35720887245841},\n", " {y: 0.6786044041487285 - 1.175377306225602*I,\n", " x: -0.6786044041487247 + 1.175377306225595*I},\n", " {y: 0.6786044041487285 + 1.175377306225602*I,\n", " x: -0.6786044041487247 - 1.175377306225595*I}]" ] }, "execution_count": 17, "metadata": { }, "output_type": "execute_result" } ], "source": [ "\n", "solve([x^2 == y^2, x^3 == y^3 + 5], [x,y], solution_dict=True)" ] }, { "cell_type": "code", "execution_count": 26, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "55fe15f56ca8056e6d67e07cbd84d7f7f5b389af" } } ], "source": [ "g = implicit_plot(x^3 == y^3 + 5, (x, -3, 3), (y,-3,3))\n", "g += implicit_plot(x^2 == y^2, (x, -3, 3), (y,-3,3),color='red')\n", "g" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "" ] }, "execution_count": 18, "metadata": { }, "output_type": "execute_result" } ], "source": [ "show(v[0][x])" ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "e^(3*log(1/2*I*5^(1/3)*sqrt(3) - 1/2*5^(1/3)))" ] }, "execution_count": 20, "metadata": { }, "output_type": "execute_result" } ], "source": [ "e^(v[0][x]*3)" ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "5" ] }, "execution_count": 21, "metadata": { }, "output_type": "execute_result" } ], "source": [ "(e^(v[0][x]*3)).simplify_full()" ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "[(log(1/2*I*5^(1/3)*sqrt(3) - 1/2*5^(1/3)), 1),\n", " (log(-1/2*I*5^(1/3)*sqrt(3) - 1/2*5^(1/3)), 1),\n", " (1/3*log(5), 1)]" ] }, "execution_count": 22, "metadata": { }, "output_type": "execute_result" } ], "source": [ "# or use roots:\n", "v = (e^(3*x) - 5).roots()\n", "v" ] }, { "cell_type": "code", "execution_count": 30, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "" ] }, "execution_count": 30, "metadata": { }, "output_type": "execute_result" } ], "source": [ "show(v[0][0])" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "Here we numerically find a *single* zero of a function in a given interval." ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "776fbd131676e677dca065e97ecc9997e4fcf991" } } ], "source": [ "f(x) = e^(3*x) - 5\n", "plot(f, -2, 2, ymax=1)" ] }, { "cell_type": "code", "execution_count": 32, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "0.5364793041447001" ] }, "execution_count": 32, "metadata": { }, "output_type": "execute_result" } ], "source": [ "f.find_root(0, 1)" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "One can also turn a symbolic expression into a numerical approximation, to any desired accuracy (with some caution)." ] }, { "cell_type": "code", "execution_count": 24, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "500000000000000000000000000001/1000000000000000000000000000000*pi - 1/1000000000000000000000000000000*e" ] }, "execution_count": 24, "metadata": { }, "output_type": "execute_result" } ], "source": [ "f = pi/10^30 + pi - e/10^30 - acos(0)\n", "f" ] }, { "cell_type": "code", "execution_count": 25, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "1.57079632679490" ] }, "execution_count": 25, "metadata": { }, "output_type": "execute_result" } ], "source": [ "numerical_approx(f)" ] }, { "cell_type": "code", "execution_count": 27, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "" ] }, "execution_count": 27, "metadata": { }, "output_type": "execute_result" } ], "source": [ "alpha = log(1/2*I*5^(1/3)*sqrt(3) - 1/2*5^(1/3))\n", "show(alpha)" ] }, { "cell_type": "code", "execution_count": 28, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "0.536479304144700 + 2.09439510239320*I" ] }, "execution_count": 28, "metadata": { }, "output_type": "execute_result" } ], "source": [ "numerical_approx(alpha)" ] }, { "cell_type": "code", "execution_count": 29, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "0.53647930414470012486691977774206254650853378475617257397088 + 2.0943951023931954923084289221863352561314462662500705473166*I" ] }, "execution_count": 29, "metadata": { }, "output_type": "execute_result" } ], "source": [ "numerical_approx(alpha, prec=200)" ] }, { "cell_type": "code", "execution_count": 30, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "0.53647930414470012487 + 2.0943951023931954923*I" ] }, "execution_count": 30, "metadata": { }, "output_type": "execute_result" } ], "source": [ "numerical_approx(alpha, digits=20)" ] }, { "cell_type": "code", "execution_count": 41, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "0.536 + 2.09*I" ] }, "execution_count": 41, "metadata": { }, "output_type": "execute_result" } ], "source": [ "# This is the number of digits used in computing the result, NOT the number of correct digits in output!\n", "numerical_approx(alpha, digits=3)" ] }, { "cell_type": "code", "execution_count": 31, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "4.5000000000000000000000000000" ] }, "execution_count": 31, "metadata": { }, "output_type": "execute_result" } ], "source": [ "N(9/2, prec=100)" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "In order to obtain an answer which is *guaranteed* to have all reported digits be accurate (assuming no bugs in CoCalc, cosmic rays, or other irreproducible phenomena), one can use *interval arithmetic*. This is like keeping track of \"error bars\". It is somewhat less efficient, and so not commonly used in \"real-life\" applications; but it absolutely essentially if you want to use real-number arithmetic as the basis of a mathematical proof, as in the solution of the [Kepler conjecture](https://en.wikipedia.org/wiki/Kepler_conjecture)." ] }, { "cell_type": "code", "execution_count": 43, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 43, "metadata": { }, "output_type": "execute_result" } ], "source": [ "N is numerical_approx # Handy shorthand!" ] }, { "cell_type": "code", "execution_count": 0, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "alpha.N()" ] }, { "cell_type": "code", "execution_count": 0, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "alpha.n()" ] }, { "cell_type": "code", "execution_count": 0, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "alpha.numerical_approx()" ] }, { "cell_type": "code", "execution_count": 32, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "0.53648? + 2.09440?*I" ] }, "execution_count": 32, "metadata": { }, "output_type": "execute_result" } ], "source": [ "# Interval arithmetic --\n", "# Every displayed digit except last in the output is definitely right:\n", "ComplexIntervalField(20)(alpha)" ] }, { "cell_type": "code", "execution_count": 36, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "RIF = RealIntervalField()" ] }, { "cell_type": "code", "execution_count": 37, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "u = RIF(sqrt(2))+RIF(sqrt(3))" ] }, { "cell_type": "code", "execution_count": 38, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "(3.14626436994197, 3.14626436994198)" ] }, "execution_count": 38, "metadata": { }, "output_type": "execute_result" } ], "source": [ "u.lower(), u.upper()" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "## More about 2d plots\n", "\n", "You can do much more than just `plot(function...)`. E.g.,\n", "\n", "- a point\n", "- a bunch of points\n", "- text\n", "- \"line\" through a bunch of points\n", "- polygon\n", "- ellipse\n", "- implicit plot\n", "- contour plot\n", "- vector field" ] }, { "cell_type": "code", "execution_count": 44, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "a4d3465bf522b81cbe5408b8e2ead721a257d470" } } ], "source": [ "point((1,2))" ] }, { "cell_type": "code", "execution_count": 45, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "45f062c5831bec5cd8de4779d6d6a05bb2b40d7e" } } ], "source": [ "point2d((1,2), pointsize=150, color='red')" ] }, { "cell_type": "code", "execution_count": 40, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "64b911044c36659ab0b2eba3222c27fe4b93ad99" } } ], "source": [ "point([(random(), random()) for i in range(100)])" ] }, { "cell_type": "code", "execution_count": 47, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "1f0be5f9a6daa51a14e60d4b0af772a9cacd2c25" } } ], "source": [ "point([(x,x^2) for x in range(50)])" ] }, { "cell_type": "code", "execution_count": 50, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "02ef06f60b5017ab29f18f1243faffd3a3c4826e" } } ], "source": [ "g = point((1,2), pointsize=300, color='red')\n", "g += text(r\"Get $\\int_0^{\\pi} \\sin(x) dx$ points!\", (1,2), alpha=0.5, fontsize=30, fontweight='bold', color='black', rotation=20, zorder=1)\n", "g.show(frame=True, axes=False)" ] }, { "cell_type": "code", "execution_count": 51, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "c6c035164658f8c110c030a49e56a669e63fb61d" } } ], "source": [ "line([(0,0), (1,1), (1,2), (2,0), (3,1)], thickness=3, color='blue')" ] }, { "cell_type": "code", "execution_count": 52, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "685f5ac125f8ac38119bdafc97232e090779b78b" } } ], "source": [ "v = [(0,0), (1,1), (1,2), (2,0), (3,1)]\n", "line(v, thickness=3, color='blue', zorder=-1) + points(v, pointsize=200, color='grey', zorder=1)" ] }, { "cell_type": "code", "execution_count": 42, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "var('x, y')\n", "pl = implicit_plot(y^2 + cos(y*e^x) == x^3 - 2*x + 3, (x,-2,2), (y,-3,3))" ] }, { "cell_type": "code", "execution_count": 43, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "pl.save('test.pdf') # To preserve your handiwork!" ] }, { "cell_type": "code", "execution_count": 45, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": "\n\n" }, "execution_count": 45, "metadata": { }, "output_type": "execute_result" } ], "source": [ "var('x,y')\n", "plot3d(x^2+y^2, (x,-2,2),(y,-2,2))" ] }, { "cell_type": "code", "execution_count": 0, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ ] } ], "metadata": { "kernelspec": { "display_name": "SageMath 8.1", "name": "sage-8.1" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.14" } }, "nbformat": 4, "nbformat_minor": 0 }