數值計算（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)