Skip to content

Commit

Permalink
Add linear equations serial solution
Browse files Browse the repository at this point in the history
  • Loading branch information
GeliezaK committed Aug 22, 2023
1 parent 5c03b6c commit 904d4ca
Showing 1 changed file with 163 additions and 0 deletions.
163 changes: 163 additions & 0 deletions notebooks/LEQ.ipynb
Original file line number Diff line number Diff line change
@@ -0,0 +1,163 @@
{
"cells": [
{
"cell_type": "markdown",
"id": "599913a3",
"metadata": {},
"source": [
"# Solving Linear Equations\n",
"\n",
"## Serial Algorithm\n",
"First, we construct a linear equations system $Ax=b$."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "fe3b5b02",
"metadata": {},
"outputs": [],
"source": [
"n = 5\n",
"A = rand(-10.0:10.0, (n, n))"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "57c912fb",
"metadata": {},
"outputs": [],
"source": [
"x = rand(-10.0:10.0, n)\n",
"b = A * x"
]
},
{
"cell_type": "markdown",
"id": "4baa0681",
"metadata": {},
"source": [
"The code in the following cell converts the problem $Ax=b$ to the upper triangular equation system $Ux=y$. "
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "11d255e0",
"metadata": {},
"outputs": [],
"source": [
"function convert_to_upper_triangular(A,b)\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": "02c47593",
"metadata": {},
"source": [
"The function in the following cell solves the upper triangular equation system using backwards substitution. "
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "3c19497c",
"metadata": {},
"outputs": [],
"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": null,
"id": "2c24d85e",
"metadata": {},
"outputs": [],
"source": [
"U,y = convert_to_upper_triangular(A,b)\n",
"sol = solve_upper_triangular(U,y)"
]
},
{
"cell_type": "markdown",
"id": "1c356b5a",
"metadata": {},
"source": [
"We can test if the obtained solution is correct using `@test`:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "5c71828d",
"metadata": {},
"outputs": [],
"source": [
"using Test\n",
"@test sol ≈ x"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "fe3c5374",
"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
}

0 comments on commit 904d4ca

Please sign in to comment.