{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Scipy\n", "\n", "[Scipy](https://docs.scipy.org/doc/scipy/tutorial/general.html) provides many more computational algorithms and is built on Numpy." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%config InlineBackend.figure_format = 'svg'\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import scipy.linalg as la" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Use LU decomposition to solve $Ax=b$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Create a random $5 \\times 5$ matrix" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[0.29145244 0.35884437 0.55784979 0.46054394 0.32163474]\n", " [0.84136821 0.83746299 0.24546883 0.5487101 0.67319653]\n", " [0.26649574 0.57778323 0.48906678 0.22667758 0.06552584]\n", " [0.96896256 0.83054189 0.37846023 0.02328101 0.5339853 ]\n", " [0.42419362 0.16210889 0.44353171 0.55697708 0.9857654 ]]\n" ] } ], "source": [ "n = 5\n", "A = np.random.rand(n,n)\n", "print(A)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Compute its LU factorization with pivoting using [scipy.linalg.lu_factor](https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.lu_factor.html)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "lu,piv = la.lu_factor(A)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Another option is to use [scipy.linalg.lu](https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.lu.html) which returns different arguments.\n", "\n", "Create a right hand side vector" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "b = np.random.rand(n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Solve $A x = b$" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[ 0.80575593 -0.8105171 1.30349815 -0.95036498 0.71113587]\n" ] } ], "source": [ "x = la.lu_solve((lu,piv),b)\n", "print(x)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Check that $x$ solves the problem by computing $Ax-b$" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[-5.55111512e-17 0.00000000e+00 -5.55111512e-17 -1.11022302e-16\n", " -1.11022302e-16]\n" ] } ], "source": [ "print(A@x-b)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If we do not want the LU decomposition, we can directly solve, using [scipy.linalg.solve](https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.solve.html)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[ 0.80575593 -0.8105171 1.30349815 -0.95036498 0.71113587]\n" ] } ], "source": [ "y = la.solve(A,b)\n", "print(y)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Note**: There are similar methods in Numpy, see [numpy.linalg.solve](https://numpy.org/doc/stable/reference/generated/numpy.linalg.solve.html)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Sparse matrix" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Scipy provides methods that can work on sparse matrices, see [scipy.sparse.linalg](https://docs.scipy.org/doc/scipy/reference/sparse.linalg.html)\n", "\n", "Sparse matrix formats are provided by\n", "\n", " * [csc_matrix](https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.csc_matrix.html)\n", " * [csr_matrix](https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.csr_matrix.html)\n", " * [coo_matrix](https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.coo_matrix.html)\n", " \n", " csc_matrix stores the entries column-wise." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", " Coords\tValues\n", " (0, 0)\t3.0\n", " (1, 0)\t1.0\n", " (0, 1)\t2.0\n", " (1, 1)\t-1.0\n", " (2, 1)\t5.0\n", " (2, 2)\t1.0\n", "[[ 3. 2. 0.]\n", " [ 1. -1. 0.]\n", " [ 0. 5. 1.]]\n" ] } ], "source": [ "from scipy.sparse import csc_matrix, csr_matrix\n", "A = csc_matrix([[3, 2, 0], \n", " [1, -1, 0], \n", " [0, 5, 1]], dtype=float)\n", "print(A)\n", "print(A.todense())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "csr_matrix stores them row-wise." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", " Coords\tValues\n", " (0, 0)\t3.0\n", " (0, 1)\t2.0\n", " (1, 0)\t1.0\n", " (1, 1)\t-1.0\n", " (2, 1)\t5.0\n", " (2, 2)\t1.0\n", "[[ 3. 2. 0.]\n", " [ 1. -1. 0.]\n", " [ 0. 5. 1.]]\n" ] } ], "source": [ "A = csr_matrix([[3, 2, 0], \n", " [1, -1, 0], \n", " [0, 5, 1]], dtype=float)\n", "print(A)\n", "print(A.todense())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can construct sparse matrix by giving indices and values of non-zero entries. Indices not specified are assumed to be zero and need not be stored." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", " Coords\tValues\n", " (0, 0)\t3.0\n", " (1, 0)\t1.0\n", " (0, 1)\t2.0\n", " (1, 1)\t-1.0\n", " (2, 1)\t5.0\n", " (2, 2)\t1.0\n", "[[ 3. 2. 0.]\n", " [ 1. -1. 0.]\n", " [ 0. 5. 1.]]\n" ] } ], "source": [ "row = np.array([0, 0, 1, 1, 2, 2])\n", "col = np.array([0, 1, 0, 1, 1, 2])\n", "data = np.array([3.0, 2.0, 1.0, -1.0, 5.0, 1.0])\n", "A = csc_matrix((data, (row, col)), shape=(3, 3))\n", "print(A)\n", "print(A.todense())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Curve fitting\n", "\n", "Suppose we are given some data set \n", "\n", "$$\n", "(x_i,y_i), \\qquad i=0,1,...,n-1\n", "$$\n", "\n", "and we suspect there is a linear relation ship between them\n", "\n", "$$\n", "y = a + b x\n", "$$\n", "\n", "But possibly the data also has some noise. Let us generate such a data set." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", " \n", " 2025-01-05T14:50:43.805237\n", " image/svg+xml\n", " \n", " \n", " Matplotlib v3.10.0, https://matplotlib.org/\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n" ], "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Sample points\n", "n = 10\n", "x = np.linspace(0,1,n)\n", "\n", "# Exact data\n", "a = 1.0\n", "b = 2.0\n", "ye = a + b*x\n", "\n", "# Add some noise\n", "y = ye + 0.1*(2*np.random.rand(n)-1)\n", "plt.plot(x,y,'o',label='Noisy Data')\n", "plt.plot(x,ye,label='True function')\n", "plt.legend();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Define the function we want to fit. The first argument is the independent variable, and remaining are the parameters that we want to find." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "def f(x,a,b):\n", " return a + b*x" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can use [curve_fit](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html) function from scipy to do the fitting." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from scipy.optimize import curve_fit\n", "popt, pcov = curve_fit(f, x, y)\n", "print('Fitted parameters = ',*popt)\n", "plt.plot(x,y,'o',label='Noisy Data')\n", "plt.plot(x,f(x,*popt),label='Fitted function')\n", "plt.plot(x,ye,label='True function')\n", "plt.legend();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can also use the [polyfit](https://numpy.org/doc/stable/reference/generated/numpy.polyfit.html) function or [fit](https://numpy.org/doc/stable/reference/generated/numpy.polynomial.polynomial.Polynomial.fit.html#numpy.polynomial.polynomial.Polynomial.fit) function from numpy to do this. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "c = np.polyfit(x, y, deg=1)\n", "print('Polynomial coefficients = ',c)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The polynomials coefficients are returned in $c$ in the following order\n", "\n", "$$\n", "p(x) = c[0] * x^{deg} + c[1] * x^{deg-1} + \\ldots + c[-1]\n", "$$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "p = np.poly1d(c)\n", "plt.plot(x,y,'o',label='Noisy Data')\n", "plt.plot(x,p(x),label='Fitted function')\n", "plt.plot(x,ye,label='True function')\n", "plt.legend();" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "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.12.8" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": {}, "toc_section_display": true, "toc_window_display": false }, "vscode": { "interpreter": { "hash": "c6e4e9f98eb68ad3b7c296f83d20e6de614cb42e90992a65aa266555a3137d0d" } } }, "nbformat": 4, "nbformat_minor": 4 }