{ "cells": [ { "cell_type": "markdown", "id": "6e55e2ba-804b-4a8c-a389-f96d71c4aa55", "metadata": {}, "source": [ "# Imports" ] }, { "cell_type": "code", "execution_count": 1, "id": "9d29038d-e09d-4f1e-8811-aeea5df85484", "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "\n", "import numpy as np\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "id": "5704d354-4e9f-420b-a085-244853e0cddc", "metadata": {}, "source": [ "# Home-made methods" ] }, { "cell_type": "markdown", "id": "adaf17d9-18d3-499c-81e9-56bfa5eaa4c7", "metadata": {}, "source": [ "## Gauss elimination" ] }, { "cell_type": "markdown", "id": "5c7ecd8f-b7f2-4759-afe5-bb259d15e496", "metadata": {}, "source": [ "It is illustrative to start writing our code working on a specific example, and then translate it into a function that accepts any matrix and vector" ] }, { "cell_type": "code", "execution_count": 11, "id": "22e345c8-f04b-46bf-941b-58866f08c4f6", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[1., 2., 3.],\n", " [0., 3., 4.],\n", " [0., 8., 9.]])" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a=np.array([\n", " [1,2,3],\n", " [1,3,4],\n", " [7,8,9]\n", "])\n", "b=np.array([9, 3, 4])\n", "\n", "# recast to avoid problems if a and b are integer ndarrays (as in this case)\n", "a = a.astype(np.float64)\n", "b = b.astype(np.float64)\n", "\n", "# find the coefficient to set to 0 the first element of the second row\n", "i=0\n", "a[i+1:,0] -= a[i+1:, i]/a[i, i] #the 0 in a[i+1:,0] is needed cause the RHS is just a vector, for now\n", "\n", "a" ] }, { "cell_type": "code", "execution_count": 12, "id": "d834b331-8712-4ec4-942d-74847dc993b3", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[1., 2., 3.],\n", " [0., 1., 1.],\n", " [7., 8., 9.]])" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a=np.array([\n", " [1,2,3],\n", " [1,3,4],\n", " [7,8,9]\n", "])\n", "b=np.array([9, 3, 4])\n", "\n", "a = a.astype(np.float64)\n", "b = b.astype(np.float64)\n", "\n", "# use that coefficient to manipulate the second row of the matrix\n", "# notice we dropped the aforementioned 0\n", "i=0\n", "a[i+1] -= a[i+1, i]/a[i, i] * a[i]\n", "\n", "a" ] }, { "cell_type": "code", "execution_count": 14, "id": "0498a8b8-404a-4b6a-9e15-810604a34805", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[ 1., 2., 3.],\n", " [ 0., 1., 1.],\n", " [ 0., -6., -12.]])" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a=np.array([\n", " [1,2,3],\n", " [1,3,4],\n", " [7,8,9]\n", "])\n", "b=np.array([9, 3, 4])\n", "\n", "a = a.astype(np.float64)\n", "b = b.astype(np.float64)\n", "\n", "# do the same to all the rows of the matrix after the first\n", "# by building a matrix from the outer product of an array of coefficients\n", "# and the first row of the matrix\n", "i=0\n", "a[i+1:] -= np.outer(a[i+1:, i]/a[i, i], a[i])\n", "\n", "a" ] }, { "cell_type": "code", "execution_count": 15, "id": "febed8ae-47ec-4c41-9deb-b15b401b7127", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "i= 0\n", "a=\n", " [[ 1. 2. 3.]\n", " [ 0. 1. 1.]\n", " [ 0. -6. -12.]]\n", "b=\n", " [ 9. -6. -59.]\n", "\n", "i= 1\n", "a=\n", " [[ 1. 2. 3.]\n", " [ 0. 1. 1.]\n", " [ 0. 0. -6.]]\n", "b=\n", " [ 9. -6. -95.]\n", "\n" ] }, { "data": { "text/plain": [ "array([[ 1., 2., 3.],\n", " [ 0., 1., 1.],\n", " [ 0., 0., -6.]])" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a=np.array([\n", " [1,2,3],\n", " [1,3,4],\n", " [7,8,9]\n", "])\n", "b=np.array([9, 3, 4])\n", "\n", "a = a.astype(np.float64)\n", "b = b.astype(np.float64)\n", "\n", "# repeat for all the rows of the matrix until the second to last\n", "# also apply the same transformations to the b vector\n", "for i in range(len(a)-1):\n", " b[i+1:] -= a[i+1:,i]/a[i,i] * b[i]\n", " a[i+1:] -= np.outer(a[i+1:,i]/a[i,i], a[i])\n", " print(\"i=\", i)\n", " print(\"a=\\n\", a)\n", " print(\"b=\\n\", b)\n", " print()\n", " \n", "a" ] }, { "cell_type": "code", "execution_count": 52, "id": "c52c61b0-bf62-45bd-9b7c-4619f5b052f6", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "i= 0\n", "a=\n", " [[ 1. 2. 3.]\n", " [ 0. 1. 1.]\n", " [ 0. -6. -12.]]\n", "b=\n", " [ 9. -6. -59.]\n", "\n", "i= 1\n", "a=\n", " [[ 1. 2. 3.]\n", " [ 0. 1. 1.]\n", " [ 0. 0. -6.]]\n", "b=\n", " [ 9. -6. -95.]\n", "\n" ] } ], "source": [ "a=np.array([\n", " [1,2,3],\n", " [1,3,4],\n", " [7,8,9]\n", "])\n", "b=np.array([9, 3, 4])\n", "\n", "a = a.astype(np.float64)\n", "b = b.astype(np.float64)\n", "\n", "# Notice that we did not bother to normalize all rows in such a way \n", "# that the first non-zero element is 1.\n", "# We will have to be mindful of this in the following.\n", "for i in range(len(a)-1):\n", " b[i+1:] -= a[i+1:,i]/a[i,i] * b[i]\n", " a[i+1:] -= np.outer(a[i+1:,i]/a[i,i], a[i])\n", " print(\"i=\", i)\n", " print(\"a=\\n\", a)\n", " print(\"b=\\n\", b)\n", " print()" ] }, { "cell_type": "code", "execution_count": 55, "id": "ce2b7476-447c-4a12-8ff3-947cdb41d8c4", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "5\n", "4\n", "3\n", "2\n", "1\n", "0\n", "\n", "2\n", "1\n", "0\n" ] } ], "source": [ "# Small digression on looping indeces backward, as required in the back substitution\n", "for i in range(5, -1, -1):\n", " print(i)\n", "\n", "print()\n", " \n", "for i in range(len(b)-1, -1, -1):\n", " print(i)" ] }, { "cell_type": "code", "execution_count": 57, "id": "a8652c57-b5f2-4f12-8ce8-c3f3746c8a42", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "starting backsub\n", "[0. 0. 0. 0.]\n", "[0. 0. 0. 0.]\n", "[0. 0. 1.66666667 0. ]\n", "[0. 0. 1.66666667 0. ]\n", "\n", "result\n", "[1.66666667 0. 1.66666667 0. ]\n", "[ 5.16666667 -21.83333333 15.83333333]\n" ] } ], "source": [ "# Implement the back substitution\n", "print(\"starting backsub\")\n", "x = np.zeros(len(b))\n", "for i in range(len(b)-1, -1, -1):\n", " print(x)\n", " x[i] = (b[i] - np.sum(x*a[i]))/a[i,i]\n", "print()\n", "\n", "# Compare our solution with a different implementation\n", "print(\"result\")\n", "print(x)\n", "\n", "a=np.array([\n", " [1,2,3],\n", " [1,3,4],\n", " [7,8,9]\n", "])\n", "b=np.array([9, 3, 4])\n", "print(np.linalg.solve(a,b))" ] }, { "cell_type": "code", "execution_count": 58, "id": "3c520202-2c6b-4cd5-af8a-b66d2064a23f", "metadata": {}, "outputs": [], "source": [ "# Start packing what we just did in a function\n", "def gauss_solve(a, b):\n", " a = a.astype(np.float64)\n", " b = b.astype(np.float64)\n", "\n", " for i in range(len(a)-1): \n", " b[i+1:] -= a[i+1:,i]/a[i,i] * b[i]\n", " a[i+1:] -= np.outer(a[i+1:,i]/a[i,i], a[i])\n", " print(\"i=\", i)\n", " print(\"a=\\n\", a)\n", " print(\"b=\\n\", b)\n", " print()\n", " \n", " print(\"starting backsub\")\n", " x = np.zeros(len(b))\n", " for i in range(len(b)-1, -1, -1):\n", " print(x)\n", " x[i] = (b[i] - np.sum(x*a[i]))/a[i,i]\n", " print()\n", "\n", " print(\"result\")\n", " return x" ] }, { "cell_type": "code", "execution_count": 59, "id": "064e8b78-4228-41bf-9843-c68e81af738b", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "i= 0\n", "a=\n", " [[ 1. 2. 3.]\n", " [ 0. 1. 1.]\n", " [ 0. -6. -12.]]\n", "b=\n", " [ 9. -6. -59.]\n", "\n", "i= 1\n", "a=\n", " [[ 1. 2. 3.]\n", " [ 0. 1. 1.]\n", " [ 0. 0. -6.]]\n", "b=\n", " [ 9. -6. -95.]\n", "\n", "starting backsub\n", "[0. 0. 0.]\n", "[ 0. 0. 15.83333333]\n", "[ 0. -21.83333333 15.83333333]\n", "\n", "result\n", "[ 5.16666667 -21.83333333 15.83333333]\n" ] } ], "source": [ "a=np.array([\n", " [1,2,3],\n", " [1,3,4],\n", " [7,8,9]\n", "])\n", "b=np.array([9, 3, 4])\n", "\n", "print(gauss_solve(a, b))" ] }, { "cell_type": "code", "execution_count": 60, "id": "f719ad51-b548-479f-b81e-57780495321b", "metadata": {}, "outputs": [], "source": [ "#Add some extra flags, and error management\n", "def gauss_solve(a, b, dbg=False):\n", " if(dbg):\n", " np_x = np.linalg.solve(a,b)\n", " \n", " a = a.astype(np.float64)\n", " b = b.astype(np.float64)\n", "\n", " for i in range(len(a)-1):\n", " if(a[i,i]==0):\n", " raise ValueError(\"0 encoutered on the diagonal. Use another solver.\")\n", " \n", " b[i+1:] -= a[i+1:,i]/a[i,i] * b[i]\n", " a[i+1:] -= np.outer(a[i+1:,i]/a[i,i], a[i])\n", " \n", " if(dbg):\n", " print(\"i=\", i)\n", " print(\"a=\\n\", a)\n", " print(\"b=\\n\", b)\n", " print()\n", " \n", " if(dbg):\n", " print(\"starting backsub\")\n", " print()\n", " \n", " x = np.zeros(len(b))\n", " for i in range(len(b)-1, -1, -1):\n", " if(dbg):\n", " print(x)\n", " x[i] = (b[i] - np.sum(x*a[i]))/a[i,i]\n", " if(dbg):\n", " print(x)\n", " print()\n", " return x, np_x\n", "\n", " return x" ] }, { "cell_type": "code", "execution_count": 61, "id": "31b135ef-acfd-4ad3-8bfa-3ae6212c508a", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[ 5.16666667 -21.83333333 15.83333333]\n", "[ 5.16666667 -21.83333333 15.83333333]\n", "[ 5.16666667 -21.83333333 15.83333333]\n" ] } ], "source": [ "a=np.array([\n", " [1,2,3],\n", " [1,3,4],\n", " [7,8,9]\n", "])\n", "b=np.array([9, 3, 4])\n", "\n", "print(gauss_solve(a, b))\n", "print(gauss_solve(a, b, False))\n", "print(gauss_solve(a, b, dbg=False))" ] }, { "cell_type": "code", "execution_count": 46, "id": "f08745ed-ea3c-4fa0-8302-e08328aef49d", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "i= 0\n", "a=\n", " [[ 1. 2. 3.]\n", " [ 0. 1. 1.]\n", " [ 0. -6. -12.]]\n", "b=\n", " [ 9. -6. -59.]\n", "\n", "i= 1\n", "a=\n", " [[ 1. 2. 3.]\n", " [ 0. 1. 1.]\n", " [ 0. 0. -6.]]\n", "b=\n", " [ 9. -6. -95.]\n", "\n", "starting backsub\n", "\n", "[0. 0. 0.]\n", "[ 0. 0. 15.83333333]\n", "[ 0. -21.83333333 15.83333333]\n", "[ 5.16666667 -21.83333333 15.83333333]\n", "\n" ] }, { "data": { "text/plain": [ "(array([ 5.16666667, -21.83333333, 15.83333333]),\n", " array([ 5.16666667, -21.83333333, 15.83333333]))" ] }, "execution_count": 46, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a=np.array([\n", " [1,2,3],\n", " [1,3,4],\n", " [7,8,9]\n", "])\n", "b=np.array([9, 3, 4])\n", "\n", "gauss_solve(a, b, True)\n" ] }, { "cell_type": "markdown", "id": "4c22dac6-aa28-4244-8a20-69a73c7a36b0", "metadata": {}, "source": [ "Not all matrices can be treated with this simple algorithm" ] }, { "cell_type": "code", "execution_count": 49, "id": "0bc405b6-982b-4d89-8991-15b7856061da", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "i= 0\n", "a=\n", " [[ 1. 2. 3.]\n", " [ 0. 0. 1.]\n", " [ 0. -6. -12.]]\n", "b=\n", " [ 9. -6. -59.]\n", "\n" ] }, { "ename": "ValueError", "evalue": "0 encoutered on the diagonal. Use another solver.", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mValueError\u001b[0m Traceback (most recent call last)", "Cell \u001b[0;32mIn[49], line 8\u001b[0m\n\u001b[1;32m 1\u001b[0m a\u001b[38;5;241m=\u001b[39mnp\u001b[38;5;241m.\u001b[39marray([\n\u001b[1;32m 2\u001b[0m [\u001b[38;5;241m1\u001b[39m,\u001b[38;5;241m2\u001b[39m,\u001b[38;5;241m3\u001b[39m],\n\u001b[1;32m 3\u001b[0m [\u001b[38;5;241m1\u001b[39m,\u001b[38;5;241m2\u001b[39m,\u001b[38;5;241m4\u001b[39m],\n\u001b[1;32m 4\u001b[0m [\u001b[38;5;241m7\u001b[39m,\u001b[38;5;241m8\u001b[39m,\u001b[38;5;241m9\u001b[39m]\n\u001b[1;32m 5\u001b[0m ])\n\u001b[1;32m 6\u001b[0m b\u001b[38;5;241m=\u001b[39mnp\u001b[38;5;241m.\u001b[39marray([\u001b[38;5;241m9\u001b[39m, \u001b[38;5;241m3\u001b[39m, \u001b[38;5;241m4\u001b[39m])\n\u001b[0;32m----> 8\u001b[0m \u001b[43mgauss_solve\u001b[49m\u001b[43m(\u001b[49m\u001b[43ma\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mb\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43;01mTrue\u001b[39;49;00m\u001b[43m)\u001b[49m\n", "Cell \u001b[0;32mIn[44], line 11\u001b[0m, in \u001b[0;36mgauss_solve\u001b[0;34m(a, b, dbg)\u001b[0m\n\u001b[1;32m 9\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m i \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mrange\u001b[39m(\u001b[38;5;28mlen\u001b[39m(a)\u001b[38;5;241m-\u001b[39m\u001b[38;5;241m1\u001b[39m):\n\u001b[1;32m 10\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m(a[i,i]\u001b[38;5;241m==\u001b[39m\u001b[38;5;241m0\u001b[39m):\n\u001b[0;32m---> 11\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mValueError\u001b[39;00m(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m0 encoutered on the diagonal. Use another solver.\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n\u001b[1;32m 13\u001b[0m b[i\u001b[38;5;241m+\u001b[39m\u001b[38;5;241m1\u001b[39m:] \u001b[38;5;241m-\u001b[39m\u001b[38;5;241m=\u001b[39m a[i\u001b[38;5;241m+\u001b[39m\u001b[38;5;241m1\u001b[39m:,i]\u001b[38;5;241m/\u001b[39ma[i,i] \u001b[38;5;241m*\u001b[39m b[i]\n\u001b[1;32m 14\u001b[0m a[i\u001b[38;5;241m+\u001b[39m\u001b[38;5;241m1\u001b[39m:] \u001b[38;5;241m-\u001b[39m\u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39mouter(a[i\u001b[38;5;241m+\u001b[39m\u001b[38;5;241m1\u001b[39m:,i]\u001b[38;5;241m/\u001b[39ma[i,i], a[i])\n", "\u001b[0;31mValueError\u001b[0m: 0 encoutered on the diagonal. Use another solver." ] } ], "source": [ "a=np.array([\n", " [1,2,3],\n", " [1,2,4],\n", " [7,8,9]\n", "])\n", "b=np.array([9, 3, 4])\n", "\n", "gauss_solve(a, b, True)" ] }, { "cell_type": "markdown", "id": "c700bc2d-6a81-4cc6-8726-6fc5e08eeea2", "metadata": {}, "source": [ "## Gauss with partial pivoting" ] }, { "cell_type": "code", "execution_count": 65, "id": "616b73cf-8ffd-44b3-b40f-7f641d9f6b74", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[1, 2, 3],\n", " [7, 8, 9],\n", " [1, 2, 4]])" ] }, "execution_count": 65, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# small digression on advanced indexing, which we use to swap the position of rows\n", "# when we implement the pivoting\n", "# https://numpy.org/doc/stable/user/basics.indexing.html\n", "a=np.array([\n", " [1,2,3],\n", " [1,2,4],\n", " [7,8,9]\n", "])\n", "\n", "a[[1, 2]] = a[[2, 1]]\n", "\n", "a" ] }, { "cell_type": "code", "execution_count": 68, "id": "91caa538-5416-4af5-a1f1-e52603875d92", "metadata": {}, "outputs": [], "source": [ "def gauss_solve_pivot(a, b, dbg=False):\n", " if(dbg):\n", " np_x = np.linalg.solve(a,b)\n", " \n", " a = a.astype(np.float64)\n", " b = b.astype(np.float64)\n", "\n", " for i in range(len(a)-1):\n", " # Notice the \"i + \". It is there cause len(a[i:,i]) = len(a) - i\n", " new_i = i + np.argmax(a[i:,i]) \n", " a[[i, new_i]] = a[[new_i, i]]\n", " b[[i, new_i]] = b[[new_i, i]]\n", " \n", " b[i+1:] -= a[i+1:,i]/a[i,i] * b[i]\n", " a[i+1:] -= np.outer(a[i+1:,i]/a[i,i], a[i])\n", " \n", " if(dbg):\n", " print(\"i=\", i)\n", " print(\"a=\\n\", a)\n", " print(\"b=\\n\", b)\n", " print()\n", " \n", " if(dbg):\n", " print(\"starting backsub\")\n", " \n", " x = np.zeros(len(b))\n", " for i in range(len(b)-1, -1, -1):\n", " x[i] = (b[i] - np.sum(x*a[i]))/a[i,i]\n", " \n", " if(dbg):\n", " return x, np_x\n", "\n", " return x" ] }, { "cell_type": "code", "execution_count": 69, "id": "72a2163b-d74d-4291-a12e-26b256c630a0", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "i= 0\n", "a=\n", " [[7. 8. 9. ]\n", " [0. 0.85714286 2.71428571]\n", " [0. 0.85714286 1.71428571]]\n", "b=\n", " [4. 2.42857143 8.42857143]\n", "\n", "i= 1\n", "a=\n", " [[ 7. 8. 9. ]\n", " [ 0. 0.85714286 2.71428571]\n", " [ 0. 0. -1. ]]\n", "b=\n", " [4. 2.42857143 6. ]\n", "\n", "starting backsub\n" ] }, { "data": { "text/plain": [ "(array([-16.66666667, 21.83333333, -6. ]),\n", " array([-16.66666667, 21.83333333, -6. ]))" ] }, "execution_count": 69, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a=np.array([\n", " [1,2,3],\n", " [1,2,4],\n", " [7,8,9]\n", "])\n", "b=np.array([9, 3, 4])\n", "\n", "gauss_solve_pivot(a, b, True)" ] }, { "cell_type": "markdown", "id": "cb202f29-6a3e-4c24-b226-8eb59ccdce45", "metadata": {}, "source": [ "It is now trivial to solve the exercise on slide 44" ] }, { "cell_type": "code", "execution_count": 71, "id": "84639ce6-3488-44f5-8b3e-59ac5ce08e13", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([3. , 1.66666667, 3.33333333, 2. ])" ] }, "execution_count": 71, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a=np.array([\n", " [ 4,-1,-1,-1],\n", " [-1, 3, 0,-1],\n", " [-1, 0, 3,-1],\n", " [-1,-1,-1, 4],\n", "])\n", "b=np.array([5, 0, 5, 0])\n", "\n", "gauss_solve_pivot(a, b)" ] }, { "cell_type": "markdown", "id": "e6a3d5ea-ebb8-4703-b318-030d6e070541", "metadata": {}, "source": [ "## numpy ndarray and linalg" ] }, { "cell_type": "code", "execution_count": 85, "id": "eb8339b3-57bd-4fbd-b646-0021741e2729", "metadata": {}, "outputs": [], "source": [ "#ndarrays are defined as\n", "a=np.array([\n", " [1,2,3],\n", " [1,2,4],\n", " [7,8,9]\n", "])\n", "b=np.array([5, 0, 5,])" ] }, { "cell_type": "code", "execution_count": 74, "id": "5907c843-8979-4e77-b51d-b1ffd3a8b842", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[1, 2, 3],\n", " [1, 2, 4],\n", " [7, 8, 9]])" ] }, "execution_count": 74, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a" ] }, { "cell_type": "code", "execution_count": 75, "id": "e37b3290-4806-498e-bfb8-0bc24c16dab0", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([2, 2, 8])" ] }, "execution_count": 75, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# you can slice them\n", "a[:,1]" ] }, { "cell_type": "code", "execution_count": 82, "id": "619e1fa6-82e5-4560-b2b3-089e2d364581", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[1, 1, 7],\n", " [2, 2, 8],\n", " [3, 4, 9]])" ] }, "execution_count": 82, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Transpose matrix\n", "a.T" ] }, { "cell_type": "code", "execution_count": 87, "id": "25e5ac6f-6cfe-467a-b4b0-a9b500d4a172", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[5 0 5]\n", "[5 0 5]\n" ] } ], "source": [ "# you can't transpose a 1d vector\n", "print(b)\n", "print(b.T)" ] }, { "cell_type": "code", "execution_count": 92, "id": "bc027058-9b21-4696-969b-a1e512d90649", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1 2 3]\n", "[5 0 5]\n" ] } ], "source": [ "# a few types of multiplications are defined\n", "# make sure to use the correct one\n", "\n", "a = np.array([1, 2, 3])\n", "print(a)\n", "print(b)" ] }, { "cell_type": "code", "execution_count": 93, "id": "f74a711f-ec92-433d-8b0a-74b8c50103dc", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[ 5 0 15]\n" ] } ], "source": [ "# element-wise\n", "print(b*a)" ] }, { "cell_type": "code", "execution_count": 94, "id": "bfdf330b-0a16-4bb3-a795-a49e6524477c", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "20" ] }, "execution_count": 94, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# inner\n", "np.dot(b,a)" ] }, { "cell_type": "code", "execution_count": 95, "id": "e230148e-81fb-499e-ae5a-199052997d79", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[ 5, 10, 15],\n", " [ 0, 0, 0],\n", " [ 5, 10, 15]])" ] }, "execution_count": 95, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# inner\n", "np.outer(b,a)" ] }, { "cell_type": "code", "execution_count": 97, "id": "320dc4a7-ae83-4aad-b241-6beea5553ab0", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([20, 25, 80])" ] }, "execution_count": 97, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# matrix multiplication\n", "a=np.array([\n", " [1,2,3],\n", " [1,2,4],\n", " [7,8,9]\n", "])\n", "b=np.array([5, 0, 5,])\n", "\n", "a@b" ] }, { "cell_type": "code", "execution_count": 105, "id": "1a89dcc6-41c1-4666-9122-8be6013abc6d", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "20\n", "20\n", "20\n", "\n", "20\n", "[[ 5 0 5]\n", " [10 0 10]\n", " [15 0 15]]\n" ] } ], "source": [ "# since .T doens't work as one might naively expect on 1d arrays\n", "# the matrix product of two arrays might not be what one expects\n", "\n", "a = np.array([1, 2, 3])\n", "\n", "print(a@b)\n", "print(a.T@b)\n", "print(a@b.T)\n", "\n", "print()\n", "\n", "# some functions are defined for you\n", "print(np.dot(a,b))\n", "print(np.outer(a,b))" ] }, { "cell_type": "code", "execution_count": 114, "id": "b9b14912-7b29-4820-929e-4dd24cdbd0d5", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1d\n", " [1 2 3]\n", "column\n", " [[1]\n", " [2]\n", " [3]]\n", "row\n", " [[1 2 3]]\n" ] } ], "source": [ "# or you can recast the 1d arrays into 2d arrays\n", "print(\"1d\\n\", a)\n", "print(\"column\\n\", np.reshape(a, (3,1)))\n", "print(\"row\\n\", np.reshape(a, (1,3)))" ] }, { "cell_type": "code", "execution_count": 121, "id": "2ecf5444-b3bd-4a74-8ff8-389015b14db1", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[ 5 0 5]\n", " [10 0 10]\n", " [15 0 15]]\n", "[[20]]\n" ] } ], "source": [ "# and manipulate them as \"expected\"\n", "print(np.reshape(a, (3,1)) @ np.reshape(b, (1,3)))\n", "print(np.reshape(a, (1,3)) @ np.reshape(b, (3,1)))" ] }, { "cell_type": "code", "execution_count": 125, "id": "bfb465f9-bff8-423c-a048-8b7025fe295f", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[-2.33333333, 1. , 0.33333333],\n", " [ 3.16666667, -2. , -0.16666667],\n", " [-1. , 1. , -0. ]])" ] }, "execution_count": 125, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# linear algebra operations can be performed with linalg\n", "a=np.array([\n", " [1,2,3],\n", " [1,2,4],\n", " [7,8,9]\n", "])\n", "b=np.array([5, 0, 5,])\n", "\n", "\n", "# inverse\n", "np.linalg.inv(a)" ] }, { "cell_type": "code", "execution_count": 126, "id": "e22b5c10-dc77-42a3-9ab9-7b82622625c1", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "6.0" ] }, "execution_count": 126, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#determinant\n", "np.linalg.det(a)" ] }, { "cell_type": "code", "execution_count": 127, "id": "1a982662-2648-41b6-b1c2-e18155d45356", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(array([13.90136774, -0.26352478, -1.63784296]),\n", " array([[-0.26184976, -0.47726026, -0.26247798],\n", " [-0.32716447, 0.810583 , -0.68029833],\n", " [-0.90796372, -0.33937861, 0.68432412]]))" ] }, "execution_count": 127, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# extract eigenvectors and eigenvalues\n", "np.linalg.eig(a)" ] }, { "cell_type": "code", "execution_count": 128, "id": "5d9d4133-8b0f-4f27-ae40-00bf21e18e19", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[13.90136774 -0.26352478 -1.63784296]\n" ] }, { "data": { "text/plain": [ "6.00000000000001" ] }, "execution_count": 128, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# or just the eigenvalues, and try to get the det from them\n", "print(np.linalg.eigvals(a))\n", "np.prod(np.linalg.eigvals(a))" ] }, { "cell_type": "code", "execution_count": 129, "id": "5aaf2870-2102-41e5-abe0-238af6458a52", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([-10., 15., -5.])" ] }, "execution_count": 129, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Solve linear problems (without re-implementing the algorithms yourself)\n", "np.linalg.solve(a,b)" ] } ], "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.8.10" } }, "nbformat": 4, "nbformat_minor": 5 }