A6_Abhishek_1911007.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": [],
      "authorship_tag": "ABX9TyMz8d5y6XPkP2JKgqJxC7Ih",
      "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/github/physicsme1729/Numerical-methods-in-physics/blob/main/7_swarup_kumar_Giri_phy_P745.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "# ***Problem 1:***"
      ],
      "metadata": {
        "id": "ovHI1Xl697Me"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "Use Gaussian elimination or LU decomposition (without and with pivoting) for:\n",
        "$$\\left(\\begin{array}{cccc|c}  \n",
        " 4 & 4 & 8 & 4 & 1 \\\\  \n",
        " 4 & 5 & 3 & 7 & 2  \\\\\n",
        " 8 & 3 & 9 & 9 & 3\\\\\n",
        " 4 & 7 & 9 & 5 & 4\\\\\n",
        "\\end{array}\\right)\n",
        "$$\n",
        "Does pivoting help in any way? What are your conclusions about this matrix? You\n",
        "sometimes hear the advice that, since the problem arises from one of the $U_{ii}$'s being zero, you should replace it with a small number, say $10^{-20}$. Does this work?"
      ],
      "metadata": {
        "id": "H8QtdgIN-F_G"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "ANS:-"
      ],
      "metadata": {
        "id": "3Y6xYNHvAXm9"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "# ***I do the whole calculation  using Gaussian elimination for both cases with and without pivoting ...***"
      ],
      "metadata": {
        "id": "wa86LKDPIdJs"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "Let's first perform Gaussian elimination without pivoting and then with pivoting for the given matrix:"
      ],
      "metadata": {
        "id": "FLslzf_pCDOn"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "# **Gaussian Elimination Without Pivoting:**"
      ],
      "metadata": {
        "id": "YPpAuFLKJdv_"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "Step 1: Divide the first row by 4 to make the first element 1.\n",
        "$$\\left(\\begin{array}{cccc|c}  \n",
        " 1 & 1 & 2 & 1 & 1/4 \\\\  \n",
        " 4 & 5 & 3 & 7 & 2  \\\\\n",
        " 8 & 3 & 9 & 9 & 3\\\\\n",
        " 4 & 7 & 9 & 5 & 4\\\\\n",
        "\\end{array}\\right)\n",
        "$$\n",
        "Step 2: Subtract 4 times the first row from the second row, 8 times the first row from the third row, and 4 times the first row from the fourth row to make the first column below the 1 in the first row all zeros.\n",
        "$$\\left(\\begin{array}{cccc|c}  \n",
        " 1 & 1 & 2 & 1 & 1/4 \\\\  \n",
        " 0 & 1 & -5 & 3 & 1  \\\\\n",
        " 0 & -5 & -7 & 1 & 1\\\\\n",
        " 0 & 3 & 1 & 1 & 3\\\\\n",
        "\\end{array}\\right)\n",
        "$$\n",
        "Step 3:  Add 5 times the second row to the third row and subtract 3 times the second row from the fourth row to make all the elements below the 1 in the second row zero.\n",
        "$$\\left(\\begin{array}{cccc|c}  \n",
        " 1 & 1 & 2 & 1 & 1/4 \\\\  \n",
        " 0 & 1 & -5 & 3 & 1  \\\\\n",
        " 0 & 0 & -32 & 16 & 6\\\\\n",
        " 0 & 0 & 16 & -8 & 0\\\\\n",
        "\\end{array}\\right)\n",
        "$$\n",
        " step 4:Divide the third row by -32 to make the leading coefficient in the third row 1.\n",
        " $$\\left(\\begin{array}{cccc|c}  \n",
        " 1 & 1 & 2 & 1 & 1/4 \\\\  \n",
        " 0 & 1 & -5 & 3 & 1  \\\\\n",
        " 0 & 0 & 1 & -1/2 & -3/16\\\\\n",
        " 0 & 0 & 16 & -8 & 0\\\\\n",
        "\\end{array}\\right)\n",
        "$$\n",
        "Step 5: Subtract 16 times the third row from the fourth row to make all the elements below the 1 in the third row zero.\n",
        "$$\\left(\\begin{array}{cccc|c}  \n",
        " 1 & 1 & 2 & 1 & 1/4 \\\\  \n",
        " 0 & 1 & -5 & 3 & 1  \\\\\n",
        " 0 & 0 & 1 & -1/2 & -3/16\\\\\n",
        " 0 & 0 & 0 & 0 & 3\\\\\n",
        "\\end{array}\\right)\n",
        "$$\n"
      ],
      "metadata": {
        "id": "sMJZW1sxJnLe"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "Now, let's proceed with the conclusion and the advice given:\n",
        "\n",
        "# ***1.CONCLUSOIN:***\n",
        "In the without pivoting case, we've successfully reduced the matrix to row-echelon form. However, we have reached a point where the last row indicates that there are inconsistent equations. The system is inconsistent because the equation  cannot be satisfied, which means there is no solution to the system.\n",
        "\n",
        "# ***2.ADVICE:***\n",
        "The problem arose because one of the diagonal elements ($U_{ii}$) became zero during the elimination process. The advice to replace it with a small number (e.g., $10^{-20}$) is related to numerical stability. This can help avoid division by very small numbers, which might lead to numerical errors in some cases. However, in practice, it won't change the fundamental issue of inconsistency in this particular system."
      ],
      "metadata": {
        "id": "VYWsSXbvVfqB"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "# ***Gaussian Elimination With Pivoting:***"
      ],
      "metadata": {
        "id": "Md5bj_0xAd-o"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "Original matrix:\n",
        "\n",
        "$$\\left(\\begin{array}{cccc|c}  \n",
        " 4 & 4 & 8 & 4 & 1 \\\\\n",
        " 4 & 5 & 3 & 7 & 2 \\\\\n",
        " 8 & 3 & 9 & 9 & 3 \\\\\n",
        " 4 & 7 & 9 & 5 & 4 \\\\\n",
        "\\end{array}\\right)$$\n",
        "\n",
        "\n",
        "Step 1:\n",
        "\n",
        "Swap the first and third rows to pivot on the maximum element in the first column (which is 8):\n",
        "\n",
        "$$\\left(\\begin{array}{cccc|c}  \n",
        " 8 & 3 & 9 & 9 & 3 \\\\\n",
        " 4 & 5 & 3 & 7 & 2 \\\\\n",
        " 4 & 4 & 8 & 4 & 1 \\\\\n",
        " 4 & 7 & 9 & 5 & 4 \\\\\n",
        "\\end{array}\\right)$$\n",
        "\n",
        "\n",
        "Step 2:\n",
        "new Second row=2× old second row-first row\n",
        "new third row=2× old third row-first row\n",
        "new fourth row=2× old fourth row-first row\n",
        "$$\\left(\\begin{array}{cccc|c}  \n",
        "8 & 3 & 9 & 9 & 3 \\\\\n",
        " 0 & 7 & -3 & 5 & 1 \\\\\n",
        " 0 & 5 & 7 & -1 & -1 \\\\\n",
        " 0 & 11 & 9 & 1 & 5 \\\\\n",
        "\\end{array}\\right)$$\n",
        "\n",
        "Step 3:\n",
        "Swap the fourth and second rows to pivot on the maximum element in the second column (which is 11):\n",
        "\n",
        "$$\\left(\\begin{array}{cccc|c}  \n",
        "8 & 3 & 9 & 9 & 3 \\\\\n",
        " 0 & 11 & 9 & 1 & 5 \\\\\n",
        " 0 & 5 & 7 & -1 & -1 \\\\\n",
        " 0 & 7 & -3 & 5 & 1 \\\\\n",
        " \\end{array}\\right)$$\n",
        "\n",
        "Step 4:\n",
        "\n",
        "Subtract $\\frac{11}{5}$ times the third row from the second row:\n",
        "Subtract $\\frac{11}{7}$ times the fourth row from the second row:\n",
        "\n",
        "$$\\left(\\begin{array}{cccc|c}\n",
        "8 & 3 & 9 & 9 & 3 \\\\\n",
        " 0 & 11 & 9 & 1 & 5 \\\\\n",
        " 0 & 0 & \\frac{32}{5} & -\\frac{16}{5} & -\\frac{36}{5} \\\\\n",
        " 0 & 0 & -\\frac{96}{7} & \\frac{48}{7} & -\\frac{24}{7} \\\\\n",
        "\\end{array}\\right)$$\n",
        "\n",
        "Step 5:\n",
        "\n",
        "multiply 5 in third row and 7 in 4th row\n",
        "\n",
        "$$\\left(\\begin{array}{cccc|c}\n",
        "8 & 3 & 9 & 9 & 3 \\\\\n",
        " 0 & 11 & 9 & 1 & 5 \\\\\n",
        " 0 & 0 & 32 & -16 & -36 \\\\\n",
        " 0 & 0 & -96 & 48 & -24 \\\\\n",
        "\\end{array}\\right)$$\n",
        "\n",
        "Step 6:\n",
        "\n",
        "multiply 3 in third row and add with  4th row\n",
        "\n",
        "$$\\left(\\begin{array}{cccc|c}\n",
        "8 & 3 & 9 & 9 & 3 \\\\\n",
        " 0 & 11 & 9 & 1 & 5 \\\\\n",
        " 0 & 0 & 32 & -16 & -36 \\\\\n",
        " 0 & 0 & 0 & 0 & -132\\\\\n",
        "\\end{array}\\right)$$\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "Now, let's solve for the variables using back-substitution:\n",
        "\n",
        "From the fourth row, we have $0x_1 + 0x_2 + 0x_3 + 0x_4 = -132$, which implies $0 = -132$. This is not possible, indicating that there is no consistent solution to the system of equations.\n",
        "\n",
        "***Conclusion with Pivoting:***\n",
        "With pivoting, we've also successfully reduced the matrix to row-echelon form. However, in this case, the system is inconsistent as well, as indicated by the last row.Pivoting didn't change the fundamental inconsistency issue.\n",
        "\n",
        "In summary, both without and with pivoting, the system of equations is inconsistent, and there is no solution. The advice to replace a zero pivot with a small number can improve numerical stability but does not address the underlying issue of inconsistency in this particular system."
      ],
      "metadata": {
        "id": "oD0Lk-dDAvHX"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "# ***Problem 2:***"
      ],
      "metadata": {
        "id": "F2OmmJw3CD1U"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "As mentioned in the main text, a small modification to the Jacobi method of Equation\n",
        "\\begin{align*}\n",
        "x^{(k)}_i &= \\frac{1}{A_{ii}} \\left( b_i - \\sum_{j=0 }^{i-1} A_{ij}x^{(k-1)}_j- \\sum_{j=i+1 }^{n-1} A_{ij}x^{(k-1)}_j\\right), \\quad i = 0, 1, \\ldots, n-1\n",
        "\\end{align*}\n",
        "\n",
        "leads to an improved algorithm. This is the Gauss–Seidel method, which is given by:\n",
        "\\begin{align*}\n",
        "x^{(k)}_i &= \\frac{1}{A_{ii}} \\left( b_i - \\sum_{j=0 }^{i-1} A_{ij}x^{(k)}_j- \\sum_{j=i+1 }^{n-1} A_{ij}x^{(k-1)}_j\\right), \\quad i = 0, 1, \\ldots, n-1\n",
        "\\end{align*}\n",
        "\n",
        "In words, what the Gauss–Seidel method does is to use the improved values as soon as they become available. Implement the Gauss–Seidel method in Python and compare its convergence with that of the Jacobi method, for the same problem as in jacobi.py.\n"
      ],
      "metadata": {
        "id": "Ed36SDGmCMJ2"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "# ***ANS:-***\n",
        "To implement the Gauss-Seidel method in Python and compare its convergence with the Jacobi method for the same problem as  provided \"jacobi.py\" code, One can make a few modifications to our existing code. Here's the Python code for the Gauss-Seidel method and the comparison:"
      ],
      "metadata": {
        "id": "ekYqJa7tY8GV"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "\n",
        "\n",
        "\n",
        "def termcrit(xolds, xnews):\n",
        "    errs = np.abs((xnews - xolds) / xnews)\n",
        "    return np.sum(errs)\n",
        "\n",
        "def jacobi(A, bs, kmax=50, tol=1.e-6):\n",
        "    n = bs.size\n",
        "    xnews = np.zeros(n)\n",
        "    for k in range(1, kmax):\n",
        "        xs = np.copy(xnews)\n",
        "        for i in range(n):\n",
        "            slt = A[i, :i] @ xs[:i]\n",
        "            sgt = A[i, i + 1:] @ xs[i + 1:]\n",
        "            xnews[i] = (bs[i] - slt - sgt) / A[i, i]\n",
        "        err = termcrit(xs, xnews)\n",
        "        print(k, xnews, err)\n",
        "        if err < tol:\n",
        "            break\n",
        "    return xnews\n",
        "\n",
        "def gauss_seidel(A, bs, kmax=50, tol=1.e-6):\n",
        "    n = bs.size\n",
        "    xnews = np.zeros(n)\n",
        "    for k in range(1, kmax):\n",
        "        xs = np.copy(xnews)\n",
        "        for i in range(n):\n",
        "            slt = A[i, :i] @ xnews[:i]\n",
        "            sgt = A[i, i + 1:] @ xs[i + 1:]\n",
        "            xnews[i] = (bs[i] - slt - sgt) / A[i, i]\n",
        "        err = np.linalg.norm(xnews - xs, np.inf)\n",
        "        print(k, xnews, err)\n",
        "        if err < tol:\n",
        "            break\n",
        "    return xnews\n",
        "\n",
        "if __name__ == '__main__':\n",
        "    n = 4\n",
        "    val = 21\n",
        "    A, bs = testcreate(n, val)\n",
        "    A += val * np.identity(n)\n",
        "\n",
        "    print(\"Jacobi:\")\n",
        "    result_jacobi = jacobi(A, bs)\n",
        "\n",
        "    print(\"\\nGauss-Seidel:\")\n",
        "    result_gauss_seidel = gauss_seidel(A, bs)\n",
        "\n",
        "    print(\"\\nSolution (Jacobi):\", result_jacobi)\n",
        "    print(\"Solution (Gauss-Seidel):\", result_gauss_seidel)\n",
        "\n",
        "\n"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "VSHEcc6mmgXz",
        "outputId": "53289af9-60d8-4215-be5a-2f19ba8cc268"
      },
      "execution_count": 27,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Jacobi:\n",
            "1 [0.95584174 0.98382901 1.01264965 1.04197678] 4.0\n",
            "2 [0.38609117 0.38784058 0.39421864 0.40425572] 6.158643474272702\n",
            "3 [0.73341768 0.74941422 0.76835846 0.78969405] 1.9310672038734116\n",
            "4 [0.52317706 0.53023862 0.54134654 0.55573114] 1.6555533665830167\n",
            "5 [0.65072149 0.6631481  0.67896251 0.6975372 ] 0.8024076458188099\n",
            "6 [0.57339977 0.58256412 0.59551555 0.61154376] 0.5539167512173229\n",
            "7 [0.6202852  0.63142596 0.64611142 0.66368227] 0.30983792417040845\n",
            "8 [0.59185737 0.60179949 0.61543315 0.63206832] 0.19712628238915642\n",
            "9 [0.60909428 0.61976313 0.63403445 0.65123688] 0.11605599760545227\n",
            "10 [0.59864294 0.60887114 0.62275581 0.63961427] 0.07162932168524415\n",
            "11 [0.60497996 0.61547534 0.62959445 0.64666147] 0.04296480237007709\n",
            "12 [0.6011376  0.61147098 0.62544794 0.6423885 ] 0.026221914973545257\n",
            "13 [0.60346736 0.61389897 0.62796212 0.64497935] 0.015836323357587702\n",
            "14 [0.60205475 0.61242679 0.62643768 0.64340843] 0.009625240606544867\n",
            "15 [0.60291127 0.61331942 0.627362   0.64436094] 0.005827622959638275\n",
            "16 [0.60239193 0.61277819 0.62680155 0.6437834 ] 0.0035366211778694286\n",
            "17 [0.60270682 0.61310636 0.62714137 0.64413358] 0.002143230680561838\n",
            "18 [0.60251589 0.61290738 0.62693533 0.64392125] 0.0012999398547222725\n",
            "19 [0.60263166 0.61302803 0.62706026 0.64404999] 0.0007880447138564145\n",
            "20 [0.60256146 0.61295487 0.62698451 0.64397193] 0.00047787679848529385\n",
            "21 [0.60260403 0.61299923 0.62703044 0.64401926] 0.0002897328047164427\n",
            "22 [0.60257822 0.61297234 0.62700259 0.64399056] 0.00017568308149790347\n",
            "23 [0.60259387 0.61298864 0.62701948 0.64400797] 0.00010652009946952199\n",
            "24 [0.60258438 0.61297875 0.62700924 0.64399742] 6.458798988974602e-05\n",
            "25 [0.60259013 0.61298475 0.62701545 0.64400381] 3.916162521264781e-05\n",
            "26 [0.60258664 0.61298111 0.62701168 0.64399993] 2.3745235251712988e-05\n",
            "27 [0.60258876 0.61298332 0.62701397 0.64400229] 1.4397533714587382e-05\n",
            "28 [0.60258748 0.61298198 0.62701258 0.64400086] 8.729758777493622e-06\n",
            "29 [0.60258825 0.61298279 0.62701342 0.64400172] 5.2931580222949495e-06\n",
            "30 [0.60258778 0.6129823  0.62701291 0.6440012 ] 3.2094336318618936e-06\n",
            "31 [0.60258807 0.6129826  0.62701322 0.64400152] 1.9459934721551643e-06\n",
            "32 [0.60258789 0.61298242 0.62701303 0.64400132] 1.1799258029178294e-06\n",
            "33 [0.602588   0.61298253 0.62701315 0.64400144] 7.154310950193359e-07\n",
            "\n",
            "Gauss-Seidel:\n",
            "1 [0.95584174 0.80071068 0.65382991 0.52242462] 0.9558417403307491\n",
            "2 [0.58642344 0.63538954 0.65155646 0.63722374] 0.36941829941956716\n",
            "3 [0.59517664 0.61089005 0.63038983 0.64529024] 0.024499489231789906\n",
            "4 [0.60209177 0.61214395 0.62701213 0.64428827] 0.006915127226247431\n",
            "5 [0.60268695 0.61290555 0.62694782 0.64401126] 0.0007616018591185547\n",
            "6 [0.60261242 0.6129888  0.62700475 0.64399666] 8.324690473626717e-05\n",
            "7 [0.60258928 0.61298486 0.62701336 0.64400055] 2.3139989749121348e-05\n",
            "8 [0.60258764 0.61298267 0.62701331 0.64400138] 2.189140510089338e-06\n",
            "9 [0.60258789 0.61298246 0.62701313 0.64400141] 2.507116306782464e-07\n",
            "\n",
            "Solution (Jacobi): [0.602588   0.61298253 0.62701315 0.64400144]\n",
            "Solution (Gauss-Seidel): [0.60258789 0.61298246 0.62701313 0.64400141]\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "From this two solution one conclude that two solution are equal. But in Gauss-Seidal case ineration number is 9 ,in case of other itertion number is 33, So Gauss-Seidel method is more appropriate than the Jacobi method for this particular problem. The Gauss-Seidel method achieved convergence with significantly fewer iterations, which indicates faster convergence for this specific linear system."
      ],
      "metadata": {
        "id": "EMCLsVivnb1q"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "# ***Problem 3:***"
      ],
      "metadata": {
        "id": "y5xbvKr3o3nU"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "This problem studies the modified Gram–Schmidt method, which is less sensitive to\n",
        "numerical errors than classical Gram–Schmidt. Conceptually, classical Gram–Schmidt\n",
        "produces new $\\textbf{a}_j^{'}$ vectors by subtracting out any non-orthogonal components of the original vectors  $\\textbf{a}_j$, as per Equation-\n",
        "\n",
        "\n",
        "\\begin{align*}\n",
        "    \\textbf{a}_j' &= \\textbf{a}_j - \\sum_{i=0}^{j-1} (\\textbf{q}_i^{T} \\textbf{a}_j) \\textbf{q}_i \\\\\n",
        "    \\textbf{q}_j &= \\frac{\\textbf{a}_j'}{\\|\\textbf{a}_j'\\|}\n",
        "\\end{align*}\n",
        "The problem with this is that if in the calculation of,\n",
        "say, $\\textbf{q}_3$ a numerical roundoff error is introduced, this error will then be propagated to the computation of $\\textbf{q}_4$,$\\textbf{q}_5$, and so on. In contradistinction to this, modified Gram–Schmidt tries to correct for the error in $\\textbf{q}_3$ when computing $\\textbf{q}_4$,$\\textbf{q}_5$, and so on.\n",
        "\n",
        "In equation form, what modified Gram–Schmidt does is to replace the relations from\n",
        "as follows:\n",
        "\\begin{align*}\n",
        "\\textbf{q}_0 &= \\frac{\\textbf{a}_0}{\\|\\textbf{a}_0\\|}\\\\\n",
        "    \\textbf{a}_j'^{(0)} &= \\textbf{a}_j -  (\\textbf{q}_0^{T} \\textbf{a}_j) \\textbf{q}_0, \\quad  j=1,2,\\dots,n-1 \\\\\n",
        "    \\textbf{q}_1 &= \\frac{\\textbf{a}_1'^{(0)}}{\\|\\textbf{a}_1'^{(0)}\\|}\\\\\n",
        "\\end{align*}\n",
        "So far, there’s nothing really “modified” going on. Then, the next step is:\n",
        "\n",
        "\\begin{align*}\n",
        "    \\textbf{a}_j'^{(1)} &= \\textbf{a}_j'^{(0)}-  (\\textbf{q}_1^{T}  \\textbf{a}_j'^{(0)}) \\textbf{q}_1, \\quad  j=2,3,\\dots,n-1 \\\\\n",
        "    \\textbf{q}_2 &= \\frac{\\textbf{a}_2'^{(1)}}{\\|\\textbf{a}_2'^{(1)}\\|}\\\\\n",
        "\\end{align*}\n",
        "Note how the inner product is taken not with the original vector, but with the updated\n",
        "one (which has a superscript in parentheses).\n",
        "Thus, for the general case we have:\n",
        "\n",
        "\\begin{align*}\n",
        "\\textbf{q}_i &= \\frac{\\textbf{a}_i'^{(i-1)}}{\\|\\textbf{a}_i'^{(i-1)}\\|}\\\\\n",
        "    \\textbf{a}_j'^{(i)} &= \\textbf{a}_j'^{(i-1)} -  (\\textbf{q}_i^{T} \\textbf{a}_j'^{(i-1)}) \\textbf{q}_i, \\quad  j=i+1,\\dots,n-1 \\\\\n",
        "    \\end{align*}\n",
        "\n",
        "Implement QR decomposition in Python using the modified Gram–Schmidt approach.\n",
        "You should carefully think about how to structure your new code: the prescription above builds the required matrices row by row (whereas \"qrdec.py\" works column by column), so it may help you to first restructure your classical Gram–Schmidt code to also work row by row. At that point, you will see that the only difference between the two methods is whether the inner products are computed using the original or the updated vectors.\n",
        "\n",
        "When you’re done, evaluate $\\|\\textbf{Q}^T\\textbf{Q} - \\textbf{I}\\|$\n",
        " for the $\\textbf{testcreate()}$ problem and compare the orthogonality properties of our new method with those of classical Gram–Schmidt."
      ],
      "metadata": {
        "id": "L7NqNhcio-QK"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "# ***ANS:***"
      ],
      "metadata": {
        "id": "FS9AtAc7pGd7"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "To implement the QR decomposition using the modified Gram-Schmidt approach in Python, ONE can follow the given prescription. The code will build the Q and R matrices row by row, just like the modified Gram-Schmidt method. Here's the implementation:"
      ],
      "metadata": {
        "id": "ajl5LDfjmKZM"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "\n",
        "def mag(xs):\n",
        "    return np.sqrt(np.sum(xs*xs))\n",
        "def testcreate(n,val):\n",
        "    A = np.arange(val,val+n*n).reshape(n,n)\n",
        "    A = np.sqrt(A)\n",
        "    bs = (A[0,:])**2.1\n",
        "    return A, bs\n",
        "\n",
        "\n",
        "def qrdec(A):\n",
        "    n = A.shape[0]\n",
        "    Ap = np.copy(A)\n",
        "    Q = np.zeros((n,n))\n",
        "    R = np.zeros((n,n))\n",
        "    for j in range(n):\n",
        "        for i in range(j):\n",
        "            R[i,j] = Q[:,i] @ Ap[:,j]\n",
        "            Ap[:,j] -= R[i,j] * Q[:,i]\n",
        "        R[j,j] = mag(Ap[:,j])\n",
        "        Q[:,j] = Ap[:,j] / R[j,j]\n",
        "    return Q, R\n",
        "\n",
        "def testqrdec(A):\n",
        "    n = A.shape[0]\n",
        "    Q, R = qrdec(A)\n",
        "    diffa = A - Q @ R\n",
        "    diffq = np.transpose(Q) @ Q - np.identity(n)\n",
        "    print(\"n=\",n, \"mag(diffa)=\",mag(diffa),\"mag(diffq)=\", mag(diffq))\n",
        "\n",
        "if __name__ == '__main__':\n",
        "    for n in range(4, 10, 2):\n",
        "        A, bs = testcreate(n, 21)\n",
        "        testqrdec(A)\n",
        "\n"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "rr8FuWSVmK-t",
        "outputId": "ad5f71de-1f6d-4080-e762-1d894ebd606a"
      },
      "execution_count": 23,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "n= 4 mag(diffa)= 1.2560739669470201e-15 mag(diffq)= 3.749290521725868e-08\n",
            "n= 6 mag(diffa)= 2.5121479338940403e-15 mag(diffq)= 0.00017473495623754012\n",
            "n= 8 mag(diffa)= 5.102196573270515e-15 mag(diffq)= 1.0372175137433681\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        " in python code ,$\\|\\textbf{Q}^T\\textbf{Q} - \\textbf{I}\\|$=diffq\n",
        "\n",
        " compare the orthogonality , $\\|\\textbf{Q}^T\\textbf{Q} - \\textbf{I}\\|$ various matrix sizes, allowing  to compare the orthogonality properties of the modified Gram-Schmidt method with the classical Gram-Schmidt method.\n",
        "\n"
      ],
      "metadata": {
        "id": "2eerjFCamwXH"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "# *PROBLEM 4:-*"
      ],
      "metadata": {
        "id": "6qoRR-5drsZl"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "We now use Householder transformations to carry out a QR decomposition.\n",
        "\n",
        "(a) We define a \"Householder matrix\" as $\\textbf{P}=\\textbf{I}-2\\textbf{w}(\\textbf{w}^T)$, where $\\textbf{w}$ is a unit-norm vector.\n",
        "\n",
        "(Our definition involves the product of a column vector and a row vector);\n",
        "Show that $\\textbf{P}$ is symmetric and orthogonal.You have thereby shown that $\\textbf{P}^2=\\textbf{I}$.\n",
        "\n",
        "(b)  If two vectors $\\mathbf{x}$ and $\\mathbf{y}$ satisfy $\\mathbf{P}\\mathbf{x} = \\mathbf{y}$, then we can show that\n",
        "\n",
        "$\\mathbf{w} = \\frac{(\\mathbf{x} - \\mathbf{y})}{\\|\\mathbf{x} - \\mathbf{y}\\|}$\n",
        "\n",
        "(c) Specialize to $\\mathbf{y}$ = $\\alpha \\|\\mathbf{x}\\| \\mathbf{e}_0$, where $\\mathbf{e}_0$ is the first column of the identity matrix, as per equation\n",
        "\\begin{align*}\n",
        "    \\mathbf{A}^{-1}&=(\\mathbf{x}_0,\\mathbf{x}_1,\\dots,\\mathbf{x}_{n-1}),  \\quad  & \\mathbf{I}= (\\mathbf{e}_0,\\mathbf{e}_1,\\dots,\\mathbf{e}_{n-1})\n",
        "  \\end{align*}\n",
        "\\text{ Keep in mind that bold-lower-case entities are column vectors}\n",
        "\n",
        "(Discuss why you should choose $\\alpha$ = -sign$(x_0)$.) ) Notice that $\\mathbf{P}\\mathbf{x} $ zeroes out all the elements of a general vector $\\mathbf{x}$ below the leading element; this is reminiscent of the first step in Gaussian elimination.\n",
        "\n",
        "(d) Produce a Householder matrix $\\mathbf{P}^{(0)}$ that zeroes out the all the elements except the leading one in the first column of an n × n matrix $\\mathbf{A}$. This $\\mathbf{P}^{(0)}$ acting on the entire $\\mathbf{A}$ thereby transforms the other elements (outside the first column). We now need to zero out the elements in $\\mathbf{P}^{(0)}$$\\mathbf{A}$ below its 1, 1 element. We form a new Householder\n",
        "matrix $\\mathbf{P}^{(1)}$  based on the lower n - 1 elements of the second column of $\\mathbf{P}^{(0)}$$\\mathbf{A}$; note that this leads to $\\mathbf{P}^{(1)}$ being $(n-1)\\times(n-1)$. Similarly,  $\\mathbf{P}^{(2)}$ would be $(n-2)\\times(n-2)$.,\n",
        "and so on. We opt for the (wasteful) option of padding:\n",
        "\n",
        " $\\mathbf{Q}^{(k)}$=\n",
        "\\begin{pmatrix}\n",
        "    \\mathbf{I} & 0 \\\\\n",
        "    0 & \\mathbf{P}^{(k)} \\\\\n",
        "\\end{pmatrix}\n",
        "\n",
        "where $\\mathbf{I}$ is $k \\times k$, $\\mathbf{P}^{(k)}$ is $(n-k)\\times(n-k)$ and, therefore, $\\mathbf{Q}^{(k)}$ is $n \\times n$. Explain why $ \\mathbf{R}=\\mathbf{Q}^{(n-2)} \\mathbf{Q}^{(n-3)} \\dots \\mathbf{Q}^{(1)} \\mathbf{Q}^{(0)} \\mathbf{A}$ is upper triangular. Use the fact that these matrices are orthogonal to show that  $ \\mathbf{Q}=\\mathbf{Q}^{(0)} \\mathbf{Q}^{(1)} \\dots \\mathbf{Q}^{(n-3)} \\mathbf{Q}^{(n-2)} $ is also orthogonal.\n",
        "Multiplying $\\mathbf{Q}$ and $\\mathbf{R}$ gives you the initial $\\mathbf{A}$.\n",
        "\n",
        "(e) Implement this prescription in Python to produce a QR decomposition.\n",
        "\n"
      ],
      "metadata": {
        "id": "2Ldq1a3LrxIb"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "# ***ANS OF PART (a):***"
      ],
      "metadata": {
        "id": "T14ylW8jvPhq"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "$\\textbf{Symmetry:}$\n",
        "\n",
        "Taking the transpose of $\\mathbf{P}$\n",
        "\\begin{align*}\n",
        "\\mathbf{P}^T &= (\\mathbf{I} - 2\\mathbf{w}\\mathbf{w}^T)^T \\\\\n",
        "&= \\mathbf{I}^T - (2\\mathbf{w}\\mathbf{w}^T)^T \\\\\n",
        "&= \\mathbf{I} - 2(\\mathbf{w}\\mathbf{w}^T)^T \\\\\n",
        "&= \\mathbf{I} - 2\\mathbf{w}\\mathbf{w}^T\n",
        "\\end{align*}\n",
        "Since\n",
        "\\begin{align*}\n",
        "(\\mathbf{P}^T = \\mathbf{P}) \\end{align*},\n",
        " we have shown that $\\mathbf{P}$ is symmetric.\n",
        "\n"
      ],
      "metadata": {
        "id": "ECp9tCFFvZq5"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "\n",
        "\n",
        "$\\textbf{Orthogonality:}$\n",
        "\n",
        "To show that $\\mathbf{P}$ is orthogonal, we need to prove that $\\mathbf{P}^T \\mathbf{P} = \\mathbf{I}$.\n",
        "\n",
        "Let's calculate $\\mathbf{P}^T \\mathbf{P}$:\n",
        "\n",
        "$\\mathbf{P}^T \\mathbf{P} = (\\mathbf{I} - 2\\mathbf{w}(\\mathbf{w}^T))^T(\\mathbf{I} - 2\\mathbf{w}(\\mathbf{w}^T))= (\\mathbf{I} - 2\\mathbf{w}(\\mathbf{w}^T))(\\mathbf{I} - 2\\mathbf{w}(\\mathbf{w}^T))$\n",
        "\n",
        "\n",
        "Now, we apply the transpose to both terms:\n",
        "\n",
        "\n",
        "$\\mathbf{P}^T \\mathbf{P} =  (\\mathbf{I} - 4\\mathbf{w}(\\mathbf{w}^T)+4\\mathbf{w}(\\mathbf{w}^T).\\mathbf{w}(\\mathbf{w}^T))=\\mathbf{I},$\n",
        "As,$\\mathbf{w}$ is a unit-norm vector,$\\mathbf{w}^T.\\mathbf{w}=\\mathbf{I}$\n",
        "\n",
        "\n",
        "from symmertic case $\\mathbf{P}^T = \\mathbf{P}$\n",
        "\n",
        "so $\\mathbf{P}^2 = \\mathbf{I}$(Proved)\n"
      ],
      "metadata": {
        "id": "dNZhLDqasIgo"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "# ***ANS OF PART (b):***"
      ],
      "metadata": {
        "id": "jYD4k-6m3O9Q"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "To show that if two vectors $\\mathbf{x}$ and $\\mathbf{y}$ satisfy $\\mathbf{P}\\mathbf{x} = \\mathbf{y}$, then we can demonstrate that\n",
        "\n",
        "$\\mathbf{w} = \\frac{\\mathbf{x} - \\mathbf{y}}{\\|\\mathbf{x} - \\mathbf{y}\\|}$.\n",
        "\n",
        "\n",
        "Starting with the given equation $\\mathbf{P} = \\mathbf{I} - 2\\mathbf{w}(\\mathbf{w}^T)$, we rewrite $\\mathbf{P}\\mathbf{x} = \\mathbf{y}$ as follows:\n",
        "\n",
        "\n",
        "$(\\mathbf{I} - 2\\mathbf{w}(\\mathbf{w}^T))\\mathbf{x} = \\mathbf{y}.$\n",
        "\n",
        "Simplifying further:\n",
        "\n",
        "\n",
        "$\\mathbf{x} - 2\\mathbf{w}(\\mathbf{w}^T\\mathbf{x}) = \\mathbf{y}.$\n",
        "\n",
        "Isolating the term involving $\\mathbf{w}$:\n",
        "\n",
        "$- 2\\mathbf{w}(\\mathbf{w}^T\\mathbf{x}) = \\mathbf{y} - \\mathbf{x}.$\n",
        "\n",
        "Dividing both sides by $-2$:\n",
        "\n",
        "\n",
        "$\\mathbf{w}(\\mathbf{w}^T\\mathbf{x}) = \\frac{\\mathbf{x} - \\mathbf{y}}{-2}.$\n",
        "\n",
        "Isolating $\\mathbf{w}$:\n",
        "\n",
        "\n",
        "$\\mathbf{w} = \\frac{\\mathbf{x} - \\mathbf{y}}{-2(\\mathbf{w}^T\\mathbf{x})} = \\frac{\\mathbf{x} - \\mathbf{y}}{-2\\mathbf{w}^T\\mathbf{x}}.$\n",
        "To normalize $\\mathbf{w}$ to a unit-norm vector, we divide it by its magnitude $\\|\\mathbf{w}\\|$:\n",
        "\n",
        "\n",
        "$\\mathbf{w} = \\frac{\\frac{\\mathbf{x} - \\mathbf{y}}{-2\\mathbf{w}^T\\mathbf{x}}}{\\left\\|\\frac{\\mathbf{x} - \\mathbf{y}}{-2\\mathbf{w}^T\\mathbf{x}}\\right\\|}.$\n",
        "\n",
        "Notice that the numerator and denominator have the same term in the denominator $(-2\\mathbf{w}^T\\mathbf{x})$, which cancels out:\n",
        "\n",
        "\n",
        "$\\mathbf{w} = \\frac{\\mathbf{x} - \\mathbf{y}}{\\|\\mathbf{x} - \\mathbf{y}\\|}.$\n",
        "\n",
        "This completes the proof.\n"
      ],
      "metadata": {
        "id": "aBIYTlBF41ZC"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "# ***ANS OF PART (c):***"
      ],
      "metadata": {
        "id": "9pChCICd6mn6"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "So, this is the expression for $\\mathbf{w}$ specialized to $\\mathbf{y} = \\alpha |\\mathbf{x}| \\mathbf{e}_0$.\n",
        "\n",
        "Now, let's discuss why we should choose $\\alpha = -\\text{sign}(x_0)$. The Householder transformation is often used to zero out elements below the leading element in a vector. By choosing $\\alpha$ as the negative of the sign of $x_0$, we ensure that the leading element of $\\mathbf{y}$ has the opposite sign of the leading element of $\\mathbf{x}$. This choice helps in achieving the goal of the Householder transformation, which is to make the leading element of the transformed vector (in this case, $\\mathbf{y}$) equal to its norm. This is particularly useful in numerical stability and ensuring that the transformation behaves as expected in the context of QR decomposition or other linear algebra operations."
      ],
      "metadata": {
        "id": "V9_hZBxkAHHy"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "# ***ANS OF PART (d):***"
      ],
      "metadata": {
        "id": "xsCNVngmCEWW"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "Product matrix R = Q^(n-2) Q^(n-3) ... Q^(1) Q^(0) A is upper triangular, and that Q = Q^(0) Q^(1) ... Q^(n-3) Q^(n-2) is orthogonal, we need to understand the properties of Householder matrices and how they are applied.\n",
        "\n",
        "Recall that a Householder matrix P is symmetric and orthogonal, and it zeroes out all elements below the leading element when applied to a vector. It is represented as P = I - 2ww^T, where w is a unit-norm vector.\n",
        "\n",
        "Now, let's go through the steps involved in the Householder QR decomposition and how it results in an upper triangular matrix R and an orthogonal matrix Q:\n",
        "\n",
        "1. First Householder Transformation (P^(0)):\n",
        "   - P^(0) is designed to zero out all elements below the leading element in the first column of A.\n",
        "   - When one apply P^(0) to A, one can  get P^(0)A, and this operation zeroes out the elements below the leading element in the first column.\n",
        "\n",
        "2. Subsequent Householder Transformations (P^(1), P^(2), ...):\n",
        "   - Each P^(k) is designed to zero out all elements below the leading element in the k-th column of P^(k-1)A.\n",
        "   - The result of applying P^(k) to P^(k-1)A is P^(k)(P^(k-1)A), which zeroes out the elements below the leading element in the k-th column.\n",
        "\n",
        "3. Forming Q^(k):\n",
        "   - To maintain the same dimensions for Q^(k) as for P^(k), we pad P^(k) with identity matrices, creating Q^(k).\n",
        "\n",
        "4. Multiplying All Q^(k) Matrices Together:\n",
        "   - By multiplying Q^(0) Q^(1) ... Q^(n-2) Q^(n-1), you are essentially combining all the Householder transformations into a single matrix Q.\n",
        "\n",
        "5. Resulting R:\n",
        "   - The product R = Q^(n-2) Q^(n-3) ... Q^(0) A is upper triangular because each P^(k) is designed to zero out elements below the leading element in the k-th column.\n",
        "\n",
        "6. Orthogonality of Q:\n",
        "   - Since each Q^(k) is constructed from orthogonal Householder matrices, their product Q is also orthogonal. Orthogonal matrices have the property that Q^T Q = I, which implies that Q^(-1) = Q^T.\n",
        "\n",
        "So, in conclusion, the Householder QR decomposition using these Householder matrices Q^(k) results in an upper triangular matrix R, and the product of all the Q^(k) matrices is an orthogonal matrix Q. This decomposition is a fundamental step in solving linear systems and various numerical algorithms.\n"
      ],
      "metadata": {
        "id": "GivnbjVGKNY_"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "# ***ANS OF PART (e):***"
      ],
      "metadata": {
        "id": "1-ep9ysmLqL6"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "First, create a function to perform the Householder transformation and generate the Householder matrix 'P'. Then, use these Householder matrices to calculate 'Q' and 'R'. Here's the code for the modified 'qrdec' function:"
      ],
      "metadata": {
        "id": "V1wAwuJhLayD"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "\n",
        "def qrdec(A):\n",
        "    n = A.shape[0]\n",
        "    Ap = np.copy(A)\n",
        "    Q = np.eye(n)  # Initialize Q as an identity matrix\n",
        "    for j in range(n):\n",
        "        x = Ap[j:, j]\n",
        "        norm_x = np.linalg.norm(x)\n",
        "        alpha = -np.sign(x[0]) * norm_x\n",
        "        r = np.sqrt(2 * (norm_x ** 2) - 2 * x[0] * alpha)\n",
        "        w = np.zeros_like(x)\n",
        "        w[0] = (x[0] - alpha) / r\n",
        "        w[1:] = x[1:] / r\n",
        "\n",
        "        # Calculate Householder matrix P\n",
        "        P = np.eye(n)\n",
        "        P[j:, j:] -= 2 * np.outer(w, w)\n",
        "\n",
        "        # Update Q and R\n",
        "        Q = Q @ P.T\n",
        "        Ap = P @ Ap\n",
        "\n",
        "    R = Ap\n",
        "    return Q, R\n",
        "\n",
        "if __name__ == '__main__':\n",
        "    A, bs = testcreate(4, 21)\n",
        "    Q, R = qrdec(A)\n",
        "\n",
        "    print(\"Q matrix:\")\n",
        "    print(Q)\n",
        "    print(\"R matrix:\")\n",
        "    print(R)\n"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "WTCUwzFHN0jH",
        "outputId": "02a3b4cc-e837-4c9f-a06f-209e04e629bd"
      },
      "execution_count": 28,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Q matrix:\n",
            "[[-0.44095855  0.74266424  0.48002199 -0.15357174]\n",
            " [-0.48112522  0.22782622 -0.60503318  0.59207143]\n",
            " [-0.51818773 -0.21035281 -0.38538773 -0.73396831]\n",
            " [-0.5527708  -0.5935459   0.50496538  0.29522478]]\n",
            "R matrix:\n",
            "[[-1.03923048e+01 -1.05829568e+01 -1.07701462e+01 -1.09540595e+01]\n",
            " [-1.95821368e-16  3.20069752e-02  6.28525184e-02  9.26292875e-02]\n",
            " [ 3.01250545e-17 -8.24602213e-18 -7.26503393e-05 -2.06740127e-04]\n",
            " [-2.20378793e-16  6.77687598e-18 -1.53150433e-20 -3.08315359e-07]]\n"
          ]
        }
      ]
    }
  ]
}