数值计算方法(高斯消元以及LU分解)
2017-10-03 17:57
309 查看
基础方法没啥好叙述的,书上随处可见,给出了具体的代码实现以及注释,错误之处请留言。谢谢
Mathematical expression of gaussian elimination
elimination-step and get a upper triangular matrix
for k=0 to n-1mki=a(k)ika(k)kk(i=k+1,…,n−1)
ak+1ij=aij(k)−mk∗a(k)kj(i,j=k+1,…,n−1)
bk+1i=bi(k)−mik∗b(k)k(i=k+1,…,n−1)
back-substitution-step
xn−1=b(n−1)n−1an−1(n−1)(n−1)for i=n-2 to 0
xi=b(n)i−∑n−1j=i+1a(n)ijxjanii
Mathematical expression of LU decomposition
1.implement A=LU
u1i=a1i(i=1,2,…,n)li1=ai1u11(i=1,2,…,n)
for k=1 to n-1
uki=aki−∑j=0k−1lkjuji(i=k,k+1,…,n−1)
lik=aik−∑k−1j=0lijujkukk(i=k+1,…,n−1,k!=n)
2.计算LY=b
y1=b1for k=1 to n-1
yk=bk−∑j=0k−1lkjyj(k=2,3,…,n−1)
计算UX=Y
xn−1=bn−1u(n−1)(n−1)for k= n-2 to 0
xk=yk−∑n−1j=k+1ukjxjukk
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-01 2.87032703e-03 6.17144914e-01] [ 3.26644902e-01 1.52705810e+00 8.85942099e-01 ..., 2.34362086e-01 6.16778357e-01 9.49016321e-01] [ 9.50176119e-01 5.56653188e-01 1.91560635e+00 ..., 7.96260777e-02 9.82817114e-01 1.81612851e-01] ..., [ 2.47870273e-01 1.11841381e-01 5.13540776e-01 ..., 1.83527618e+00 2.39285522e-01 7.30797255e-02] [ 9.52966404e-01 1.12326974e-01 8.02396496e-01 ..., 4.95854311e-01 1.50837092e+00 8.47333803e-02] [ 4.37268153e-01 8.63201246e-01 6.80236396e-01 ..., 5.95336949e-02 1.08043656e-01 1.76279378e+00]]"""
gengerate a vector x=(1,2,3,+⋯+100)T
x=np.array(range(1,101)).T
generate a vector b as b=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 matrix for k in range(0,n-1): for i in range(k+1,n): m=A[i][k]/A[k][k]#calculate multiplier for j in range(k,n): A[i][j]=A[i][j]-m*A[k][j] #elimination of ith row b[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-2 while(i>=0): sum2=0 for j in range(i+1,n): sum2+=A[i][j]*newx[j] newx[i]=(b[i]-sum2)/A[i][i] i-=1 return 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=k maxA=A[k][k] for t in range(k+1,n): if(A[t][k]>maxA): maxk=t maxA=A[t][k] for j in range(k,n): tmp=A[k][j] A[k][j]=A[maxk][j] A[maxk][j]=tmp tmp=b[k] b[k]=b[maxk] b[maxk]=tmp #exchange maxk with k for 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 Gauss newx=np.zeros((n,1)) newx[n-1]=(b[n-1]/A[n-1][n-1]) i=n-2 while(i>=0): sum1=0 for j in range(i+1,n): sum1+=A[i][j]*newx[j] newx[i]=(b[i]-sum1)/A[i][i] i-=1 return 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=LU for 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=0 for j in range(0,k): sum3+=L[k][j]*U[j][i] U[k][i]=A[k][i]-sum3 for i in range(k+1,n): sum3=0 for 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=b Y=np.zeros((n,1)) Y[0]=b[0] for k in range(1,n): sum3=0 for j in range(0,k): sum3+=L[k][j]*Y[j] Y[k]=b[k]-sum3 #calculate UX=Y newx=np.zeros((n,1)) newx[n-1]=(Y[n-1]/U[n-1][n-1]) i=n-2 while(i>=0): sum3=0 for j in range(i+1,n): sum3+=U[i][j]*newx[j] newx[i]=(Y[i]-sum3)/U[i][i] i-=1 return 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.T the 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: cbf 17.10.3
相关文章推荐
- 计算方法 -- 解线性方程组直接法(LU分解、列主元高斯消元、追赶法)
- 计算方法实验三 高斯消元
- 整数行列式计算 高斯方法 无精度损失
- 【书单】matlab 科学计算、数值分析以及数学物理问题
- iOS NSDecimalNumber精确数值计算以及小数点后精确保留2位数字
- hdu 3915 Game 求N个数中取若干个数字使得它们的异或值为0的方法数 高斯消元(mod2)
- ArcGis中投影的方法以及计算面积的方法
- 期末了,好久没上了,传一个最近写的矩阵类的原型(目前只有乘法,求行列式以及高斯全主元消元)待完善
- js数值四舍五入的方法以及其中潜在bug的解决方案
- 数值计算方法 数值积分(伪代码 c/c++ python)
- 聊一聊PV和并发、以及计算web服务器的数量的方法【转】
- poj3185(开关问题一般解法,以及高斯消元解法)
- linux shell 时间运算以及时间差计算方法
- Linux Shell 时间运算以及时间差计算方法
- 计算出的多小数位的数值控制小数位的方法
- 采用数值方法计算最速曲线
- linux shell 时间运算以及时间差计算方法
- 用PostgreSQL了解一些统计学术语以及计算方法和表示方法 - 1
- 14章类型信息-之类型转换前先做检查--之使用类字面常量--类名.class--以及动态instanceof(isInstance方法)----递归计数(计算各个类的个数)
- 作业(数值计算方法)