Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

补充一下自己写的python代码 #1

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
Open

补充一下自己写的python代码 #1

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
363 changes: 363 additions & 0 deletions 三角分解,高斯列主元,Jacobi迭代法,Gauss-Seidal迭代法.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,363 @@
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 三角分解法解线性方程组"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"def LUFact(A):# LU Factorization，分解为L,U矩阵\n",
" n = len(A)\n",
" L = np.eye(n)\n",
" U = np.zeros(A.shape)\n",
" U[0,:] = A[0,:]\n",
" L[1:,0] = A[1:,0]/U[0,0]\n",
" for r in range(1,n):\n",
" lt = L[r,:r].reshape(r,1)\n",
" U[r,r:] = A[r,r:] - np.sum(lt*U[:r,r:],axis=0)\n",
"# print('\\n#U%d:'%(r))\n",
"# print(U)\n",
" if r==n-1:\n",
" continue\n",
" ut = U[:r,r].reshape(r,1)\n",
" L[r+1:,r] = (A[r+1:,r] - np.sum(ut*L[r+1:,:r].T,axis=0))/U[r,r] \n",
"# print('\\n#L%d:'%(r))\n",
"# print(L)\n",
" return L,U\n",
" \n",
"def solveLineq(L,U,b):#根据分解后的矩阵计算\n",
" # Ly = b\n",
" rows = len(b)\n",
" y = np.zeros(rows)\n",
" y[0] = b[0]/L[0,0]\n",
" for k in range(1,rows):\n",
" y[k] = (b[k] - np.sum(L[k,:k]*y[:k]))/L[k,k] \n",
" # Ux = y (back substitution) \n",
" x = np.zeros(rows)\n",
" k = rows-1\n",
" x[k] = y[k]/U[k,k] \n",
" for k in range(rows-2,-1,-1):\n",
" x[k] = (y[k] - np.sum(x[k+1:]*U[k,k+1:]))/U[k,k] \n",
" return x\n",
"\n",
"def SanJiao(A, b):#封装三角分解法\n",
" \"\"\"\n",
" 输入：A：numpy矩阵\n",
" b：numpy矩阵，无需转置\n",
" 输出：x：numpy矩阵\n",
" \"\"\"\n",
" L,U = LUFact(A)\n",
" print(\"分解后的L矩阵：\",L)\n",
" print(\"分解后的U矩阵：\",U)\n",
" \n",
" x = solveLineq(L,U,b)\n",
" print(\"计算结果：\",x)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 高斯列主元"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"\n",
"def getInput(A,b):\n",
" matrix_a = np.mat(A, dtype=float)\n",
" matrix_b = np.mat(b)\n",
" # 答案：-2 0 1 1\n",
" return matrix_a, matrix_b\n",
" \n",
" \n",
" \n",
"# 交换\n",
"def swap(mat, num):\n",
" print(num)\n",
" print(\"调换前\")\n",
" print(mat)\n",
" maxid = num\n",
" for j in range(num, mat.shape[0]):\n",
" if mat[j, num] > mat[num, num]:\n",
" maxid = j\n",
" if maxid is not num:\n",
" mat[[maxid,num],:] = mat[[num,maxid],:]\n",
" else:pass\n",
" print(\"调换后\")\n",
" print(maxid)\n",
" print(mat)\n",
" return mat\n",
" \n",
"def SequentialGauss(mat_a):\n",
" for i in range(0, (mat_a.shape[0])-1):\n",
" swap(mat_a, i)\n",
" if mat_a[i, i] == 0:\n",
" print(\"终断运算：\")\n",
" print(mat_a)\n",
" break\n",
" else:\n",
" for j in range(i+1, mat_a.shape[0]):\n",
" mat_a[j:j+1 , :] = mat_a[j:j+1,:] - \\\n",
" (mat_a[j,i]/mat_a[i,i])*mat_a[i, :]\n",
" return mat_a\n",
" \n",
"# 回带过程\n",
"def revert(new_mat):\n",
" # 创建矩阵存放答案 初始化为0\n",
" x = np.mat(np.zeros(new_mat.shape[0], dtype=float))\n",
" number = x.shape[1]-1\n",
" # print(number)\n",
" b = number+1\n",
" x[0,number] = new_mat[number,b]/new_mat[number, number]\n",
" for i in range(number-1,-1,-1):\n",
" try:\n",
" x[0, i] = (new_mat[i,b]-np.sum(np.multiply(new_mat[i,i+1:b],x[0,i+1:b])))/(new_mat[i,i])\n",
" except:\n",
" print(\"错误\")\n",
" print(x)\n",
" \n",
" \n",
"def GSLZY(A, b):\n",
" A,b = getInput(A, b)\n",
" # 合并两个矩阵\n",
" print(\"原矩阵\")\n",
" print(np.hstack((A, b.T)))\n",
" new_mat = SequentialGauss(np.hstack((A, b.T)))\n",
" print(\"三角矩阵\")\n",
" print(new_mat)\n",
" print(\"方程的解\")\n",
" revert(new_mat)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Jacobi迭代法"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"#Jacobi迭代法 输入系数矩阵mx、值矩阵mr、迭代次数n、误差c(以list模拟矩阵 行优先)\n",
" \n",
"def Jacobi(mx,mr,n=100,c=0.0001):\n",
" mr = [[i] for i in mr]\n",
" print(\"mr:\",mr)\n",
" if len(mx) == len(mr): #若mx和mr长度相等则开始迭代 否则方程无解\n",
" x = [] #迭代初值 初始化为单行全0矩阵\n",
" for i in range(len(mr)):\n",
" x.append([0])\n",
" count = 0 #迭代次数计数\n",
" while count < n:\n",
" nx = [] #保存单次迭代后的值的集合\n",
" for i in range(len(x)):\n",
" nxi = mr[i][0]\n",
" for j in range(len(mx[i])):\n",
" if j!=i:\n",
" nxi = nxi+(-mx[i][j])*x[j][0]\n",
" nxi = nxi/mx[i][i]\n",
" nx.append([nxi]) #迭代计算得到的下一个xi值\n",
" lc = [] #存储两次迭代结果之间的误差的集合\n",
" for i in range(len(x)):\n",
" lc.append(abs(x[i][0]-nx[i][0]))\n",
" if max(lc) < c:\n",
" print(\"满足误差要求的结果：\")\n",
" print(nx)\n",
" print(\"迭代次数：\",count)\n",
" return nx #当误差满足要求时 返回计算结果\n",
" x = nx\n",
" count = count + 1\n",
" return False #若达到设定的迭代结果仍不满足精度要求 则方程无解\n",
" else:\n",
" return False"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Gauss-Seidal迭代法"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"\n",
"def G_S(a, b, x=np.array([0,0,0]), g=1e-6): # a为系数矩阵 b增广的一列 x迭代初始值 g计算精度\n",
" x = x.astype(float) #设置x的精度，让x计算中能显示多位小数\n",
" m, n = a.shape\n",
" times = 0 #迭代次数\n",
" if (m < n):\n",
" print(\"There is a 解空间。\") # 保证方程个数大于未知数个数\n",
" else:\n",
" while True:\n",
" for i in range(n):\n",
" s1 = 0\n",
" tempx = x.copy() #记录上一次的迭代答案\n",
" for j in range(n):\n",
" if i != j:\n",
" s1 += x[j] * a[i][j]\n",
" x[i] = (b[i] - s1) / a[i][i]\n",
" times += 1 #迭代次数加一\n",
" gap = max(abs(x - tempx)) #与上一次答案模的差\n",
"\n",
" if gap < g: #精度满足要求，结束\n",
" break\n",
"\n",
" elif times > 10000: #如果迭代超过10000次，结束\n",
" break\n",
" print(\"10000次迭代仍不收敛\")\n",
"\n",
" print(\"迭代次数：\",times)\n",
" print(\"计算结果：\",x)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 数值结果"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"A矩阵：\n",
"[[10 -2 -1]\n",
" [-2 10 -1]\n",
" [-1 -2 5]]\n",
"b矩阵：\n",
"[ 3 15 10]\n",
"===== 三角分解法 =====\n",
"\n",
"#U1:\n",
"[[10. -2. -1. ]\n",
" [ 0. 9.6 -1.2]\n",
" [ 0. 0. 0. ]]\n",
"\n",
"#L1:\n",
"[[ 1. 0. 0. ]\n",
" [-0.2 1. 0. ]\n",
" [-0.1 -0.22916667 1. ]]\n",
"\n",
"#U2:\n",
"[[10. -2. -1. ]\n",
" [ 0. 9.6 -1.2 ]\n",
" [ 0. 0. 4.625]]\n",
"===== 高斯列主元 =====\n",
"原矩阵\n",
"[[10. -2. -1. 3.]\n",
" [-2. 10. -1. 15.]\n",
" [-1. -2. 5. 10.]]\n",
"0\n",
"调换前\n",
"[[10. -2. -1. 3.]\n",
" [-2. 10. -1. 15.]\n",
" [-1. -2. 5. 10.]]\n",
"调换后\n",
"0\n",
"[[10. -2. -1. 3.]\n",
" [-2. 10. -1. 15.]\n",
" [-1. -2. 5. 10.]]\n",
"1\n",
"调换前\n",
"[[10. -2. -1. 3. ]\n",
" [ 0. 9.6 -1.2 15.6]\n",
" [ 0. -2.2 4.9 10.3]]\n",
"调换后\n",
"1\n",
"[[10. -2. -1. 3. ]\n",
" [ 0. 9.6 -1.2 15.6]\n",
" [ 0. -2.2 4.9 10.3]]\n",
"三角矩阵\n",
"[[10. -2. -1. 3. ]\n",
" [ 0. 9.6 -1.2 15.6 ]\n",
" [ 0. 0. 4.625 13.875]]\n",
"方程的解\n",
"[[1. 2. 3.]]\n",
"===== Jacobi迭代法 =====\n",
"mr: [[3], [15], [10]]\n",
"满足误差要求的结果：\n",
"[[0.9999967155648001], [1.9999967163839998], [2.999994593216]]\n",
"迭代次数： 12\n",
"===== Gauss-Seidal迭代法 =====\n",
"迭代次数： 27\n",
"计算结果： [0.99999989 1.99999995 2.99999996]\n"
]
}
],
"source": [
"A = np.array([[10,-2,-1],[-2,10,-1],[-1,-2,5]])\n",
"b = np.array([3,15,10])\n",
"print(\"A矩阵：\")\n",
"print(A)\n",
"print(\"b矩阵：\")\n",
"print(b)\n",
"print(\"=\"*5,\"三角分解法\",\"=\"*5)\n",
"SanJiao(A,b)\n",
"print(\"=\"*5,\"高斯列主元\",\"=\"*5)\n",
"GSLZY(A, b)\n",
"print(\"=\"*5,\"Jacobi迭代法\",\"=\"*5)\n",
"Jacobi(A, b, 999, 0.00001)\n",
"print(\"=\"*5,\"Gauss-Seidal迭代法\",\"=\"*5)\n",
"G_S(A, b)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"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.8.5"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
Loading