數值計算（Python實現）（一）

本篇內容簡介：

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

一、解線性方程組

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] #常數項向量也做同樣的操作#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為線性方程組的解

　　一般情形下并不需要自己動手編寫一個計算線性方程組的高斯消元函數，直接調用出現自帶的求解函數即可。

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?區間半分法

def Half(a,b,to1):# a - 左邊端點# b - 右邊端點# 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()為所需求解的函數#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?不動點迭代法

def Picard(x0, to1, N):# x0 - 初始值# to1 - Epsilon 誤差# N - 最大的迭代次數# I- 最終是否收斂for k in range(N): #print("x(%d)=%.10f"%(k,x0))x1=g(x0) # g(x)為所需求解的函數 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 - 最大迭代數# Y - 所求解的函數# Y1 - 所求解函數的導數# I - 最終是否收斂；此函數下，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)

?

轉載于:https://www.cnblogs.com/Negan-ZW/p/8407067.html

總結

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

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。