{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "836f12a9f8fd45c489bd987d0f18b087", "version_major": 2, "version_minor": 0 }, "text/plain": [ "SW50ZXJhY3RpdmUgZnVuY3Rpb24gPGZ1bmN0aW9uIGZsb3dlciBhdCAweDEwNTFiMDg0OD4gd2l0aCAxIHdpZGdldAogIGVuZHQ6IFRyYW5zZm9ybUZsb2F0U2xpZGVyKHZhbHVlPTAuMTU3MDfigKY=\n" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# 1. Consider the curve defined by r = sin(8*t).\n", "# Make an interact to plot this curve.\n", "# The range for t always starts at zero.\n", "# The end of the range for t is controlled by a slider.\n", "# The initial value for the end is pi/2.\n", "# The increment for the end value is pi/40.\n", "@interact\n", "def flower(endt = slider(0,2*pi,pi/40,pi/20,label='angle')):\n", " t = var('t')\n", " P = polar_plot(sin(8*t),(t,0,endt))\n", " P.show(xmin=-1,xmax=1,ymin=-1,ymax=1,figsize=4)" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[5 1 0 0 0]\n", "[1 5 1 0 0]\n", "[0 1 5 1 0]\n", "[0 0 1 5 1]\n", "[0 0 0 1 5] \n", "[[5 1 0 0 0]\n", " [1 5 1 0 0]\n", " [0 1 5 1 0]\n", " [0 0 1 5 1]\n", " [0 0 0 1 5]] \n" ] } ], "source": [ "# 2. Use numpy to solve a 5-by-5 tridiagonal system A*x = b.\n", "# (1) The diagonal element of A is 5, the elements just above \n", "# and below the diagonal are one.\n", "# Everywhere else the matrix is zero.\n", "# (2) Define a 5-dimensional right hand side vector b of ones.\n", "# Solve the system A*x = b and compute the residual.\n", "# Answer to (1)\n", "reset()\n", "a = lambda i,j: 5 if i == j else (1 if abs(i-j) == 1 else 0)\n", "m = matrix(5,5,a)\n", "print m, type(m)\n", "import numpy as np\n", "A = np.matrix(m)\n", "print A, type(A)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[1]\n", " [1]\n", " [1]\n", " [1]\n", " [1]]\n", "[[ 0.17272727]\n", " [ 0.13636364]\n", " [ 0.14545455]\n", " [ 0.13636364]\n", " [ 0.17272727]]\n", "2.48253415325e-16\n" ] } ], "source": [ "# Answer to (2)\n", "b = np.matrix([1 for k in range(5)])\n", "b = b.transpose()\n", "print b\n", "x = np.linalg.solve(A,b)\n", "print x\n", "r = b - A*x\n", "print np.linalg.norm(r)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "-2.2704074859237847541179172835354872824\n", "-2988163778/638512875\n" ] } ], "source": [ "# 3. Use the Cauchy integral formula to compute the number of\n", "# complex roots in a disk centered at 0 and with radius 1.1\n", "# of (x+1)*sin(2 x).\n", "# Give the number of roots in that disk of the complex plane.\n", "# Write all relevant PARI/GP commands.\n", "# ANSWER: the sin() function is represented by a power series\n", "# and we have to truncate the series for the subspol() command.\n", "gp('f(x) = (x+1)*sin(2*x)')\n", "gp('df = deriv(f(z),z)')\n", "gp('g(x) = substpol(truncate(df),z, x)')\n", "print gp('f(2)')\n", "print gp('g(2)')" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1.9999869691649256131596720954138006077 - 4.583818152745797615 E-40*I\n", "the number of roots : 2\n" ] } ], "source": [ "n = gp('intcirc(z=0,1.1,g(z)/f(z))')\n", "print n\n", "print 'the number of roots :', round(abs(n.sage()))" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Ideal (x^2*y - 2*x + 3, x*y^2 - 2*y + 3) of Multivariate Polynomial Ring in x, y over Rational Field\n", "y^3-2*y+3,\n", "x-y\n" ] } ], "source": [ "# 4. Consider the polynomial system defined by the polynomials\n", "# p = x^2*y - 2*x + 3 and q = x*y^2 - 2*y + 3.\n", "# (1) Bring the system in triangular form. \n", "# Use this triangular form to determine the number of complex solutions.\n", "# (2) If N is the number of solutions, compute the companion matrix of the system.\n", "# The rows of this matrix are the reductions of\n", "# the products of y with y^k for k ranging between 0 and N-1.\n", "# Show that the y-coordinates of the solutions are the eigenvalues of this companion matrix.\n", "reset()\n", "R. = PolynomialRing(QQ, order='lex')\n", "p = x^2*y - 2*x + 3\n", "q = x*y^2 - 2*y + 3\n", "J = Ideal(p,q)\n", "print J\n", "g = singular.groebner(J)\n", "print g" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[y, y^2, 2*y-3]\n" ] } ], "source": [ "# We see that there are three solutions and compute the reductions with respect to g.\n", "r = [singular.reduce(y*y^k,g) for k in range(3)]\n", "print r" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[y, y^2, 2*y - 3]\n", "[[0, 1, 0], [0, 0, 1], [-3, 2, 0]]\n", "[0.000000000000000 1.00000000000000 0.000000000000000]\n", "[0.000000000000000 0.000000000000000 1.00000000000000]\n", "[-3.00000000000000 2.00000000000000 0.000000000000000]\n", "[-1.89328919630450, 0.946644598152249 - 0.829703552862405*I, 0.946644598152249 + 0.829703552862405*I]\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "/Users/jan/sage-8.3/local/lib/python2.7/site-packages/sage/all_cmdline.py:9: UserWarning: Using generic algorithm for an inexact ring, which will probably give incorrect results due to numerical precision issues.\n", " #\n" ] } ], "source": [ "# For the companion matrix we use Sage.\n", "sr = [p.sage() for p in r]\n", "print sr\n", "m = [[e.coefficient({y:k}) for k in range(3)] for e in sr]\n", "print m\n", "import numpy as np\n", "M = Matrix(RR,m)\n", "print M\n", "print M.eigenvalues()" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[y^3-2*y+3, x-y]\n", "y^3 - 2*y + 3 \n", "[(-1.89328919630450, 1), (0.946644598152249 - 0.829703552862405*I, 1), (0.946644598152249 + 0.829703552862405*I, 1)]\n" ] } ], "source": [ "# Let us verify with the first element of the Groebner basis,\n", "# which is a polynomial in y. Notice the type conversions...\n", "lg = list(g)\n", "print lg\n", "g0 = lg[0].sage()\n", "print g0, type(g0)\n", "g0y = CC['y'](g0)\n", "print g0y.roots()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "SageMath 8.3", "language": "", "name": "sagemath" }, "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.15" } }, "nbformat": 4, "nbformat_minor": 2 }