{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"toc": "true"
},
"source": [
"# Table of Contents\n",
" "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 3x3行列のLU分解と解"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"array([[ 1., 4., 3.],\n",
" [ 1., -2., 1.],\n",
" [ 2., -2., -1.]])\n",
"array([[ 11.],\n",
" [ 11.],\n",
" [ 11.]])\n"
]
}
],
"source": [
"from pprint import pprint\n",
"import scipy.linalg as linalg # SciPy Linear Algebra Library\n",
"import numpy as np\n",
"\n",
"# A = np.array([[1,1,-2,-4],[1,-2,1,5],[2,-2,-1,2]])\n",
"A = np.array([[1.0,4,3],[1,-2,1],[2,-2,-1]])\n",
"b = np.array([[11.0],[11],[11]])\n",
"pprint(A)\n",
"pprint(b)"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[ 0.18181818 -0.09090909 0.45454545]\n",
" [ 0.13636364 -0.31818182 0.09090909]\n",
" [ 0.09090909 0.45454545 -0.27272727]]\n"
]
},
{
"data": {
"text/plain": [
"array([[ 6.],\n",
" [-1.],\n",
" [ 3.]])"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"inv_A = linalg.inv(A)\n",
"print(inv_A)\n",
"np.dot(inv_A,b)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"array([[ 1., 4., 3., 11.],\n",
" [ 1., -2., 1., 11.],\n",
" [ 2., -2., -1., 11.]])\n",
"array([[ 1., 4., 3., 11.],\n",
" [ 1., -2., 1., 11.],\n",
" [ 2., -2., -1., 11.]])\n"
]
}
],
"source": [
"e_A = np.hstack((A, b))\n",
"pprint(e_A)\n",
"pprint(np.column_stack((A,b)))"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"P:\n",
"array([[ 0., 1., 0.],\n",
" [ 0., 0., 1.],\n",
" [ 1., 0., 0.]])\n",
"L:\n",
"array([[ 1. , 0. , 0. ],\n",
" [ 0.5, 1. , 0. ],\n",
" [ 0.5, -0.2, 1. ]])\n",
"U:\n",
"array([[ 2. , -2. , -1. , 11. ],\n",
" [ 0. , 5. , 3.5, 5.5],\n",
" [ 0. , 0. , 2.2, 6.6]])\n"
]
}
],
"source": [
"P, L, U = linalg.lu(e_A)\n",
"\n",
"print(\"P:\")\n",
"pprint(P)\n",
"\n",
"print(\"L:\")\n",
"pprint(L)\n",
"\n",
"print(\"U:\")\n",
"pprint(U)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([[ 2., -2., -1., 11.],\n",
" [ 1., 4., 3., 11.],\n",
" [ 1., -2., 1., 11.]])"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"np.dot(L,U)"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"array([[ 1. , 4. , 3. ],\n",
" [ 0. , -6. , -2. ],\n",
" [ 0. , 0. , -3.66666667]])\n",
"array([[ 11.],\n",
" [ 0.],\n",
" [-11.]])\n"
]
}
],
"source": [
"from pprint import pprint\n",
"import scipy.linalg as linalg # SciPy Linear Algebra Library\n",
"import numpy as np\n",
"\n",
"\n",
"A = np.array([[1.0,4,3],[1,-2,1],[2,-2,-1]])\n",
"A0 = np.array([[1.0,4,3],[1,-2,1],[2,-2,-1]])\n",
"b = np.array([[11.0],[11],[11]])\n",
"\n",
"\n",
"n = 3\n",
"L = np.identity(n)\n",
"T = []\n",
"for i in range(n): #i行目\n",
" T.append(np.identity(n))\n",
" for j in range(i+1, n):\n",
" am = A[j,i]/A[i,i] #i行の要素を使って，i+1行目の先頭を消す係数を求める\n",
" T[i][j,i]=-am #i番目の消去行列に要素を入れる\n",
" L[j,i]=am #LTMの要素\n",
" for k in range(n):\n",
" A[j,k] -= am*A[i,k] #もとの行列をUTMにしていく\n",
" b[j] -= b[i]*am #数値ベクトルも操作\n",
"pprint(A)\n",
"pprint(b)"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([[ 1., 4., 3.],\n",
" [ 1., -2., 1.],\n",
" [ 2., -2., -1.]])"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"np.dot(L,A)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"array([[ 1. , 0. , 0. ],\n",
" [ 1. , 1. , 0. ],\n",
" [ 2. , 1.66666667, 1. ]])\n"
]
}
],
"source": [
"pprint(L)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# pivot付きのLU分解"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"array([[ 0., 0., 0., 1.],\n",
" [ 0., 0., 1., 0.],\n",
" [ 0., 1., 0., 0.],\n",
" [ 1., 0., 0., 0.]])\n",
"array([[ 1. , 0. , 0. , 0. ],\n",
" [ 0.2 , 1. , 0. , 0. ],\n",
" [ 0.6 , -0.076923, 1. , 0. ],\n",
" [ 0.6 , -0.076923, 0.75 , 1. ]])\n",
"array([[ 5. , 3. , -2. , 5. , 2. ],\n",
" [ 0. , -2.6 , -2.6 , -0. , -9.4 ],\n",
" [ 0. , 0. , 4. , -2. , 0.076923],\n",
" [ 0. , 0. , 0. , -0.5 , -7.980769]])\n",
"array([[ 5., 3., -2., 5.],\n",
" [ 1., -2., -3., 1.],\n",
" [ 3., 2., 3., 1.],\n",
" [ 3., 2., 2., 1.]])\n",
"array([[ 5., 3., -2., 5., 2.],\n",
" [ 1., -2., -3., 1., -9.],\n",
" [ 3., 2., 3., 1., 2.],\n",
" [ 3., 2., 2., 1., -6.]])\n"
]
}
],
"source": [
"from pprint import pprint\n",
"import scipy.linalg as linalg # SciPy Linear Algebra Library\n",
"import numpy as np\n",
"\n",
"A=np.array([[3,2,2,1],[3,2,3,1],[1,-2,-3,1],[5,3,-2,5]])\n",
"b=np.array([-6,2,-9,2])\n",
"ex_A = np.column_stack((A,b))\n",
"P,L,U = linalg.lu(ex_A)\n",
"pprint(P)\n",
"pprint(L)\n",
"pprint(U)\n",
"\n",
"pprint(np.dot(P,A))\n",
"pprint(np.dot(L,U))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Jacobi vs Gauss-Seidel\n",
"\n",
"```\n",
"[-1.2 0.666667 0.777778 0.6 ]\n",
"[-1.608889 0.607407 0.562963 0.751111]\n",
"[-1.584296 0.764938 0.54749 0.782519]\n",
"[-1.618989 0.751429 0.518705 0.676892]\n",
"[-1.589405 0.807797 0.506116 0.680422]\n",
"```"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[-1.2 1.066667 0.407407 0.362963]\n",
"[-1.567407 0.932346 0.436763 0.528779]\n",
"[-1.579578 0.871345 0.46739 0.580064]\n",
"[-1.58376 0.845435 0.478382 0.600844]\n",
"[-1.584932 0.835236 0.482827 0.608976]\n"
]
}
],
"source": [
"# Jacobi iterative method for solving Ax=b by bob\n",
"import numpy as np\n",
"np.set_printoptions(precision=6, suppress=True)\n",
"\n",
"A=np.array([[5,1,1,1],[1,3,1,1],[1,-2,-9,1],[1,3,-2,5]])\n",
"b=np.array([-6,2,-7,3])\n",
"n=4\n",
"x0=np.zeros(n)\n",
"x1=np.zeros(n)\n",
"\n",
"for iter in range(0, 5):\n",
" for i in range(0, n):\n",
" x1[i]=b[i]\n",
" for j in range(0, n):\n",
" x1[i]=x1[i]-A[i][j]*x0[j]\n",
" x1[i]=x1[i]+A[i][i]*x0[i]\n",
" x1[i]=x1[i]/A[i][i]\n",
" x0[i]=x1[i]\n",
" print(x0)"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[-1.2 1.066667 0.407407 0.362963]\n",
"[-1.567407 0.932346 0.436763 0.528779]\n",
"[-1.579578 0.871345 0.46739 0.580064]\n",
"[-1.58376 0.845435 0.478382 0.600844]\n",
"[-1.584932 0.835236 0.482827 0.608976]\n"
]
}
],
"source": [
"# Jacobi iterative method for solving Ax=b by bob\n",
"import numpy as np\n",
"np.set_printoptions(precision=6, suppress=True)\n",
"\n",
"A=np.array([[5,1,1,1],[1,3,1,1],[1,-2,-9,1],[1,3,-2,5]])\n",
"b=np.array([-6,2,-7,3])\n",
"n=4\n",
"x0=np.zeros(n)\n",
"\n",
"for iter in range(0, 5):\n",
" for i in range(0, n):\n",
" x1 = b[i]\n",
" for j in range(0, n):\n",
" x1 -= A[i][j]*x0[j]\n",
" x1 += A[i][i]*x0[i]\n",
" x1 /= A[i][i]\n",
" x0[i] = x1\n",
" print(x0)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
}
],
"metadata": {
"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": "30px",
"width": "252px"
},
"navigate_menu": true,
"number_sections": true,
"sideBar": true,
"threshold": 4,
"toc_cell": true,
"toc_section_display": "block",
"toc_window_display": true,
"widenNotebook": false
}
},
"nbformat": 4,
"nbformat_minor": 2
}