数值计算（Python实现）（一）-CSDN博客

数值计算（Python实现）（一）

本篇内容简介：

解线性方程组：高斯消元法和高斯列主元消去法
解线性方程组的迭代方法：雅克比（Jacobi）迭代法与高斯－赛德尔迭代法
拉格朗日插值法
解非线性方程的迭代方法：区间半分法、不动点迭代法和牛顿法

一、解线性方程组

1.1 高斯消元法

def gauss(a,b):
    # a - a is an N * N nonsigular matrix
    # b - b is an N * 1 matrix
    n=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 matrix
    n=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=b
    return 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 iterations
    Output - 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)+f
    for i in range(max1):
        if np.linalg.norm(X-P)>delta:
            P=X
            X=np.dot(G,P)+f
            #print "i=%d"%i,X,np.linalg.norm(X-P)
            I=1
    I=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 iterations
    Output - 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)+f
    for i in range(max1):
        if np.linalg.norm(X-P)>delta:
            P=X
            X=np.dot(G,P)+f
            #print "i=%d"%i,X,np.linalg.norm(X-P)
            I=1
    I=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 estimate
    Output:
    The Lagrange result
    '''
    result=0.0
    for 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 +=t
    return result

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

4.1 区间半分法

def Half(a,b,to1):
    # a - 左边端点
    # b - 右边端点
    # tol - Epsilon 误差
    c=(a+b)/2
    k=1
    n=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 c
        elif f(a)*f(c)<0:
             b=(a+b)/2
        else:
             a=(a+b)/2
        c=(a+b)/2
        k=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=x1
    
    I=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=x1
    
    I=1    
    return (x1, Y(x1), k, I)