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为线性方程组的解

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为线性方程组的解

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

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

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        

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

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)