Skip to content

Instantly share code, notes, and snippets.

@DianSano
Last active March 8, 2022 12:18
Show Gist options
  • Select an option

  • Save DianSano/752ab300eaea7bde88dd793004432e86 to your computer and use it in GitHub Desktop.

Select an option

Save DianSano/752ab300eaea7bde88dd793004432e86 to your computer and use it in GitHub Desktop.
Lab_SC02_2022_Gauss-Seidel Method.ipynb
Display the source blob
Display the rendered blob
Raw
{
"nbformat": 4,
"nbformat_minor": 0,
"metadata": {
"colab": {
"name": "Lab_SC02_2022_Gauss-Seidel Method.ipynb",
"provenance": [],
"authorship_tag": "ABX9TyPvRGUNtvbLK6w5LMAWg94+",
"include_colab_link": true
},
"kernelspec": {
"name": "python3",
"display_name": "Python 3"
},
"language_info": {
"name": "python"
}
},
"cells": [
{
"cell_type": "markdown",
"metadata": {
"id": "view-in-github",
"colab_type": "text"
},
"source": [
"<a href=\"https://colab.research.google.com/gist/DianSano/752ab300eaea7bde88dd793004432e86/lab_sc02_2022_gauss-seidel-method.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"id": "E-mbzJVTEtT1"
},
"outputs": [],
"source": [
"import numpy as np"
]
},
{
"cell_type": "markdown",
"source": [
"**# Let us first check if the coefficient matrix is diagonally dominant or not.**"
],
"metadata": {
"id": "pXcM5KHMOIwO"
}
},
{
"cell_type": "code",
"source": [
"a = [[8, 3, -3], [-2, -8, 5], [3, 5, 10]]\n",
"\n",
"# Find diagonal coefficients\n",
"diag = np.diag(np.abs(a)) \n",
"\n",
"# Find row sum without diagonal\n",
"off_diag = np.sum(np.abs(a), axis=1) - diag \n",
"\n",
"if np.all(diag > off_diag):\n",
" print('matrix is diagonally dominant')\n",
"else:\n",
" print('NOT diagonally dominant')"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "FbZ7rmp-FGkX",
"outputId": "09d07d75-5524-4751-d39d-421c77199465"
},
"execution_count": null,
"outputs": [
{
"output_type": "stream",
"name": "stdout",
"text": [
"matrix is diagonally dominant\n"
]
}
]
},
{
"cell_type": "markdown",
"source": [
"**Since it is guaranteed to converge, we can use Gauss-Seidel method to solve it.**"
],
"metadata": {
"id": "jDnw1k0sOW6P"
}
},
{
"cell_type": "code",
"source": [
"x1 = 0\n",
"x2 = 0\n",
"x3 = 0\n",
"epsilon = 0.01\n",
"converged = False\n",
"\n",
"x_old = np.array([x1, x2, x3])\n",
"\n",
"print('Iteration results')\n",
"print(' k, x1, x2, x3 ')\n",
"for k in range(1, 50):\n",
" x1 = (14-3*x2+3*x3)/8\n",
" x2 = (5+2*x1-5*x3)/(-8)\n",
" x3 = (-8-3*x1-5*x2)/(-5)\n",
" x = np.array([x1, x2, x3])\n",
" # check if it is smaller than threshold\n",
" dx = np.sqrt(np.dot(x-x_old, x-x_old))\n",
" \n",
" print(\"%d, %.4f, %.4f, %.4f\"%(k, x1, x2, x3))\n",
" if dx < epsilon:\n",
" converged = True\n",
" print('Converged!')\n",
" break\n",
" \n",
" # assign the latest x value to the old value\n",
" x_old = x\n",
"\n",
"if not converged:\n",
" print('Not converge, increase the # of iterations')"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "-SHj1lGHN75v",
"outputId": "b1a99f11-3c30-4112-cf3c-071ce720d2c8"
},
"execution_count": null,
"outputs": [
{
"output_type": "stream",
"name": "stdout",
"text": [
"Iteration results\n",
" k, x1, x2, x3 \n",
"1, 1.7500, -1.0625, 1.5875\n",
"2, 2.7437, -0.3188, 2.9275\n",
"3, 2.9673, 0.4629, 3.8433\n",
"4, 3.0177, 1.0226, 4.4332\n",
"5, 3.0290, 1.3885, 4.8059\n",
"6, 3.0315, 1.6208, 5.0397\n",
"7, 3.0321, 1.7668, 5.1861\n",
"8, 3.0322, 1.8582, 5.2776\n",
"9, 3.0322, 1.9154, 5.3348\n",
"10, 3.0323, 1.9512, 5.3705\n",
"11, 3.0323, 1.9735, 5.3929\n",
"12, 3.0323, 1.9875, 5.4068\n",
"13, 3.0323, 1.9962, 5.4156\n",
"14, 3.0323, 2.0017, 5.4210\n",
"Converged!\n"
]
}
]
},
{
"cell_type": "markdown",
"source": [
"# Solve Systems of Linear Equations in Python¶\n",
"Though we discussed various methods to solve the systems of linear equations, it is actually very easy to do it in Python. In this section, we will use Python to solve the systems of equations. The easiest way to get a solution is via the solve function in Numpy."
],
"metadata": {
"id": "nTHBodm5WmyU"
}
},
{
"cell_type": "markdown",
"source": [
"## TRY IT! Use numpy.linalg.solve to solve the following equations.\n",
"\n",
"# 4𝑥1+3𝑥2−5𝑥3=2\n",
"# -2𝑥1−4𝑥2+5𝑥3=5\n",
"# 8𝑥1+8𝑥2=−3"
],
"metadata": {
"id": "_lm-xhw1Wmdm"
}
},
{
"cell_type": "code",
"source": [
"import numpy as np\n",
"\n",
"A = np.array([[4, 3, -5], \n",
" [-2, -4, 5], \n",
" [8, 8, 0]])\n",
"y = np.array([2, 5, -3])\n",
"\n",
"x = np.linalg.solve(A, y)\n",
"print(x)"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "014vihbbW3Kq",
"outputId": "e378f538-0153-46e0-a651-788c858b12d1"
},
"execution_count": 1,
"outputs": [
{
"output_type": "stream",
"name": "stdout",
"text": [
"[ 2.20833333 -2.58333333 -0.18333333]\n"
]
}
]
},
{
"cell_type": "markdown",
"source": [
"We can see we get the same results as that in the previous section when we calculated by hand. Under the hood, the solver is actually doing a LU decomposition to get the results. You can check the help of the function, it needs the input matrix to be square and of full-rank, i.e., all rows (or, equivalently, columns) must be linearly independent."
],
"metadata": {
"id": "44ftAu33XSY9"
}
},
{
"cell_type": "markdown",
"source": [
"TRY IT! Try to solve the above equations using the matrix inversion approach.\n",
"\n"
],
"metadata": {
"id": "q1NWAWe2XTaS"
}
},
{
"cell_type": "code",
"source": [
"A_inv = np.linalg.inv(A)\n",
"\n",
"x = np.dot(A_inv, y)\n",
"print(x)"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "zTKDCN3gXYi-",
"outputId": "3d9e96b4-fc4a-4746-d526-d2adf527212f"
},
"execution_count": 2,
"outputs": [
{
"output_type": "stream",
"name": "stdout",
"text": [
"[ 2.20833333 -2.58333333 -0.18333333]\n"
]
}
]
}
]
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment