用编程实现计算方法/数值分析的概念知识 （二）

线性方程组的数值解法

这一章是我觉得全书最难的一章，变量多，变量的下标也多，时间一久很难分清。能想出追赶法的人真的是个鬼才！！太巧妙了！！

1. 用列主元消去法求解如下方程组：

   $$ x_1+2x_2+3x_3=14 $$

   $$ 2x_1+5x_2+2x_3=18 $$

   $$ 3x_1+x_2+5x_3=20 $$

   解：

   import numpy as np
   
   def swap(a, b, k, n):           # 找到主元并交换，这仅是一个仅用来交换的函数
    ans = 0
    for i in range(k, n):
        if ans < np.fabs(a[i][k]):    #fabs是绝对值，将a中绝对值最大的找出来
            ans = a[i][k]
            maxn = i
    a[[k, maxn], :] = a[[maxn, k], :]     #交换
    b[k], b[maxn] = b[maxn], b[k]
   
   #主算法
   def gaussin(a, b):
    cnt = 0                         # 统计计算次数
    m, n = a.shape                  # m:a的行数    n:a的列数
   
    if ( m < n ):
        print("方程有多个解")        # 保证方程个数大于未知数个数
        
    else:
        l = np.zeros((m,n))         # 0矩阵
        for i in range(n):
            # 限制条件
            if (a[i][i] == 0):
                print("no answer")  # 对角线不可为0
        # j表示列
        for k in range(n - 1):          # k表示第一层循环，(0，n-1)行
            swap(a, b, k, n)            #在每次计算前，找到最大主元，进行换行
            for i in range(k + 1, n):   # i表示第二层循环,(k+1,n)行,计算该行消元的系数
                l[i][k] = a[i][k] / a[k][k]     #计算
                cnt += 1
                for j in range(m):      # j表示列，对每一列进行运算
                    a[i][j] = a[i][j] - l[i][k] * a[k][j]
                    cnt += 1
                b[i] = b[i] - l[i][k] * b[k]
        # 回代求出方程解
        x = np.zeros(n)
        x[n - 1] = b[n - 1] / a[n - 1][n - 1] #先算最后一位的x解
   
        for i in range(n - 2, -1, -1):      #依次回代倒着算每一个解
            for j in range(i + 1, n):
                b[i] -= a[i][j] * x[j]
                               #自增自减
            x[i] = b[i] / a[i][i]
        for i in range(n):
            print("x" + str(i + 1) + " = ", x[i])
        print("x" " = ", x)
        print("计算次数","=",cnt)
   
   
   if __name__ == '__main__':
    a = np.array([
                    [1, 2, 3], 
                    [2, 5, 2], 
                    [3, 1, 5]
                ])
    b = np.array([14, 18, 20])
    gaussin(a, b)

   插值与逼近

2. 已知$y=f(x)=sinx$的函数表如下，试用拉格朗日插值多项式，求解$sin(0.45)$、$sin(0.55)$、$sin(0.75)$的近似值。

   $x_i$	0.4	0.5	0.6	0.7	0.8
   $f(x_i)$	0.38942	0.47943	0.56464	0.64422	0.71736

   解：

   x_   = [0.4, 0.5, 0.6, 0.7, 0.8]
   f_x_ = [0.38942, 0.47943, 0.56464, 0.64422, 0.71736]
   
   # 插值基函数
   def l0_(x):
    return (x-x_[1])*(x-x_[2])*(x-x_[3])*(x-x_[4])/((x_[0]-x_[1])*(x_[0]-x_[2])*(x_[0]-x_[3])*(x_[0]-x_[4]))
   
   def l1_(x):
    return (x-x_[0])*(x-x_[2])*(x-x_[3])*(x-x_[4])/((x_[1]-x_[0])*(x_[1]-x_[2])*(x_[1]-x_[3])*(x_[1]-x_[4]))
   
   def l2_(x):
    return (x-x_[0])*(x-x_[1])*(x-x_[3])*(x-x_[4])/((x_[2]-x_[0])*(x_[2]-x_[1])*(x_[2]-x_[3])*(x_[2]-x_[4]))
   
   def l3_(x):
    return (x-x_[0])*(x-x_[1])*(x-x_[2])*(x-x_[4])/((x_[3]-x_[0])*(x_[3]-x_[1])*(x_[3]-x_[2])*(x_[3]-x_[4]))
   
   def l4_(x):
    return (x-x_[0])*(x-x_[1])*(x-x_[2])*(x-x_[3])/((x_[4]-x_[0])*(x_[4]-x_[1])*(x_[4]-x_[2])*(x_[4]-x_[3]))
   
   # ============================================================================================
   # 拉格朗日插值多项式
   def L4_(x):
    return f_x_[0]*l0_(x)+f_x_[1]*l1_(x)+f_x_[2]*l2_(x)+f_x_[3]*l3_(x)+f_x_[4]*l4_(x)
   
   # ============================================================================================
   def main():
    print(0.45, "\t", L4_(0.45))
    print(0.55, "\t", L4_(0.55))
    print(0.75, "\t", L4_(0.75))
   
   
   if __name__ == '__main__':
    main()
   

解得近似值为：

2. 利用第1题给出的函数表，构造差商表，并用牛顿差商插值多项式，求解$sin(0.45)$、$sin(0.55)$、$sin(0.75)$的近似值。

解：

x_   = [0.4, 0.5, 0.6, 0.7, 0.8]
f_x_ = [0.38942, 0.47943, 0.56464, 0.64422, 0.71736]
diff_table = [0]

# ==============================================================================
# 生成差商表
def fun(x, f):
    rows = len(x)  # 行数
    cols = len(x) + 1  # 列数
    table = [[0 for col in range(cols)] for row in range(rows)]

    for i in range(rows):
        table[i][0] = x[i]
        table[i][1] = f[i]

    for i in range(1, rows):
        for j in range(2, i + 2):
            table[i][j] = (table[i][j - 1] - table[i - 1][j - 1]) / (table[i][0] - table[i - j + 1][0])

    print("差商表为")
    for col in range(rows):
        print(table[col])
    
    return table

# ==============================================================================
# 生成牛顿差值多项式
def N3_(x):
    diff_table = fun(x_, f_x_)
    return diff_table[0][1]+diff_table[1][2]*(x-x_[0])+diff_table[2][3]*(x-x_[0])*(x-x_[1])+diff_table[3][4]*(x-x_[0])*(x-x_[1])*(x-x_[2])+diff_table[4][5]*(x-x_[0])*(x-x_[1])*(x-x_[2])*(x-x_[3])

# ==============================================================================
def main():
    print(0.45, "\t", N3_(0.45))
    print(0.55, "\t", N3_(0.55))
    print(0.75, "\t", N3_(0.75))

if __name__ == '__main__':
    main()