diff --git a/notebooks/LEQ.ipynb b/notebooks/LEQ.ipynb new file mode 100644 index 0000000..a7e4fca --- /dev/null +++ b/notebooks/LEQ.ipynb @@ -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 +}