Commit 4e1c149b authored by James Sutherland's avatar James Sutherland

Initial commit of a few resources

parents
.ipynb_checkpoints
This diff is collapsed.
This diff is collapsed.
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"---\n",
"# Sparse Linear Systems\n",
"\n",
"Contributors:\n",
" * [James C. Sutherland](sutherland.che.utah.edu)\n",
"\n",
"---"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Introduction to Sparse Linear Systems\n",
"We frequently encounter sparse linear systems. For example, when solving the Poisson equation, $$\\nabla^2\\phi=f$$ we end up with a sparse system.\n",
"\n",
"For a one-dimensional solution to the above system if ($f$ is not a function of $\\phi$), we can approximate the Poisson equation as \n",
"$$\n",
"\\frac{\\partial^2 \\phi_{i}}{\\partial x^2} = \\frac{\\phi_{i+1}-2\\phi_{i}+\\phi_{i-1}}{\\Delta x^2} = f_i\n",
"$$\n",
"This results in a large, sparse (lots of zeros in the matrix) system of linear equations to be solved to determine $\\phi$ given $f(x,y)$.\n",
"For Dirichlet boundary conditions ($\\phi$ given on the boundary), this gives a linear system of the form:\n",
"$$\n",
" \\left[\n",
" \\begin{array}{cccccc}\n",
" 1 & 0 & 0 & 0 & 0 & \\cdots & 0 \\\\\n",
" \\frac{1}{\\Delta x^2} & -\\frac{1}{2\\Delta x^2} & \\frac{1}{\\Delta x^2} & 0 & \\cdots & 0 \\\\\n",
" 0 & \\frac{1}{\\Delta x^2} & -\\frac{1}{2\\Delta x^2} & \\frac{1}{\\Delta x^2} & 0 & \\vdots \\\\\n",
" \\vdots & \\ddots & \\ddots & \\ddots & \\ddots & \\vdots \\\\\n",
" 0 & \\cdots & 0 & \\frac{1}{\\Delta x^2} & -\\frac{1}{2\\Delta x^2} & \\frac{1}{\\Delta x^2} \\\\\n",
" 0 & \\cdots & 0 & 0 & 0 & 1 \\\\\n",
" \\end{array}\n",
" \\right]\n",
" \\left(\n",
" \\begin{array}{c}\n",
" \\phi_1 \\\\ \\phi_2 \\\\ \\phi_3 \\\\ \\vdots \\\\ \\phi_{n-1} \\\\ \\phi_n \n",
" \\end{array}\n",
" \\right)\n",
" =\n",
" \\left(\n",
" \\begin{array}{c}\n",
" \\phi_\\mathrm{BC}^\\mathrm{left} \\\\ f_2 \\\\ f_3 \\\\ \\vdots \\\\ f_{n-1} \\\\ \\phi_\\mathrm{BC}^\\mathrm{right }\n",
" \\end{array}\n",
" \\right)\n",
"$$\n",
"\n",
"For large systems, it can be very costly to store the matrix as a dense array. Fortunately, Python can handle [sparse matrices through SciPy](https://docs.scipy.org/doc/scipy/reference/sparse.html)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Key SciPy Functions for Sparse Systems\n",
"There are a few key functions to be aware of:\n",
" * [diags](https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.diags.html#scipy.sparse.diags) and [spdiags](https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.spdiags.html#scipy.sparse.spdiags) which create sparse matrices given diagonal vectors.\n",
" * [Sparse solvers](https://docs.scipy.org/doc/scipy/reference/sparse.linalg.html#solving-linear-problems)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Example\n",
"First we create a function that will set up and solve the linear system:"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def setup_and_solve( nx ):\n",
" import numpy as np\n",
"\n",
" Lx = 2*np.pi # domain extent\n",
" dx = Lx/(nx-1) # grid spacing\n",
"\n",
" x = np.linspace( 0, Lx, nx )\n",
" f = -300*np.exp( -(x-Lx/3)**2 / 0.05 ) - np.exp( -(x-2*Lx/3)**2 / 0.05 )\n",
"\n",
" ld = np.ones(nx-1) / dx**2 # lower diagonal\n",
" md = -2 * np.ones(nx ) / dx**2 # main diagonal \n",
" ud = ld.copy()\n",
" \n",
" # repair boundary condition entries\n",
" ld[nx-2] = 0\n",
" md[0 ] = 1\n",
" md[nx-1] = 1\n",
" ud[0 ] = 0\n",
" f [0 ] = 300\n",
" f [nx-1] = 300\n",
"\n",
" from scipy.sparse import diags\n",
" A = diags(ld,-1) + diags(md,0) + diags(ud,1)\n",
"\n",
"# print(diags(ld,-1))\n",
"# print( 'rank of A = ',np.linalg.matrix_rank(A) )\n",
"# from scipy.sparse.linalg import matrix_rank\n",
"# print( 'rank of A = ',matrix_rank(A) )\n",
"# print(A)\n",
"# print(ld,md,ud,f)\n",
"\n",
" from scipy.sparse.linalg import spsolve\n",
" soln = spsolve( A, f )\n",
" return x, f, soln"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"Next create a function to plot the results:"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"%matplotlib inline\n",
"def plot_solution( x, f, phi ):\n",
" import matplotlib.pyplot as plt\n",
" fig, axes = plt.subplots(nrows=1, ncols=2)\n",
" fig.subplots_adjust(right=2)\n",
" p0 = axes[0].plot( x, phi )\n",
" p1 = axes[1].plot( x, f )\n",
" axes[0].grid()\n",
" axes[1].grid()\n",
" for ax in axes:\n",
" ax.set_xlabel('x');\n",
" axes[0].set_ylabel('$\\phi$')\n",
" axes[1].set_ylabel('f')\n",
" plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally, let's put some interactivity here. Define a \"driver\" function and a widget to interact with. In this case, we will control the number of grid points."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"scrolled": false
},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "8129e04e6d8d47369723127642ed45a2"
}
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"<function __main__.driver>"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"def driver(nx):\n",
" x, f, phi = setup_and_solve(nx)\n",
" plot_solution(x,f,phi)\n",
"\n",
"from ipywidgets import interact\n",
"interact( driver, nx=(5,70,5) )"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(array([ 0. , 1.57079633, 3.14159265, 4.71238898, 6.28318531]),\n",
" array([ 3.00000000e+02, -1.24691053e+00, -8.98302589e-08,\n",
" -4.15636844e-03, 3.00000000e+02]),\n",
" array([ 300. , 302.31003528, 301.54344215, 300.77684879, 300. ]))"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"setup_and_solve(5)\n"
]
}
],
"metadata": {
"hide_input": false,
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"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.6.1"
},
"latex_envs": {
"LaTeX_envs_menu_present": true,
"autocomplete": true,
"bibliofile": "biblio.bib",
"cite_by": "apalike",
"current_citInitial": 1,
"eqLabelWithNumbers": true,
"eqNumInitial": 1,
"hotkeys": {
"equation": "Ctrl-E",
"itemize": "Ctrl-I"
},
"labels_anchors": false,
"latex_user_defs": false,
"report_style_numbering": false,
"user_envs_cfg": false
},
"toc": {
"colors": {
"hover_highlight": "#DAA520",
"navigate_num": "#000000",
"navigate_text": "#333333",
"running_highlight": "#FF0000",
"selected_highlight": "#FFD700",
"sidebar_border": "#EEEEEE",
"wrapper_background": "#FFFFFF"
},
"moveMenuLeft": true,
"nav_menu": {
"height": "66px",
"width": "252px"
},
"navigate_menu": true,
"number_sections": true,
"sideBar": true,
"threshold": 4,
"toc_cell": false,
"toc_section_display": "block",
"toc_window_display": true,
"widenNotebook": false
},
"widgets": {
"state": {
"50738212ca8c4e6facfefb403cdffbe0": {
"views": [
{
"cell_index": 7
}
]
}
},
"version": "1.2.0"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment