數(shù)值計(jì)算（Python實(shí)現(xiàn)）（一）

本篇內(nèi)容簡介：

• 解線性方程組：高斯消元法和高斯列主元消去法
• 解線性方程組的迭代方法：雅克比（Jacobi）迭代法與高斯－賽德爾迭代法
• 拉格朗日插值法
• 解非線性方程的迭代方法：區(qū)間半分法、不動(dòng)點(diǎn)迭代法和牛頓法

一、解線性方程組

1.1?高斯消元法

def gauss(a,b):# a - a is an N * N nonsigular matrix# b - b is an N * 1 matrixn=len(b)for i in range(0,n-1):for j in range(i+1,n):if a[j,i]!=0.0:lam=float(a[j,i]/a[i,i]) #乘子a[j,i:n]=a[j,i:n]-lam*a[i,i:n] #第2行起每一行減去第一行的lam倍 b[j]=b[j]-lam*b[i] #常數(shù)項(xiàng)向量也做同樣的操作#print("第%d行"%(j+1),[a,b])for k in range(n-1,-1,-1):b[k]=(b[k]-np.dot(a[k,(k+1):n],b[(k+1):n]))/a[k,k] #可在教科書上找到該公式 result=b return result #result為線性方程組的解

　　一般情形下并不需要自己動(dòng)手編寫一個(gè)計(jì)算線性方程組的高斯消元函數(shù)，直接調(diào)用出現(xiàn)自帶的求解函數(shù)即可。

1.2?高斯列主元消去法

def Pivot_Gauss(a,b):# a - a is an N * N nonsigular matrix# b - b is an N * 1 matrixn=len(b)for i in range(0,n-1):for j in range(i+1,n):if a[j,i]>a[i,i]:a[[i,j],:]=a[[j,i],:]b[[i,j]]=b[[j,i]]#print("換行",[a,b])if a[j,i]!=0.0:lam=float(a[j,i]/a[i,i])a[j,i:n]=a[j,i:n]-lam*a[i,i:n]b[j]=b[j]-lam*b[i]#print("第%d行"%(j+1),[a,b])for k in range(n-1,-1,-1):b[k]=(b[k]-np.dot(a[k,(k+1):n],b[(k+1):n]))/a[k,k] #可在教科書上找到該公式result=breturn result #result為線性方程組的解

二、解線性方程組的迭代方法

2.1?雅克比（Jacobi）迭代法

import numpy as np def Jocobi(A,b,P,delta,max1):'''Input - A is an N * N nonsigular matrix- b is an N * 1 matrix- P is an N * 1 matrix; the initial guess- delta is the tolerance for P- max1 is the maximum number of iterationsOutput - X is an N * 1 matrix; the jacobi approximation to the solution of AX = b- I the indicator of result'''D=np.diag(np.diag(A))E=np.tril(A,-1)F=np.triu(A,1)d=np.linalg.inv(D)G=-np.dot(d,E+F)f=np.dot(d,b)X=np.dot(G,P)+ffor i in range(max1):if np.linalg.norm(X-P)>delta:P=XX=np.dot(G,P)+f#print "i=%d"%i,X,np.linalg.norm(X-P)I=1I=0 return X,I

2.2?高斯－賽德爾迭代法

import numpy as np def Seidel(A,b,P,delta,max1):'''Input - A is an N * N nonsigular matrix- b is an N * 1 matrix- P is an N * 1 matrix; the initial guess- delta is the tolerance for P- max1 is the maximum number of iterationsOutput - X is an N * 1 matrix; the jacobi approximation to the solution of AX = b- I the indicator of result'''D=np.diag(np.diag(A))E=np.tril(A,-1)F=np.triu(A,1)d=np.linalg.inv(D+E)G=-np.dot(d,F)f=np.dot(d,b)X=np.dot(G,P)+ffor i in range(max1):if np.linalg.norm(X-P)>delta:P=XX=np.dot(G,P)+f#print "i=%d"%i,X,np.linalg.norm(X-P)I=1I=0 return X,I

三、拉格朗日插值法

def Lagrange(X, Y, Z):'''Input:X - X[0],X[1],...,X[N]Y - Y[0],Y[1],...,Y[N]Z - The point to estimateOutput:The Lagrange result'''result=0.0for i in range(len(Y)):t=Y[i]for j in range(len(Y)):if i !=j:t*=(Z-X[j])/(X[i]-X[j])result +=treturn result

四、解非線性方程的迭代方法

4.1?區(qū)間半分法

def Half(a,b,to1):# a - 左邊端點(diǎn)# b - 右邊端點(diǎn)# tol - Epsilon 誤差c=(a+b)/2k=1n=1+round((np.log10(b-a)-np.log10(to1)/np.log10(2)))#print ("Total n=%d" % n)while k<=n:#print ("k=%d,a=%f,b=%f,c=%f,f(c)=%f" % (k,a,b,c,f(c)))if f(c)==0: # f()為所需求解的函數(shù)#print("k=%d,c=%c"%(k,c))return celif f(a)*f(c)<0:b=(a+b)/2else:a=(a+b)/2c=(a+b)/2k=k+1#print("x=%f"%c)return c

4.2?不動(dòng)點(diǎn)迭代法

def Picard(x0, to1, N):# x0 - 初始值# to1 - Epsilon 誤差# N - 最大的迭代次數(shù)# I- 最終是否收斂for k in range(N): #print("x(%d)=%.10f"%(k,x0))x1=g(x0) # g(x)為所需求解的函數(shù) if abs(x1-x0)<to1:I=0#print("x(%d)=%.10f"%(k+1,x1))return x1,k+1,I # 最終的x1為所求的解x0=x1I=1 return x1,k+1,I

4.3?牛頓法

def newton_iteration(x0, to1, N, Y, Y1):# x0 - 初始值# to1 - 誤差容限# N - 最大迭代數(shù)# Y - 所求解的函數(shù)# Y1 - 所求解函數(shù)的導(dǎo)數(shù)# I - 最終是否收斂；此函數(shù)下，I=0為收斂for k in range(N): #print("x(%d)=%.10f"%(k,x0))x1=x0-Y(x0)/Y1(x0)if abs(x1-x0)<to1:I=0#print("x(%d)=%.10f"%(k+1,x1))return (x1, Y(x1), k, I) # x1為所求的解x0=x1I=1 return (x1, Y(x1), k, I)

?

轉(zhuǎn)載于:https://www.cnblogs.com/Negan-ZW/p/8407067.html

總結(jié)

以上是生活随笔為你收集整理的数值计算（Python实现）（一）的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。