数值计算方法(高斯消元以及LU分解)
基础方法没啥好叙述的,书上随处可见,给出了具体的代码实现以及注释,错误之处请留言。谢谢
Mathematical expression of gaussian elimination
elimination-step and get a upper triangular matrix
for k=0 to n-1
m_i^k= \frac {a_{ik}^{(k)}} {a_{kk}^{(k)}} (i=k+1,\dots,n-1)
a_{ij}^{k+1}=a_{ij}{(k)}-m_k*a_{kj}^{(k)} (i,j=k+1,\dots,n-1)
b_{i}^{k+1}=b_{i}{(k)}-m_{ik}*b_k^{(k)} (i=k+1,\dots,n-1)
back-substitution-step
x_{n-1}=\frac{b_{n-1}^{(n-1)}}{a_{(n-1)(n-1)}^{n-1}}
for i=n-2 to 0
x_i=\frac{b_i^{(n)}-\sum_{j=i+1}^{n-1}a_{ij}^{(n)}x_j} {a_{ii}^{n}}
Mathematical expression of LU decomposition
1.implement A=LU
u_{1i}=a_{1i} (i=1,2,\dots,n)
l_{i1}=\frac{a_{i1}}{u_{11}} (i=1,2,\dots,n)
for k=1 to n-1
u_{ki}=a_{ki}-\sum_{j=0}^{k-1}l_{kj}u_{ji} (i=k,k+1,\dots,n-1)
l_{ik}=\frac{a_{ik}-\sum_{j=0}^{k-1}l_{ij}u_{jk}}{u_kk} (i=k+1,\dots,n-1,k!=n)
2.计算LY=b
y1=b1
for k=1 to n-1
y_k=b_k- \sum_{j=0}^{k-1}l_{kj}y_j (k=2,3,\dots,n-1)
计算UX=Y
x_{n-1}=\frac{b_{n-1}}{u_{(n-1)(n-1)}}
for k= n-2 to 0
x_k=\frac{y_k-\sum_{j=k+1}^{n-1}u_{kj}x_j} {u_{kk}}
import numpy as np #a package to calculate
from scipy.sparse import identity# to get a identity matrix
import time#calculate time of diffient method
create a random matrix M with dimension 100*100 ,generate A by adding an identify matrix
np.random.seed(1) #set seed to make s unchanged
s=np.random.rand(100,100)
A=s+identity(100).toarray()
#add a identity to make sure has a inverse matrix
print A
"""[[ 1.41702200e+00 7.20324493e-01 1.14374817e-04 ..., 5.73679487e-012.87032703e-03 6.17144914e-01][ 3.26644902e-01 1.52705810e+00 8.85942099e-01 ..., 2.34362086e-016.16778357e-01 9.49016321e-01][ 9.50176119e-01 5.56653188e-01 1.91560635e+00 ..., 7.96260777e-029.82817114e-01 1.81612851e-01]..., [ 2.47870273e-01 1.11841381e-01 5.13540776e-01 ..., 1.83527618e+002.39285522e-01 7.30797255e-02][ 9.52966404e-01 1.12326974e-01 8.02396496e-01 ..., 4.95854311e-011.50837092e+00 8.47333803e-02][ 4.37268153e-01 8.63201246e-01 6.80236396e-01 ..., 5.95336949e-021.08043656e-01 1.76279378e+00]]"""
gengerate a vector x=(1,2,3,+⋯+100)Tx=(1,2,3,+ \cdots + 100)^T
x=np.array(range(1,101)).T
generate a vector b as b=AXb=AX
b=np.dot(A,x)
#get result by calling function ,also you can calculate by yourself
#for i in range(100):
# b[i]=0
# for j in range(100):
# b[i]+=A[i][j]*x[j][0]
calculate x
#Gauss-Jordan elimination
def Gauss(n,A,b):#elimination-step and get a upper triangular matrixfor k in range(0,n-1):for i in range(k+1,n):m=A[i][k]/A[k][k]#calculate multiplierfor j in range(k,n):A[i][j]=A[i][j]-m*A[k][j] #elimination of ith rowb[i]=b[i]-m*b[k]"""back-substitution-step,easy to implement according to arithmetic expression"""newx=np.zeros((n,1))newx[n-1]=(b[n-1]/A[n-1][n-1])i=n-2while(i>=0):sum2=0for j in range(i+1,n):sum2+=A[i][j]*newx[j]newx[i]=(b[i]-sum2)/A[i][i]i-=1return newx"""this method is based on Gauss,just exchange kth row with maxkth
row before calculate m to avoid some mistakes caused by float"""
def Gauss_exchange(n,A,b):for k in range(0,n-1):maxk=kmaxA=A[k][k]for t in range(k+1,n):if(A[t][k]>maxA):maxk=tmaxA=A[t][k]for j in range(k,n):tmp=A[k][j]A[k][j]=A[maxk][j]A[maxk][j]=tmptmp=b[k] b[k]=b[maxk]b[maxk]=tmp #exchange maxk with kfor i in range(k+1,n):m=A[i][k]/A[k][k]for j in range(k,n):A[i][j]=A[i][j]-m*A[k][j]b[i]=b[i]-m*b[k]#same as Gaussnewx=np.zeros((n,1))newx[n-1]=(b[n-1]/A[n-1][n-1])i=n-2while(i>=0):sum1=0for j in range(i+1,n):sum1+=A[i][j]*newx[j]newx[i]=(b[i]-sum1)/A[i][i]i-=1return newx
def dolittle(n,A,b):L=np.zeros(A.shape,dtype="float")+identity(100).toarray()U=np.zeros(A.shape,dtype="float")#implement A=LUfor i in range(0,n):U[0][i]=A[0][i]for i in range(0,n):L[i][0]=A[i][0]/U[0][0]for k in range(1,n):for i in range(k,n):sum3=0for j in range(0,k):sum3+=L[k][j]*U[j][i]U[k][i]=A[k][i]-sum3for i in range(k+1,n):sum3=0for j in range(0,k):sum3=sum3+L[i][j]*U[j][k]L[i][k]=(A[i][k]-sum3)/U[k][k]#calculate LY=bY=np.zeros((n,1))Y[0]=b[0]for k in range(1,n):sum3=0for j in range(0,k):sum3+=L[k][j]*Y[j]Y[k]=b[k]-sum3#calculate UX=Ynewx=np.zeros((n,1))newx[n-1]=(Y[n-1]/U[n-1][n-1])i=n-2while(i>=0):sum3=0for j in range(i+1,n):sum3+=U[i][j]*newx[j]newx[i]=(Y[i]-sum3)/U[i][i]i-=1return newx"""copy() function is to make sure do not change A
and pass correct parameters in follow functions"""
t0=time.time() #start time
newx2=Gauss(100,A.copy(),b.copy())
print "the time of Gaussian elimination is ",time.time()-t0
t0=time.time()
newx1=Gauss_exchange(100,A.copy(),b.copy())
print "the time of column principal element elimination is ",time.time()-t0
t0=time.time()
newx3=dolittle(100,A.copy(),b.copy())
print "the time of LU decomposition method is ",time.time()-t0
print "first ",newx1.T, " second",newx2.T,"\n third",newx3.Tthe time of Gaussian elimination is 0.805999994278
the time of column principal element elimination is 0.884999990463
the time of LU decomposition method is 0.444000005722""" first [[ 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14.15. 16. 17. 18. 19. 20. 21. 22. 23. 24. 25. 26. 27. 28.29. 30. 31. 32. 33. 34. 35. 36. 37. 38. 39. 40. 41. 42.43. 44. 45. 46. 47. 48. 49. 50. 51. 52. 53. 54. 55. 56.57. 58. 59. 60. 61. 62. 63. 64. 65. 66. 67. 68. 69. 70.71. 72. 73. 74. 75. 76. 77. 78. 79. 80. 81. 82. 83. 84.85. 86. 87. 88. 89. 90. 91. 92. 93. 94. 95. 96. 97. 98.99. 100.]] second [[ 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14.15. 16. 17. 18. 19. 20. 21. 22. 23. 24. 25. 26. 27. 28.29. 30. 31. 32. 33. 34. 35. 36. 37. 38. 39. 40. 41. 42.43. 44. 45. 46. 47. 48. 49. 50. 51. 52. 53. 54. 55. 56.57. 58. 59. 60. 61. 62. 63. 64. 65. 66. 67. 68. 69. 70.71. 72. 73. 74. 75. 76. 77. 78. 79. 80. 81. 82. 83. 84.85. 86. 87. 88. 89. 90. 91. 92. 93. 94. 95. 96. 97. 98.99. 100.]] third [[ 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14.15. 16. 17. 18. 19. 20. 21. 22. 23. 24. 25. 26. 27. 28.29. 30. 31. 32. 33. 34. 35. 36. 37. 38. 39. 40. 41. 42.43. 44. 45. 46. 47. 48. 49. 50. 51. 52. 53. 54. 55. 56.57. 58. 59. 60. 61. 62. 63. 64. 65. 66. 67. 68. 69. 70.71. 72. 73. 74. 75. 76. 77. 78. 79. 80. 81. 82. 83. 84.85. 86. 87. 88. 89. 90. 91. 92. 93. 94. 95. 96. 97. 98.99. 100.]]
conclusion:the data is small ,It’s hard to compare performance
author: cbf17.10.3
数值计算方法(高斯消元以及LU分解)相关推荐
- 计算矩阵的逆和行列式的值(高斯消元+LU分解)
计算矩阵的逆 选主元的高斯消元法 朴素的高斯消元法是将矩阵A和单位矩阵放在一起,通过行操作(或者列操作)将A变为单位矩阵,这个时候单位矩阵就是矩阵A的逆矩阵.从上到下将A变为上三角矩阵的复杂度为O(n ...
- 【算法设计与分析基础】高斯消元
参考<算法导论>第七部分 算法问题选编中的第28章 矩阵运算. #include<iostream> using namespace std;/*** 一些矩阵例子*//* d ...
- 并行程序设计实验——高斯消元
并行程序设计实验--高斯消元 一.问题描述 熟悉高斯消元法解线性方程组的过程,然后实现SSE算法编程.过程中,自行构造合适的线性方程组,并选取至少2个角度,讨论不同算法策略对性能的影响. 可选角度包括 ...
- 【bzoj3601】一个人的数论 莫比乌斯反演+莫比乌斯函数性质+高斯消元
Description Sol 这题好难啊QAQ 反正不看题解我对自然数幂求和那里是一点思路都没有qwq 先推出一个可做一点的式子: \(f(n)=\sum_{k=1}^{n}[(n,k)=1]k^d ...
- hihoCoder#1196 : 高斯消元·二(开关灯问题)
传送门 高斯消元解异或方程组 小Ho在游戏板上忙碌了30分钟,任然没有办法完成,于是他只好求助于小Hi. 小Ho:小Hi,这次又该怎么办呢? 小Hi:让我们来分析一下吧. 首先对于每一个格子的状态,可 ...
- 【SDOI2017】硬币游戏【KMP】【概率期望】【高斯消元】
题意:给 nnn 个长度为 mmm 的 01 串,一个 01 串初始为空,不断随机一个字符加在后面,当出现给定的 nnn 个串中的一个时停止.分别求在 nnn 个串处停止的概率. 考场思路历程: 显然 ...
- Matrix 高斯消元Gaussian elimination 中的complete pivoting和partial pivoting
首先科普下Pivoting的含义 一般翻译为"主元",在对矩阵做某种算法时,首先进行的部分元素.在线性规划的单纯形法中常见. wiki的解释如下: Pivot element (t ...
- [算法] 高斯消元详解
原文链接(强烈建议安利) 0.前置知识 知道如何解三元一次方程组 有手,有脑子 1.答案的表示与存储 先解一个方程组: 2x+3y+5z=31 x-4y -z=-6 4x+2y-5z=9 我们把这个方 ...
- 解线性方程组——高斯消元の板子
ATP记得它在很久以前看过一点点高斯消元的东西然后做过一点点题目..但是当时实在是太zz了所以本来就没有很懂这个东西现在更是忘得差不多了.. 所以现在就当重新学一遍了QwQ 一点口胡的解释 高斯消元. ...
最新文章
- 鸿蒙系统搁置,华为:我们将坚定的支持安卓生态,鸿蒙系统没有明确上市时间...
- 如何自学python数据分析-Python学习干货 |如何用Python进行数据分析?
- [LeetCode]#13 3sum
- Php7实现文件下载,PHP7 SFTP下载文件并重命名该下载文件
- how to request a curl operation from alibaba cloud
- mysql如何开启远程链接_mysql怎么开启远程连接
- 如何查询服务器是否安装系统时间,如何查看系统当前的NTP配置?
- NetBeans 8.0的五个新性能提示
- 清华大学-曾鸣-《ARM微控制器与嵌入式系统》I2C总线(一)
- 重磅!央行启动企业信息联网核查系统
- 100行代码搞定抖音短视频App,终于可以和美女合唱了。
- 基于matlab的车牌识别
- 用HTML和CSS制作简单的静态网页
- 手把手BC26模组OpenCPU开发之旅-1.简介
- php写浏览记录,php 浏览历史记录的实现方法
- android q mix3,小米MIX3成首款适配Android Q的5G手机
- iOS设备的越狱方法
- 北邮C++——破解简单密码
- WPF快速学习--一布局
- 华为应用市场名称问题
热门文章
- Windows动态定义模板类对象
- [有限元] 四结点三角形单元和五结点三角形单元的形函数
- 计算机教师自媒体方向,教师和自媒体,我该选择哪个深耕?
- layui响应式:隐藏与显示(class 类名后缀)
- phpcmsV9支付: 支付宝支付配置 (资源汇总)
- css文字溢出部分在另一个div显示(代码篇)
- Vue.js项目去除url中的#/ - 解决篇
- Bootstrap警告框、弹出提示层、模态框的js插件效果总结
- python怎么用第三方库_python中第三方库的下载方法
- deepnode处理过的图片_教你用PS快速修复图片脏乱和瑕疵,快来一起学习吧!