{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"id": "0a394fca-3535-4541-80aa-dc0953aee499",
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"import pandas as pd"
]
},
{
"cell_type": "markdown",
"id": "7083ba64-eda7-4ab7-a388-4ed048cc47bf",
"metadata": {},
"source": [
"# Exercise sheet 7"
]
},
{
"cell_type": "markdown",
"id": "7686a72f-ecd2-469b-9cc6-2754d8cda4ba",
"metadata": {},
"source": [
"## Free training"
]
},
{
"cell_type": "markdown",
"id": "8377cf5d-cfbb-417f-bc0c-dc25d5b31d82",
"metadata": {},
"source": [
"### Naive Gaussian Elemination method:"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "59b454e3-1d11-4935-b3da-2691782fe42d",
"metadata": {},
"outputs": [],
"source": [
"# first we define a function that switches all the rows of a matrix in a cyclic permutation:\n",
"\n",
"def switch (B): # B is the concatenated A and b\n",
" \n",
" dim = len(B)\n",
" \n",
" B1 = np.zeros((dim, dim+1))\n",
" \n",
" B[dim-1] = B[0]\n",
" \n",
" for i in range(dim-1):\n",
" B[i] = B[i+1]\n",
" \n",
" return A, b"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "955114df-d784-4f41-b571-9e99878d57e3",
"metadata": {},
"outputs": [],
"source": [
"# define the method which takes a NxN-matrix A and a vector b\n",
"# and returns the upper triangular matrix A1 and the corresponding vector b1\n",
"# this method is implemented so that also division by zero cannot happen\n",
"\n",
"\n",
"def triangular (A, b): # insert A and b as numpy-arrays (float!!)\n",
" \n",
" dim = len(A)\n",
" \n",
" B = np.concatenate((A.T, [b])).T # use the transformation also for b\n",
" \n",
" L = np.zeros((dim, dim))\n",
" \n",
" # we want L to have 1s on the diagonal:\n",
" \n",
" for i in range(dim):\n",
" L[i][i] = 1\n",
" \n",
" # the rest of L will be filled later!\n",
" \n",
" i = 0\n",
" \n",
" while (i < dim): # this loops over the elements that get eliminated\n",
" \n",
" # np.any(A) returns True if it exists one non-zero element in the last (dim-i'th) elements of the i'th column \n",
"\n",
" if (np.any(A[i:], 0)[i]):\n",
"\n",
" # If there are non-zero elements in the i'th column but the i'th element is zero, switch the rows until\n",
" # the first element is not zero anymore\n",
"\n",
" while (B[i][i] == 0):\n",
" B = switch(B)\n",
"\n",
" # then the first element is non-zero --> eliminate it in all rows:\n",
"\n",
" j = i+1\n",
"\n",
"\n",
" while (j < dim): # this loops over the different rows\n",
"\n",
" # write A and b in an array so that the multiplication is done at the same time\n",
" # with this we avoid rewriting the variabels of A for different steps\n",
" \n",
" \n",
" L[j][i] = B[j][i]/B[i][i] # multipliers\n",
" \n",
" B[j] = B[j] - L[j][i] * B[i]\n",
" \n",
" \n",
" j += 1\n",
" \n",
" i += 1\n",
"\n",
"\n",
" # if all the elements in the first column are zero, go ahead with the next step:\n",
"\n",
" else:\n",
" \n",
" i += 1\n",
" \n",
" A1 = B.T[:-1].T\n",
" b1 = B.T[-1].T\n",
" \n",
" return A1, b1, L"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "ef42a2a2-fc22-4ddc-964f-b7f4870b29e9",
"metadata": {},
"outputs": [],
"source": [
"# we define the backward substitution:\n",
"\n",
"def backwards (A, b): # takes a upper triangular matrix A and its corresponding vector b\n",
" \n",
" dim = len(A)\n",
" \n",
" x = np.zeros(dim)\n",
" \n",
" # the last element is already known:\n",
" \n",
" x[dim-1] = b[dim-1]/A[dim-1][dim-1]\n",
" \n",
" # calculate all the other elements:\n",
" \n",
" \n",
" \n",
" for i in range (dim-2, -1, -1): # go from the second last element up to the zeroth element\n",
" \n",
" #calculate the sum needed for calculation of x_i:\n",
" \n",
" sum0 = 0\n",
" \n",
" for j in range(i+1, dim):\n",
" \n",
" sum0 += A[i][j] * x[j]\n",
" \n",
" x[i] = (b[i] - sum0)/A[i][i]\n",
" \n",
" return x\n",
"\n",
" "
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "6dcad87f-cb7f-4f26-ac31-fb9b82f8e543",
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"id": "3c9c0f90-4280-4091-89cd-1b4cdb8b1d11",
"metadata": {},
"source": [
"### Naive Gaussian Elemination with Partial Pivoting:"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "1c280044-2733-4914-9592-46ce8b0857d4",
"metadata": {},
"outputs": [],
"source": [
"# we do the same as before, but now always switch the rows so that the biggest absolute value in the k'th column is at the top\n",
"# therefore we define a new switch function:\n",
"\n",
"def switch2 (B, a, b): # this function switches the two rows a and b in matrix B\n",
" \n",
" dim1 = len(B[0]) # number if columns\n",
" dim2 = len(B.T[0]) # number of rows\n",
" \n",
" B1 = np.zeros((dim2, dim1))\n",
" \n",
" for i in range(dim2):\n",
" if (i == a):\n",
" B1[i] = B[b]\n",
" elif (i == b):\n",
" B1[i] = B[a]\n",
" else: \n",
" B1[i] = B[i]\n",
" \n",
" return B1\n",
" \n",
"\n",
"def switch_max (B, k): # input: the concatenated B = (A, b) and the number of the current step k, hence we consider the biggest absolute value in the k'th column\n",
" # and change that two rows\n",
" \n",
" dim = len(B) # number of rows\n",
"\n",
" column = B.T[k]\n",
" \n",
" maximum = np.max(np.abs(column)) # maximum of the absolute values in the given column\n",
" \n",
"\n",
" for i in range(k, dim):\n",
" if (maximum == np.abs(column[i])): # if the element of the column is the maximum, switch this row:\n",
" B = switch2(B, k, i)\n",
" \n",
" return B\n",
"\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "2f071c71-9028-44ac-8bb0-99db14157562",
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": 6,
"id": "16b94d52-b5e1-4cc6-a7b4-b9fbd6a89714",
"metadata": {},
"outputs": [],
"source": [
"# define the method which takes a NxN-matrix A and a vector b\n",
"# and returns the upper triangular matrix A1 and the corresponding vector b1\n",
"# this method is implemented so that also division by zero cannot happen\n",
"# and we always shift the rows so that the biggest absolute value is at the top\n",
"\n",
"# note that we already implement some parts of the LU decomposition:\n",
"# we save the multipliers in an matrix L and also return this\n",
"\n",
"\n",
"def triangular2 (A, b): # insert A and b as numpy-arrays (float!!)\n",
" \n",
" dim = len(A)\n",
" \n",
" B = np.concatenate((A.T, [b])).T # use the transformation also for b, easiest if one concatenates A and b into a new matrix\n",
" \n",
" L = np.zeros((dim, dim))\n",
" \n",
" # we want L to have 1s on the diagonal:\n",
" \n",
" for i in range(dim):\n",
" L[i][i] = 1\n",
" \n",
" # the rest of L will be filled later!\n",
" \n",
" i = 0\n",
" \n",
" while (i < dim): # this loops over the elements that get eliminated\n",
" \n",
" # np.any(A) returns True if it exists one non-zero element in the last (dim-i'th) elements of the i'th column \n",
"\n",
" if (np.any(A[i:], 0)[i]):\n",
" \n",
"\n",
" # If there are non-zero elements in the first column but the first is zero, switch the biggest element at the top!!\n",
"\n",
" B = switch_max(B, i)\n",
"\n",
" # then the first element is non-zero --> eliminate it in all rows:\n",
"\n",
" j = i+1\n",
"\n",
"\n",
" while (j < dim): # this loops over the different rows\n",
"\n",
" # write A and b in an array so that the multiplication is done at the same time\n",
" # with this we avoid rewriting the variabels of A for different steps\n",
" \n",
" #################### Cause for inaccuracies suspected here, see Nr. 1:\n",
" \n",
" L[j][i] = B[j][i]/B[i][i] # multipliers\n",
" \n",
" B[j] = B[j] - B[j][i]/B[i][i] * B[i]\n",
" \n",
" ####################\n",
" \n",
" \n",
" j += 1\n",
" \n",
" i += 1\n",
"\n",
"\n",
" # if all the elements in the first column are zero, go ahead with the next step:\n",
"\n",
" else:\n",
" \n",
" i += 1\n",
" \n",
" A1 = B.T[:-1].T\n",
" b1 = B.T[-1].T\n",
" \n",
" return A1, b1, L"
]
},
{
"cell_type": "markdown",
"id": "7b339fd7-e0e8-46a8-9ad3-d686072607db",
"metadata": {},
"source": [
"### LU decomposition:"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "94a15008-fa85-45e7-a8d9-a3535d22ae0b",
"metadata": {},
"outputs": [],
"source": [
"# this will be done in the exercise 2"
]
},
{
"cell_type": "markdown",
"id": "e0b9eed9-20aa-4001-933d-2d13cbbf3a6b",
"metadata": {
"tags": []
},
"source": [
"### Exercise 1:"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "deeb9966-53fb-4cb1-9e39-8bb7bdd36963",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"\n",
"Naive Gaussian Elimination method:\n",
"\n",
"The upper triangular matrix and its corresponding vector is:\n",
"\n",
"A:\n",
"[[ 13. 4. 7. 9. ]\n",
" [ 0. 2.92307692 -0.38461538 5.07692308]\n",
" [ 0. 0. 2.47368421 1.94736842]\n",
" [ 0. 0. 0. -25.68085106]]\n",
"\n",
"b:\n",
"[ 111. 32.61538462 19.63157895 -128.40425532]\n",
"\n",
"The solution of the given system is:\n",
"\n",
"x = [2. 3. 4. 5.]\n",
"\n",
"\n",
"Naive Gaussian Elimination with Partial Pivoting:\n",
"\n",
"The upper triangular matrix and its corresponding vector is:\n",
"\n",
"A:\n",
"[[ 1.30000000e+01 4.00000000e+00 7.00000000e+00 9.00000000e+00]\n",
" [ 0.00000000e+00 1.30769231e+01 1.33846154e+01 8.92307692e+00]\n",
" [ 0.00000000e+00 0.00000000e+00 -6.41176471e+00 1.00588235e+01]\n",
" [ 0.00000000e+00 0.00000000e+00 4.44089210e-16 -2.21467890e+00]]\n",
"\n",
"b:\n",
"[111. 137.38461538 24.64705882 -11.0733945 ]\n",
"\n",
"The solution of the given system is:\n",
"\n",
"x = [2. 3. 4. 5.]\n"
]
}
],
"source": [
"A = np.double(np.array([[13, 4, 7, 9], [10, 6, 5, 12], [1, 8, 2, 16], [3, 14, 15, 11]]))\n",
"b = np.double(np.array([111, 118, 114, 163]))\n",
"\n",
"print(f'''\n",
"\n",
"Naive Gaussian Elimination method:\n",
"\n",
"The upper triangular matrix and its corresponding vector is:\n",
"\n",
"A:\n",
"{triangular(A, b)[0]}\n",
"\n",
"b:\n",
"{triangular(A, b)[1]}\n",
"\n",
"The solution of the given system is:\n",
"\n",
"x = {backwards(triangular(A, b)[0], triangular(A, b)[1])}''')\n",
"\n",
"A = np.double(np.array([[13, 4, 7, 9], [10, 6, 5, 12], [1, 8, 2, 16], [3, 14, 15, 11]]))\n",
"b = np.double(np.array([111, 118, 114, 163]))\n",
"\n",
"print(f'''\n",
"\n",
"Naive Gaussian Elimination with Partial Pivoting:\n",
"\n",
"The upper triangular matrix and its corresponding vector is:\n",
"\n",
"A:\n",
"{triangular2(A, b)[0]}\n",
"\n",
"b:\n",
"{triangular2(A, b)[1]}\n",
"\n",
"The solution of the given system is:\n",
"\n",
"x = {backwards(triangular2(A, b)[0], triangular2(A, b)[1])}''')"
]
},
{
"cell_type": "markdown",
"id": "a7f71a03-c5ab-4309-af6d-3c786356a4e4",
"metadata": {},
"source": [
"Comparing both algorithms one can see that we indeed obtain the same solution for x. However, for the triangular matrix, the errors seem to be greater when calculated by the pivot method than when calculated via the naive method, while the opposite would be expected. This might be due to the fact that the forward elimination of the pivoting method is carried out in the exact same manner as in the naive method, which does counterproductively guarantee that the values are divided by a maximum value, causing very small results. This could be fixed by changing up the order of operations a bit in the block of code in \"triangular2\" highlighted with \"#####...\" above, so that the multiplications be carried out before the divisions for example. This was tried, but did not yield better results, it remains our best guess for the source of the error however."
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "55a88b8a-4a6c-4b98-a5ce-b09d59e084a5",
"metadata": {},
"outputs": [],
"source": [
"# with this, we can easily determine the determinante, since the determinante of an upper triangular matrix\n",
"# is just the product of its diagonal elements:\n",
"\n",
"def det(A):\n",
" \n",
" dim = len(A)\n",
" \n",
" b = np.zeros(dim) # arbitrary, just because triangular needs this input\n",
" \n",
" # note that a cyclic permutation of the rows as done in 'switch' for triangular does not effect the determinante\n",
" # for triangular2, we use switch_max, this can effect the sign of the determinante\n",
" # here we use triangular to avoid this problem\n",
" \n",
" A1 = triangular(A, b)[0]\n",
" \n",
" det = 1\n",
" \n",
" for i in range(dim):\n",
" \n",
" det *= A1[i][i]\n",
" \n",
" return det\n",
" "
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "f3114264-b96a-4f77-b070-5d821b640e85",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"The determinante of the given matrix is:\n",
"\n",
"-2414.000000000001\n"
]
}
],
"source": [
"print(f'''\n",
"The determinante of the given matrix is:\n",
"\n",
"{det(A)}''')\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "1f81505a-2de6-4587-b7e1-c35be2912d83",
"metadata": {},
"source": [
"## Exercise 2"
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "d515f0bd-0420-4c68-8448-01068a28fb66",
"metadata": {},
"outputs": [],
"source": [
"# LU decomposition:\n",
"\n",
"# we write a function that decompose a squared matrix A into an upper triangular matrix U and a lower one L:\n",
"def decomp (A):\n",
" \n",
" dim = len(A)\n",
" \n",
" b = np.zeros(dim) # arbitrary vector since not needed for this function\n",
" \n",
" # the matrix U is just the same as after the Gaussian forward elimination:\n",
" \n",
" U = triangular(A, b)[0]\n",
" \n",
" # the matrix L consists of the multipliers and is already calculated within the triangular function:\n",
" \n",
" L = triangular(A, b)[2]\n",
" \n",
" return L, U\n",
"\n",
"# in order to check wether we are right, we define a matrix multiplication for squared matrices:\n",
"\n",
"def product(A, B): # input: two NxN-matrices as numpy arrays!\n",
" \n",
" dim = len(A) # dimension N of NxN-matrices\n",
" \n",
" AB = np.zeros((dim, dim))\n",
" \n",
" for i in range (dim):\n",
" for j in range (dim):\n",
" AB[i][j] = sum(A[i]*B.T[j]) # common matrix multiplication\n",
" return AB"
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "ad96ba85-a88f-40d1-8a5e-6b4f72d8331c",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"\n",
"Decomposition of the given matrix A:\n",
"\n",
"A = L*U with:\n",
"\n",
"L:\n",
"[[1. 0. 0. 0. ]\n",
" [0.76923077 1. 0. 0. ]\n",
" [0.07692308 2.63157895 1. 0. ]\n",
" [0.23076923 4.47368421 6.10638298 1. ]]\n",
"\n",
"U:\n",
"[[ 13. 4. 7. 9. ]\n",
" [ 0. 2.92307692 -0.38461538 5.07692308]\n",
" [ 0. 0. 2.47368421 1.94736842]\n",
" [ 0. 0. 0. -25.68085106]]\n",
"\n",
"test: the product L*U gives:\n",
"\n",
"[[13. 4. 7. 9.]\n",
" [10. 6. 5. 12.]\n",
" [ 1. 8. 2. 16.]\n",
" [ 3. 14. 15. 11.]]\n",
"\n",
"This is exactly the matrix A, hence our algorithm works!\n"
]
}
],
"source": [
"# decomposition of A:\n",
"\n",
"A = np.double(np.array([[13, 4, 7, 9], [10, 6, 5, 12], [1, 8, 2, 16], [3, 14, 15, 11]]))\n",
"b = np.double(np.array([111, 118, 114, 163]))\n",
"\n",
"L = decomp(A)[0]\n",
"U = decomp(A)[1]\n",
"product1 = product(decomp(A)[0], decomp(A)[1]) # should be A\n",
"\n",
"print(f'''\n",
"\n",
"Decomposition of the given matrix A:\n",
"\n",
"A = L*U with:\n",
"\n",
"L:\n",
"{L}\n",
"\n",
"U:\n",
"{U}\n",
"\n",
"test: the product L*U gives:\n",
"\n",
"{product1}\n",
"\n",
"This is exactly the matrix A, hence our algorithm works!''')"
]
},
{
"cell_type": "code",
"execution_count": 13,
"id": "ced232fa-512a-4a2c-acec-e407d0277741",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"The solution of the given linear equation computed with LU decomposition is:\n",
"\n",
"X = [2. 3. 4. 5.]\n",
"\n",
"This gives the same as before, hence our algorithm worked again.\n",
"\n"
]
}
],
"source": [
"# we solve the linear equation with LU decomposition:\n",
"\n",
"A = np.double(np.array([[13, 4, 7, 9], [10, 6, 5, 12], [1, 8, 2, 16], [3, 14, 15, 11]]))\n",
"b = np.double(np.array([111, 118, 114, 163]))\n",
"\n",
"def sol_LU(A, b):\n",
"\n",
" # first we have to compute Z out of L*Z = C with C = b:\n",
" \n",
" L = decomp(A)[0]\n",
" U = decomp(A)[1]\n",
"\n",
" C = b\n",
"\n",
" Z = backwards(triangular(L, C)[0], triangular(L, C)[1])\n",
"\n",
" # then we have to compute X out of U*X = Z:\n",
"\n",
" X = backwards(triangular(U, Z)[0], triangular(U, Z)[1])\n",
" \n",
" return X\n",
"\n",
"print(f'''\n",
"The solution of the given linear equation computed with LU decomposition is:\n",
"\n",
"X = {sol_LU(A, b)}\n",
"\n",
"This gives the same as before, hence our algorithm worked again.\n",
"''')"
]
},
{
"cell_type": "code",
"execution_count": 14,
"id": "f1c8cccb-a944-4c05-98ae-aba1fba39531",
"metadata": {},
"outputs": [],
"source": [
"# we write a function for the inverse of a square matrix. Let's call the inverse matrix J and the unitary matrix I:\n",
"\n",
"# we decompose the equation A*J = I in N linear equations of A*J_i = I_i\n",
"# where J_i is the i'th column of J and I_e the i'th column of I.\n",
"# These equations can be solved for J_i, hence we obtain J\n",
"# the solution of the N different linear equations is done with the LU decomposition\n",
"\n",
"\n",
"def inverse(A):\n",
"\n",
" dim = len(A) # dim = N = number of rows/columns\n",
" \n",
" J = np.zeros((dim, dim)) # initialize the inverse to fill it later\n",
" \n",
" # initialize unitary matrix and fill the diagonal elements with ones\n",
" \n",
" I = np.zeros((dim, dim)) \n",
" \n",
" for i in range (dim):\n",
" I[i][i] = 1\n",
"\n",
" # solve the N different equations:\n",
" \n",
" for i in range (dim):\n",
" J.T[i] = sol_LU(A, I[i])\n",
" \n",
" return J"
]
},
{
"cell_type": "code",
"execution_count": 15,
"id": "6b72f686-2cfa-4d8c-8a3a-6aeae69c481d",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"The inverse of the matrix A is given by:\n",
"\n",
"[[-0.1971831 0.37613919 -0.16321458 -0.01159901]\n",
" [-0.78169014 1.03314002 -0.38442419 0.07166529]\n",
" [ 0.52112676 -0.70836785 0.21706711 0.03065452]\n",
" [ 0.33802817 -0.45153273 0.23777962 -0.03893952]]\n",
"\n",
"We check this again just by multiplication and we expect to obtain the unitary matrix again:\n",
"\n",
"A*J:\n",
"\n",
"[[ 1.00000000e+00 0.00000000e+00 0.00000000e+00 5.55111512e-17]\n",
" [ 0.00000000e+00 1.00000000e+00 8.88178420e-16 5.55111512e-17]\n",
" [ 0.00000000e+00 -8.88178420e-16 1.00000000e+00 0.00000000e+00]\n",
" [-4.44089210e-16 8.88178420e-16 -8.88178420e-16 1.00000000e+00]]\n",
"\n",
"Rounded we obtain:\n",
"\n",
"[[1. 0. 0. 0.]\n",
" [0. 1. 0. 0.]\n",
" [0. 0. 1. 0.]\n",
" [0. 0. 0. 1.]]\n",
"\n"
]
}
],
"source": [
"# compute the Inverse of A:\n",
"\n",
"A = np.double(np.array([[13, 4, 7, 9], [10, 6, 5, 12], [1, 8, 2, 16], [3, 14, 15, 11]]))\n",
"\n",
"J = inverse(A)\n",
"\n",
"product2 = product(A, J)\n",
"\n",
"product2_rounded = np.abs(np.round(product2))\n",
"\n",
"print(f'''\n",
"The inverse of the matrix A is given by:\n",
"\n",
"{J}\n",
"\n",
"We check this again just by multiplication and we expect to obtain the unitary matrix again:\n",
"\n",
"A*J:\n",
"\n",
"{product2}\n",
"\n",
"Rounded we obtain:\n",
"\n",
"{product2_rounded}\n",
"''')"
]
},
{
"cell_type": "markdown",
"id": "f78ede95-2a2f-47e7-b833-9d51dae3fcd8",
"metadata": {},
"source": [
"It can be seen that in the limits of truncation errors we indeed obtain the unitary matrix which means that J is really the inverse of A."
]
}
],
"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.7.3"
}
},
"nbformat": 4,
"nbformat_minor": 5
}