XM_40017/notebooks/LEQ.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "9c32b051",
   "metadata": {},
   "source": [
    "# Solving Linear Equations\n",
    "\n",
    "## Serial Algorithm\n",
    "To demonstrate the algorithm, we will consider a simple system of linear equations $Ax = b$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "id": "2b369b73",
   "metadata": {},
   "outputs": [],
   "source": [
    "A = [1.0 4.0 5.0 8.0 1.0; \n",
    "    2.0 -1.0 4.0 3.0 0.0; \n",
    "    7.0 6.0 3.0 -4.0 5.0; \n",
    "    -3.0 4.0 2.0 2.0 2.0; \n",
    "    0.0 -4.0 2.0 1.0 2.0]\n",
    "\n",
    "b = [61.0; 24.0; 37.0; 29.0; 12.0];"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "53124eb8",
   "metadata": {},
   "source": [
    "The code in the following cell converts the general problem $Ax=b$ to the upper triangular equation system $Ux=y$. Note that this function assumes that the pivots are all nonzero. This function will be erroneos if any of the diagonal entries are zero!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "id": "7a7b926a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "convert_to_upper_triangular! (generic function with 1 method)"
      ]
     },
     "execution_count": 67,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function convert_to_upper_triangular!(A,b)\n",
    "    n = size(A,1)\n",
    "    # Upper Triangularization: convert Ax=b to Ux=y\n",
    "    for k in 1:n\n",
    "        for j in k+1:n\n",
    "            # Divide by pivot\n",
    "            A[k,j] = A[k,j] / A[k,k]\n",
    "        end\n",
    "        b[k] = b[k] / A[k,k]\n",
    "        A[k,k] = 1\n",
    "        # Substract lower rows\n",
    "        for i in k+1:n \n",
    "            for j in k+1:n\n",
    "                A[i,j]=A[i,j] - A[i,k] * A[k,j]\n",
    "            end\n",
    "            b[i] = b[i] - A[i,k] * b[k]\n",
    "            A[i,k] = 0\n",
    "        end\n",
    "    end\n",
    "    return A, b #U,y\n",
    "end\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "78ef1849",
   "metadata": {},
   "source": [
    "The function in the following cell solves the upper triangular equation system using backwards substitution. Note that the function alters the input values. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "id": "5a134433",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "solve_upper_triangular! (generic function with 1 method)"
      ]
     },
     "execution_count": 68,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function solve_upper_triangular!(U,y)\n",
    "    n = size(U,1)\n",
    "    for step in reverse(1:n)\n",
    "        if U[step,step] == 0\n",
    "            if y[step] != 0\n",
    "                return \"No solution\"\n",
    "            else\n",
    "                return \"Infinity solutions\"\n",
    "            end\n",
    "        else\n",
    "        # Backwards substitution\n",
    "            y[step] = y[step] / U[step,step]\n",
    "        end\n",
    "        for row in reverse(1:step-1)\n",
    "            y[row] -= U[row,step] * y[step]\n",
    "        end\n",
    "    end\n",
    "    return y \n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "id": "b92332f7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5-element Vector{Float64}:\n",
       " 1.0000000000000009\n",
       " 1.999999999999999\n",
       " 2.9999999999999964\n",
       " 4.000000000000002\n",
       " 5.000000000000005"
      ]
     },
     "execution_count": 69,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "U,y = convert_to_upper_triangular!(A,b)\n",
    "sol = solve_upper_triangular!(U,y)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "962ec4a9",
   "metadata": {},
   "source": [
    "We can test if the obtained solution is correct using `@test`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "id": "0e336c85",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\u001b[32m\u001b[1mTest Passed\u001b[22m\u001b[39m"
      ]
     },
     "execution_count": 70,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using Test\n",
    "@test sol ≈ [1.0; 2.0; 3.0; 4.0; 5.0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "28dff449",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.9.1",
   "language": "julia",
   "name": "julia-1.9"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.9.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}