system-of-equations.ipynb

{
 "metadata": {
  "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.5-final"
  },
  "orig_nbformat": 2,
  "kernelspec": {
   "name": "venv",
   "display_name": "venv"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2,
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# System of Linear equations\n",
    "\n",
    "Iterative techniques to find solutions of `AX = b`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Jacobi method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "def jacobi(dim, A, B, tol, n):\n",
    "    \"\"\"\n",
    "    Arguments\n",
    "    ---------\n",
    "    dim: int\n",
    "        Dimension of square matrix `A`\n",
    "    A: 2D np.array\n",
    "        Coefficient matrix\n",
    "    B: np.array\n",
    "        Constant matrix\n",
    "    tol: float\n",
    "        tolerance\n",
    "    n: int\n",
    "        maximum number of iterations\n",
    "    \"\"\"\n",
    "    X0 = np.array(dim*[np.inf])\n",
    "    X1 = np.array(dim*[0.]) # Initial vector\n",
    "\n",
    "    iterations = 0\n",
    "    while np.max(abs(X0 - X1)) > tol and iterations < n:\n",
    "        X0 = X1\n",
    "        X1 = np.array(dim*[0.])\n",
    "\n",
    "        for i in range(dim):\n",
    "            for j in range(dim):\n",
    "                if j != i:\n",
    "                    X1[i] += (-A[i][j]*X0[j])\n",
    "            X1[i] += B[i]\n",
    "            X1[i] /= A[i][i]\n",
    "        \n",
    "        print(X1)\n",
    "        iterations += 1\n",
    "\n",
    "    print(\"iterations:\", iterations)\n",
    "    \n",
    "    return X1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "dim = 3\n",
    "A = np.array([\n",
    "    [15, 3, -2],\n",
    "    [2, 10, 1],\n",
    "    [1, -2, 8]\n",
    "    ], dtype=float)\n",
    "B = np.array([85, 51, 5], dtype=float)\n",
    "tol = 1e-5\n",
    "n = 100\n",
    "\n",
    "jacobi(dim, A, B, tol, n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gauss-Seidel method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "def gauss_seidel(dim, A, B, tol, n):\n",
    "    \"\"\"\n",
    "    Arguments\n",
    "    ---------\n",
    "    dim: int\n",
    "        Dimension of square matrix `A`\n",
    "    A: 2D np.array\n",
    "        Coefficient matrix\n",
    "    B: np.array\n",
    "        Constant matrix\n",
    "    tol: float\n",
    "        tolerance\n",
    "    n: int\n",
    "        maximum number of iterations\n",
    "    \"\"\"\n",
    "    X0 = np.array(dim*[np.inf])\n",
    "    X1 = np.array(dim*[0.]) # Initial vector\n",
    "\n",
    "    iterations = 0\n",
    "    while np.max(abs(X0 - X1)) > tol and iterations < n:\n",
    "        X0 = X1\n",
    "        X1 = np.array(dim*[0.])\n",
    "\n",
    "        for i in range(dim):\n",
    "            for j in range(i):\n",
    "                X1[i] += (-A[i][j]*X1[j])\n",
    "            for j in range(i+1, dim):\n",
    "                X1[i] += (-A[i][j]*X0[j])\n",
    "            \n",
    "            X1[i] += B[i]\n",
    "            X1[i] /= A[i][i]\n",
    "        \n",
    "        print(X1)\n",
    "        iterations += 1\n",
    "\n",
    "    print(\"iterations:\", iterations)\n",
    "    \n",
    "    return X1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "dim = 3\n",
    "A = np.array([\n",
    "    [8, 2, -2],\n",
    "    [1, -8, 3],\n",
    "    [2, 1, 9]\n",
    "    ], dtype=float)\n",
    "B = np.array([8, -4, 12], dtype=float)\n",
    "tol = 1e-5\n",
    "n = 100\n",
    "\n",
    "gauss_seidel(dim, A, B, tol, n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## SOR (Succesive Over Reduction) method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "def SOR(dim, A, B, ω, tol, n):\n",
    "    \"\"\"\n",
    "    Arguments\n",
    "    ---------\n",
    "    dim: int\n",
    "        Dimension of square matrix `A`\n",
    "    A: 2D np.array\n",
    "        Coefficient matrix\n",
    "    B: np.array\n",
    "        Constant matrix\n",
    "    ω: float\n",
    "        Weigth to residual part of Gauss-Siedel method\n",
    "    tol: float\n",
    "        tolerance\n",
    "    n: int\n",
    "        maximum number of iterations\n",
    "    \"\"\"\n",
    "    X0 = np.array(dim*[np.inf])\n",
    "    X1 = np.array(dim*[0.]) # Initial vector\n",
    "\n",
    "    iterations = 0\n",
    "    while np.max(abs(X0 - X1)) > tol and iterations < n:\n",
    "        X0 = X1\n",
    "        X1 = np.array(dim*[0.])\n",
    "\n",
    "        for i in range(dim):\n",
    "            for j in range(i):\n",
    "                X1[i] += (-A[i][j]*X1[j])\n",
    "            for j in range(i+1, dim):\n",
    "                X1[i] += (-A[i][j]*X0[j])\n",
    "            X1[i] += B[i]\n",
    "            X1[i] /= A[i][i]\n",
    "\n",
    "            X1[i] *= ω\n",
    "            X1[i] += (1-ω)*X0[i]\n",
    "        \n",
    "        print(X1)\n",
    "        iterations += 1\n",
    "\n",
    "    print(\"iterations:\", iterations)\n",
    "    \n",
    "    return X1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "dim = 3\n",
    "A = np.array([\n",
    "    [4, 3, 0],\n",
    "    [3, 4, -1],\n",
    "    [0, -1, 4]\n",
    "    ], dtype=float)\n",
    "B = np.array([24, 30, -24], dtype=float)\n",
    "ω = 1.25\n",
    "tol = 1e-5\n",
    "n = 100\n",
    "\n",
    "SOR(dim, A, B, ω, tol, n)"
   ]
  },
  {
   "source": [
    "# System of non-linear equations\n",
    "\n",
    "Atleast one equation is non-linear."
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "source": [
    "## Newton-Raphson method"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "def newton_raphson(F, J, X, tol, n):\n",
    "    \"\"\"\n",
    "    Arguments\n",
    "    ---------\n",
    "    F: function\n",
    "        Function returning vector of functions\n",
    "    J: function returning a matrix\n",
    "        Function returning Jacobian matrix\n",
    "    X: np.array\n",
    "        Initial approximation vector\n",
    "    tol: float\n",
    "        tolerance\n",
    "    n: int\n",
    "        maximum number of iterations\n",
    "    \"\"\"\n",
    "    E = np.array(np.inf) # Error vector\n",
    "\n",
    "    iterations = 0\n",
    "    while np.max(abs(E)) > tol and iterations < n:\n",
    "        E = -np.dot(np.linalg.inv(J(X)), F(X))\n",
    "        X = X + E\n",
    "        iterations += 1\n",
    "        print(X)\n",
    "\n",
    "    print(\"iterations:\", iterations)\n",
    "\n",
    "    return X"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from math import inf\n",
    "\n",
    "f = lambda x, y: x**2 + x*y + y**2 - 7\n",
    "g = lambda x, y: x**3 + y**3 - 9\n",
    "\n",
    "F = lambda X: np.array([f(*X), g(*X)])\n",
    "\n",
    "fx = lambda x, y: 2*x + y\n",
    "fy = lambda x, y: x + 2*y\n",
    "gx = lambda x, y: 3*x**2\n",
    "gy = lambda x, y: 3*y**2\n",
    "\n",
    "J = lambda X: np.array([\n",
    "    [fx(*X), fy(*X)],\n",
    "    [gx(*X), gy(*X)]\n",
    "])\n",
    "\n",
    "X = np.array([1.5, 0.5])\n",
    "tol = 1e-6\n",
    "n = inf\n",
    "\n",
    "newton_raphson(F, J, X, tol, n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# For 3 variables, x,y and z\n",
    "\n",
    "import numpy as np\n",
    "from math import inf, sin, cos\n",
    "\n",
    "f = lambda x, y, z: 10*x + sin(x+y) - 1\n",
    "g = lambda x, y, z: 8*y - cos(z-y)**2 - 1\n",
    "h = lambda x, y, z: 12*z + sin(z) - 1\n",
    "\n",
    "F = lambda X: np.array([f(*X), g(*X), h(*X)])\n",
    "\n",
    "fx = lambda x, y, z: 10 + cos(x+y)\n",
    "fy = lambda x, y, z: cos(x+y)\n",
    "fz = lambda x, y, z: 0\n",
    "gx = lambda x, y, z: 0\n",
    "gy = lambda x, y, z: 8 - 2*cos(z-y)*sin(z-y)\n",
    "gz = lambda x, y, z: 2*cos(z-y)*sin(z-y)\n",
    "hx = lambda x, y, z: 0\n",
    "hy = lambda x, y, z: 0\n",
    "hz = lambda x, y, z: 12 + cos(z)\n",
    "\n",
    "J = lambda X: np.array([\n",
    "    [fx(*X), fy(*X), fz(*X)],\n",
    "    [gx(*X), gy(*X), gz(*X)],\n",
    "    [hx(*X), hy(*X), hz(*X)]\n",
    "])\n",
    "\n",
    "X = np.array([1/10, 1/4, 1/12])\n",
    "tol = 1e-6\n",
    "n = inf\n",
    "\n",
    "newton_raphson(F, J, X, tol, n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Using X instead of x, y, z, ... for math function lambdas\n",
    "# Use this as input method if there are 4 or more parameters\n",
    "\n",
    "import numpy as np\n",
    "from math import inf\n",
    "\n",
    "f = lambda X: X[0]**2 + X[0]*X[1] + X[1]**2 - 7\n",
    "g = lambda X: X[0]**3 + X[1]**3 - 9\n",
    "\n",
    "F = lambda X: np.array([f(X), g(X)])\n",
    "\n",
    "fx = lambda X: 2*X[0] + X[1]\n",
    "fy = lambda X: X[0] + 2*X[1]\n",
    "gx = lambda X: 3*X[0]**2\n",
    "gy = lambda X: 3*X[1]**2\n",
    "\n",
    "J = lambda X: np.array([\n",
    "    [fx(X), fy(X)],\n",
    "    [gx(X), gy(X)]\n",
    "])\n",
    "\n",
    "X = np.array([1.5, 0.5])\n",
    "tol = 1e-6\n",
    "n = inf\n",
    "\n",
    "newton_raphson(F, J, X, tol, n)"
   ]
  }
 ]
}