{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## Linear algebra\n", "\n", "#### scipy.linalg documentation at https://docs.scipy.org/doc/scipy/reference/linalg.html" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Gaussian elimination" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# implementing my own code to solve a linear system with \n", "# Gaussian elimination (no pivoting)\n", "\n", "import numpy as np\n", "from numpy.linalg import solve\n", "\n", "def linear_solve_gauss(A,b):\n", "\n", " # Solve linear system Ax = b via \n", " # Gauss elimination\n", "\n", " demomode = False\n", " \n", " N = len(b)\n", " \n", " for i in range(N):\n", " #print(i)\n", " \n", " # Normalize diagonal to 1.\n", " fact = A[i,i]\n", " A[i,:] /= fact\n", " b[i] /= fact\n", " \n", " for j in range(i+1,N):\n", " # Subract upper rows, multiplied by suitable \n", " # factor, to generate an upper diagonal matrix\n", " # At the end of each loop, column i (set by loop above) \n", " # has become (1, 0, .... , 0) \n", " fact = A[j,i]\n", " A[j,:] -= fact*A[i,:]\n", " b[j] -= fact*b[i]\n", " \n", " if (demomode):\n", " print( ' ' )\n", " print('i = ', i)\n", " print('A=', '\\n', A)\n", " \n", " # Checking that I did it right\n", " if (demomode):\n", " print(' ')\n", " print('After Gauss elimination')\n", " print('A = ', '\\n', A)\n", " print( ' ' )\n", " print('b = ', b)\n", " print( ' ')\n", " \n", " \n", " # Final step: backsubstitution\n", " \n", " # Initialize solution vector\n", " x = np.zeros(len(b)) \n", " # Gets last solution, x_(N-1)\n", " x[N-1] = b[N-1]/A[N-1,N-1]\n", " \n", " # Proceed backward, line by line\n", " for i in range(N-2,-1,-1):\n", " v = np.dot(A[i,i+1:N],x[i+1:N])\n", " #if (i==1): print(i, i+1, N-1, A[i,i+1:N])\n", " x[i] = (b[i] - v)/A[i,i]\n", " \n", " return x" ] }, { "cell_type": "code", "execution_count": 111, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[ 2. 1. 4. 1.]\n", " [ 3. 4. -1. -1.]\n", " [ 1. -4. 1. 5.]\n", " [ 2. -2. 1. 3.]]\n" ] } ], "source": [ "# Create example for testing\n", "\n", "A = [ [2, 1, 4, 1], [3, 4, -1, -1], [1, -4, 1, 5], [2, -2, 1, 3]]\n", "A = np.array(A,dtype=float)\n", "print(A)\n", "\n", "b = [-4, 3, 9, 7]\n", "b = np.array(b,dtype=float)\n", "#print(b)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x = linear_solve_gauss(A,b)\n", "\n", "#print(A)\n", "\n", "print('The solution vector x is')\n", "print(x)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Gaussian elimination with pivoting" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# Adding pivoting to the previous code\n", "\n", "import numpy as np\n", "from numpy.linalg import solve\n", "\n", "def linsolve_gauss_pivot(A,b):\n", "\n", " # Solve linear system Ax = b via \n", " # Gauss elimination with pivoting\n", "\n", " demomode = False\n", " \n", " N = len(b)\n", " \n", " for i in range(N):\n", " \n", " # Partial pivoting\n", " ####################### \n", " # Collecting elements of column i in vector 'col'\n", " col = abs(A[:,i]) \n", " # Find max element in col (I can of course also perform this and \n", " # the previous step in a single line)\n", " ind = np.argmax(col)\n", " # If index of max is not the current i in the loop, swap rows\n", " if (ind > i):\n", " # Swap\n", " A[ [i,ind] ] = A[ [ind,i] ]\n", " b[ [i,ind] ] = b[ [ind,i] ]\n", " \n", " if demomode:\n", " print('\\n', 'After pivoting')\n", " print()\n", " print('i = ', i)\n", " print('A=', '\\n', A)\n", " print( ' ' )\n", " print('b = ', b)\n", " print( ' ')\n", " \n", " # Normalize diagonal to 1.\n", " fact = A[i,i]\n", " A[i,:] /= fact\n", " b[i] /= fact\n", " \n", " for j in range(i+1,N):\n", " # Subract upper rows, multiplied by suitable \n", " # factor, to generate an upper diagonal matrix\n", " # At the end of each loop, column i (set by loop above) \n", " # has become (1, 0, .... , 0) \n", " fact = A[j,i]\n", " A[j,:] -= fact*A[i,:]\n", " b[j] -= fact*b[i]\n", " \n", " # Checking that I did it right\n", " if demomode:\n", " print(' ')\n", " print('After Gauss elimination')\n", " print('A = ', '\\n', A)\n", " print( ' ' )\n", " print('b = ', b)\n", " print( ' ')\n", " \n", " \n", " # Final step: backsubstitution\n", " \n", " # Initialize solution vector\n", " x = np.zeros(len(b)) \n", " # Gets last solution, x_(N-1)\n", " x[N-1] = b[N-1]/A[N-1,N-1]\n", " \n", " # Proceed backward, line by line\n", " for i in range(N-2,-1,-1):\n", " v = np.dot(A[i,i+1:N],x[i+1:N])\n", " #if (i==1): print(i, i+1, N-1, A[i,i+1:N])\n", " x[i] = (b[i] - v)/A[i,i]\n", " \n", " return x\n" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# Testing example for pivoting\n", "\n", "A = [ [0, 1, 4, 1], [3, 4, -1, -1], [1, -4, 1, 5], [2, -2, 1, 3]]\n", "A = np.array(A,dtype=float)\n", "#print(A)\n", "\n", "b = [-4, 3, 9, 7]\n", "b = np.array(b,dtype=float)\n", "#print(b)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[ 1.61904762 -0.42857143 -1.23809524 1.38095238]\n" ] } ], "source": [ "x = linsolve_gauss_pivot(A,b)\n", "print(x)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print( (np.dot(A,x)-b) )" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[ 1.61904762 -0.42857143 -1.23809524 1.38095238]\n" ] } ], "source": [ "# Cross checking against python solve\n", "x = solve(A,b)\n", "print(x)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Another test\n", "\n", "A = [ [0, 1, 4, 1], [3, 0, -1, -1], [1, -4, 1, 5], [2, -2, 1, 0]]\n", "A = np.array(A,dtype=float)\n", "#print(A)\n", "\n", "b = [-4, 3, 9, 7]\n", "b = np.array(b,dtype=float)\n", "#print(b)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x = linsolve_gauss_pivot(A,b)\n", "print(x)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x = solve(A,b)\n", "print(x)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#Check with 5x5 array\n", "\n", "A = [ [0, 1, 4, 1, 7], [3, 0, -1, -1, 3], [1, -4, 1, 5, 0], [2, -2, 1, 0, -5], [-6, 3, 12, -10, 4]]\n", "A = np.array(A,dtype=float)\n", "#print(A)\n", "\n", "b = [-4, 3, 9, 7, 8]\n", "b = np.array(b,dtype=float)\n", "#print(b)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x = solve(A,b)\n", "print(x)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "y = linsolve_gauss_pivot(A,b)\n", "print(y)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### LU decomposition" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Implementing LU decomposition of an array A. \n", "\n", "A = LU with L = lower triangular, U = upper triangular.\n", "\n", "No need to restart from scratch; it can be done by modifying the previous linear solver \n", "code. L and U are generated using formulae we saw in class, in the same loop used to update \n", "A with Gauss elimination.\n", "\n", "#### Note below if you want to try and implement LU with pivoting \n", "\n", "Important: when we store L and U, we need also to record all the swaps in the pivoting, because the same swaps need to applied to b later on!\n", "\n", "Consider for example a 4 x 4 matrix (generalization to NxN is straightforward). Assume you apply partial pivoting to compute U. Represent the row swaps with permutation\n", "matrices P (for later remember that $P^2 = P$).\n", "\n", "Then we can write (with the L matrices as defined in class, but we don't normalize the diagonal of U to 1):\n", "\n", "$L_2 P_2 L_1 P_1 L_0 P_0 A = U$\n", "\n", "1) Show that you can write:\n", "\n", "$L'_2 L'_1 L'_0 P_2 P_1 P_0 = L_2 P_2 L_1 P_1 L_0 P_0$,\n", "\n", "where\n", "\n", "$L'_0 = P_2 P_1 L_0 P^{-1}_1 P^{-1}_2$\n", "$L'_1 = P_2 L_1 P^{-1}_2$\n", "$L'_2 = L_2$\n", "\n", "Note also that the L' arrays are lower triangular, hence their inverse can be obtained trivially by taking the opposite of sub-diagonal elements of $L'$, as seen in class.\n", "\n", "Defining $L = {L'}^{-1}_2 {L'}^{-1}_1 {L'}^{-1}_0$ and $P = P_2 P_1 P_0$, the desired LU decomposition is then:\n", "\n", "$LU = PA$\n", "\n", "2) Based on the previous result, write an algorithm for LU decompostion with pivoting. First get the U array via standard Gauss elimination and record all the swaps in a suitable list. As you proceed, save also the L arrays (Without swapping). At the end, apply the right swaps to get the L' arrays and get the decomposition. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "from numpy.linalg import solve\n", "from numpy.linalg import inv\n", "\n", "def LUdec(A):\n", "\n", " # Given square array A, finds its LU decomposition\n", " # Applies partial pivoting\n", " \n", " N = len(A)\n", " \n", " # Vector of diagonal elements of A\n", " # To be updated during pivoting\n", " diag = np.empty(N)\n", " \n", " L = np.identity(N)\n", " \n", " swap = [ ]\n", " \n", " for i in range(N):\n", " \n", " # Partial pivoting\n", " ####################### \n", " # Find max element in column i\n", " ind = np.argmax(abs(A[:,i])) \n", " # If index of max is not the current i in the loop, swap rows\n", " if (ind > i):\n", " # Swap\n", " A[ [i,ind] ] = A[ [ind,i] ]\n", " # Recording all permutations in a list.\n", " # Needed to: permute rows of A, when comparing it with LU, \n", " # permute rows of L, permute entries of b when solving Ax = b\n", " swap.append(ind)\n", " else:\n", " swap.append(i)\n", " #########\n", " \n", " norm = A[i,i]\n", " L[i,i] = 1.0\n", " \n", " # Normalizing diagonal of A to 1 for gauss elimination\n", " #A[i,:] /= norm\n", " \n", " for j in range(i+1,N):\n", " # Updating L\n", " fact = A[j,i]\n", " L[j,i] = fact/norm\n", " # Gauss elimination to get U \n", " A[j,:] -= fact*(A[i,:]/norm)\n", " \n", " # Applying pivoting permutations to L\n", " # For each column j in L, I have to interchange rows \n", " # for all swaps that were performed at steps j+1...N\n", " # For example. If am considering column 1, I swap rows \n", " # loking at all the permutations from step 2 onwards, but I DON'T\n", " # apply swaps made at steps 0 and 1.\n", " for j in range(N-1):\n", " for i in range(j+1,N):\n", " p = swap[i]\n", " L[i,j], L[p,j] = L[p,j], L[i,j]\n", " \n", " return L, A, swap" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Example for testing. No pivoting needed here\n", "A = [ [1, 1, 4, 1], [1, 5, -1, -1], [1, -4, 2, 5], [2, -2, 1, 0]]\n", "A = np.array(A,dtype=float)\n", "\n", "#print(A,'\\n')\n", "\n", "Ain = np.copy(A)\n", "\n", "swaps = [ ]\n", "\n", "L, U, swaps = LUdec(A)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "disp = np.round(L,2)\n", "print('L=', '\\n', disp)\n", "disp = np.round(U,2)\n", "print('U=', '\\n', disp)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print('Checking that LU = A (modulo swaps)', '\\n')\n", "LU = np.dot(L,U)\n", "# Redoing the swaps in inverse order\n", "# to \"recompose\" the original A\n", "for i in range(len(swaps)-1,-1,-1):\n", " p = swaps[i]\n", " LU[[i,p]]=LU[[p,i]]\n", "disp = np.round(LU,2)\n", "print('LU = \\n',disp)\n", "print(' ')\n", "print('Input A = \\n', Ain)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Testing example with pivoting\n", "\n", "# Testing example for pivoting\n", "A = [ [0, 1, 4, 2], [3, 7, -1, 6], [1, -4, 1, 10], [0, 11, 0, 2]]\n", "A = np.array(A,dtype=float)\n", "#print(A)\n", "\n", "Ain = np.copy(A)\n", "\n", "L2, U2, swaps = LUdec(A)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "disp = np.round(L2,2)\n", "print('L=', '\\n', disp)\n", "disp = np.round(U2,2)\n", "print('U=', '\\n', disp)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print('Checking that LU = A', '\\n')\n", "LU = np.dot(L2,U2)\n", "# Redoing the swaps in inverse order\n", "# to \"recompose\" the original A\n", "for i in range(len(swaps)-1,-1,-1):\n", " p = swaps[i]\n", " LU[[i,p]]=LU[[p,i]]\n", "disp = np.round(LU,2)\n", "print(disp)\n", "print(' ')\n", "print(Ain)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Linear solver, using LU decomposition" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plan: use the LU function implemented above to initially compute and store A = LU. \n", "\n", "Then for any vector b, we can solve Ax = b as LUx = b via double backsubstitution. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "from numpy.linalg import solve\n", "from numpy.linalg import inv\n", "\n", "\n", "class solveinput:\n", " \n", " # Defining an object that contains all the input ingredients to solve Ax = b.\n", " # That is: L, U and the list which records the swaps when doing the LU decomp\n", " \n", " def __init__(self,L,U,rec):\n", " self.lower = L\n", " self.upper = U\n", " self.swaps = rec\n", "\n", "def LUdec(A):\n", "\n", " # Given square array A, finds its LU decomposition\n", " # Applies partial pivoting\n", " \n", " N = len(A)\n", " \n", " # Vector of diagonal elements of A\n", " # To be updated during pivoting\n", " diag = np.empty(N)\n", " \n", " L = np.identity(N)\n", " \n", " swap = [ ]\n", " \n", " for i in range(N):\n", " \n", " # Partial pivoting\n", " ####################### \n", " # Collecting diagonal elements of A in vector 'diag'\n", " # Find max element in diag \n", " ind = np.argmax(abs(A[:,i])) \n", " # If index of max is not the current i in the loop, swap rows\n", " if (ind > i):\n", " # Swap\n", " A[ [i,ind] ] = A[ [ind,i] ]\n", " # Recording all permutations in a list\n", " # Needed for later\n", " swap.append(ind)\n", " else:\n", " swap.append(i)\n", " #########\n", " \n", " norm = A[i,i]\n", " L[i,i] = 1.0\n", " \n", " for j in range(i+1,N):\n", " # Updating L\n", " fact = A[j,i]\n", " L[j,i] = fact/norm\n", " # Gauss elimination to get U \n", " A[j,:] -= fact*(A[i,:]/norm)\n", " \n", " # Applying pivoting permutations to L\n", " # For each column j in L, I have to interchange rows \n", " # for all swaps that were performed at steps j+1...N\n", " # For example. If am considering column 1, I swap rows \n", " # loking at all the permutations from step 2 onwards, but I DON'T\n", " # apply swaps made at steps 0 and 1.\n", " for j in range(N-1):\n", " for i in range(j+1,N):\n", " p = swap[i]\n", " L[i,j], L[p,j] = L[p,j], L[i,j]\n", " \n", " LU = solveinput(L,A,swap)\n", " \n", " return LU\n", "\n", "\n", "def linsolve_LU(LU,b):\n", " \n", " L = LU.lower\n", " U = LU.upper\n", " swaps = LU.swaps\n", " \n", " N = len(b)\n", " \n", " # Applying all swaps to b (can be done at once)\n", " for i in range(len(swaps)):\n", " b[i], b[swaps[i]] = b[swaps[i]], b[i]\n", " \n", " ######################################################\n", " # Now solving LUx = b via double backsubstitution\n", " # First solve Ly = b, to get y. Then solve Ux = y, \n", " # to get the desired solutionx \n", " \n", " # Solving Ly = b.\n", " # Initialize solution vector\n", " y = np.zeros(N)\n", " # Gets first solution\n", " y[0] = b[0]/L[0,0]\n", " # Proceed forward, line by line\n", " for i in range(1,N):\n", " v = np.dot(L[i,0:i],y[0:i])\n", " y[i] = (b[i] - v)/L[i,i]\n", " \n", " # Solving Ux = y\n", " # Initialize solution vector\n", " x = np.zeros(N) \n", " # Gets last solution, x_(N-1)\n", " x[N-1] = y[N-1]/U[N-1,N-1]\n", " # Proceed backward, line by line\n", " for i in range(N-2,-1,-1):\n", " v = np.dot(U[i,i+1:N],x[i+1:N])\n", " #if (i==1): print(i, i+1, N-1, A[i,i+1:N])\n", " x[i] = (y[i] - v)/U[i,i]\n", " \n", " return x" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Testing example\n", "A = [ [2, 1, 4, 1], [3, 4, -1, -1], [1, -4, 1, 5], [2, -2, 1, 3]]\n", "A = np.array(A,dtype=float)\n", "#print(A)\n", "\n", "b = [-4, 3, 9, 7]\n", "b = np.array(b,dtype=float)\n", "#print(b)\n", "\n", "LU = LUdec(A)\n", "\n", "#disp = np.round(U,2)\n", "#print('U =', '\\n', disp)\n", "#disp = np.round(L,2)\n", "#print('L=', '\\n', disp, '\\n')\n", "\n", "x = linsolve_LU(LU,b)\n", "\n", "print('My solution =', x)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Cross checking against numpy solver\n", "\n", "A = [ [2, 1, 4, 1], [3, 4, -1, -1], [1, -4, 1, 5], [2, -2, 1, 3]]\n", "A = np.array(A,dtype=float)\n", "#print(A)\n", "\n", "b = [-4, 3, 9, 7]\n", "b = np.array(b,dtype=float)\n", "#print(b)\n", "\n", "x = solve(A,b)\n", "print('Numpy solution = ', x)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Testing example for pivoting\n", "\n", "A = [ [0, 1, 4, 1], [3, 4, -1, -1], [1, -4, 1, 5], [2, -2, 1, 3]]\n", "A = np.array(A,dtype=float)\n", "#print(A)\n", "\n", "b = [-4, 3, 9, 7]\n", "b = np.array(b,dtype=float)\n", "#print(b)\n", "\n", "LU = LUdec(A)\n", "\n", "x = linsolve_LU(LU,b)\n", "\n", "print('My solution =', x)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "A = [ [0, 1, 4, 1], [3, 4, -1, -1], [1, -4, 1, 5], [2, -2, 1, 3]]\n", "A = np.array(A,dtype=float)\n", "#print(A)\n", "\n", "b = [-4, 3, 9, 7]\n", "b = np.array(b,dtype=float)\n", "\n", "y = solve(A,b)\n", "print('Numpy solution =', y)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Exploiting LU to solve systems with same A and different b\n", "\n", "A_in = np.array([ [0, 1, 4, 1, 0], [3, 4, -1, -1, 7], [1, -4, 1, 5, -14], [2, -2, 1, 3, -1], [21, -6, -9, 8, 1] ], dtype=float)\n", "b_in = np.array([0, 0, 3, -6, 21],float) \n", "\n", "A = np.copy(A_in)\n", "\n", "# Do this only once for any b\n", "LU = LUdec(A)\n", "\n", "###########################\n", "b = np.copy(b_in)\n", "x = linsolve_LU(LU,b)\n", "print('My solution =', x)\n", "\n", "A = np.copy(A_in)\n", "b = np.copy(b_in) \n", "y = solve(A,b)\n", "print('Numpy solution =', y, '\\n') \n", "\n", "b_in = np.array([1, 7, 102, -6, -49.5],float) \n", "b = np.copy(b_in)\n", "x = linsolve_LU(LU,b)\n", "print('My solution =', x)\n", "\n", "A = np.copy(A_in)\n", "b = np.copy(b_in) \n", "y = solve(A,b)\n", "print('Numpy solution =', y, '\\n') " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Gauss-Seidel" ] }, { "cell_type": "code", "execution_count": 108, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "\n", "def GSeid(A,b,x0,tol):\n", "\n", " # Solving linear systems via Gauss-Seidel method\n", " \n", " demomode = False\n", " if demomode:\n", " step = 1\n", " \n", " # Array size\n", " N = len(b)\n", " # Vector containing diagonal elements\n", " d = np.diag(A)\n", " # Vector containing inverse of diagonal\n", " dm1 = (np.diag(A))**(-1)\n", " \n", " x = x0\n", " \n", " acc = tol + 1.\n", " while (acc > tol):\n", " \n", " # \\sum Aij xj - Aii xi\n", " y = np.matmul(A,x) - np.multiply(x,d) \n", " # Gauss-Seidel formula\n", " x1 = np.multiply(dm1, b - y)\n", " \n", " if demomode:\n", " print('step =', step, 'x = ', x1)\n", " step += 1\n", " \n", " # Calculating accuracy and updating x\n", " acc = np.max(abs(x1 - x))\n", " x = x1 \n", " \n", " return x" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Example: Use Gauss-Seidel to solve problem from slide 44 in the \"linear algebra\" file." ] }, { "cell_type": "code", "execution_count": 109, "metadata": {}, "outputs": [], "source": [ "A = np.array([ [ 4, -1, 1], [-1, 4, -2], [1, -2, 4] ],dtype=float)\n", "b = [12.0, -1.0, 5.0]\n", "x0 = [1, 1, 1]\n", "tol = 1e-10" ] }, { "cell_type": "code", "execution_count": 110, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The solution is [3. 1. 1.]\n" ] } ], "source": [ "x = GSeid(A,b,x0,tol)\n", "\n", "print('The solution is ', x)" ] }, { "cell_type": "code", "execution_count": 121, "metadata": {}, "outputs": [], "source": [ "A = np.array([ [ 4, -1, 1, 3], [4, 10, -2, 1], [-2, 4, 14, -5], [5, -1, 4, 20]],dtype=float)\n", "b = [12.0, -1.0, 5.0, 2.0]\n", "x0 = [1, 1, 1, 1]\n", "tol = 1e-10" ] }, { "cell_type": "code", "execution_count": 122, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The solution is [ 3.20531654 -1.12836656 0.80902413 -0.91955229]\n" ] } ], "source": [ "x = GSeid(A,b,x0,tol)\n", "\n", "print('The solution is ', x)" ] }, { "cell_type": "code", "execution_count": 123, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[12. -1. 5. 2.] [12.0, -1.0, 5.0, 2.0]\n" ] } ], "source": [ "print(np.matmul(A,x), b)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "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.9.12" } }, "nbformat": 4, "nbformat_minor": 4 }