[数值分析] python代码模板

数值分析课，作业中算法代码模板。已结课，之后随缘更新

仅供个人使用，转载请注明源地址。

更新日志：
2020.05.05 添加椭圆坐标生成函数、非线性方程求解方法
2020.06.15 添加A = LU分解和PA = LU分解方法

import numpy as np
import matplotlib.pyplot as plt
import math

# 生成椭圆坐标
def get_ellipse(e_x, e_y, a = 1, b = 1, e_angle = 0.0):
    angles_circle = np.arange(0, 2 * np.pi, 0.01)
    x = []
    y = []
    for angles in angles_circle:
        or_x = a * np.cos(angles)
        or_y = b * np.sin(angles)
        length_or = np.sqrt(or_x * or_x + or_y * or_y)
        or_theta = math.atan2(or_y, or_x)
        new_theta = or_theta + e_angle/180*np.pi
        new_x = e_x + length_or * np.cos(new_theta)
        new_y = e_y + length_or * np.sin(new_theta)
        x.append(new_x)
        y.append(new_y)
    return x, y

# 秦九韶算法(Horner Algorithm)模板
def nest(d, c, x, b = 0):
    '''秦九韶算法(Horner Algorithm)模板
    '''
    if b == 0:
        b = np.zeros((d, 1))
    y = c[d]
    for i in range(d - 1, -1, -1):
        y = y * (x - b[i]) + c[i]
    return y

# 对分法(二分法)模板
def dichotomy(F, lef, rig, eps = 1e-8, round = 0, debug = False):
    '''对分法(二分法)模板
    '''
    i = 0
    while (rig - lef) > eps:
        if round != 0 and i > round:
            break
        i += 1
        mid = (lef + rig) / 2
        val = F(mid)
        if debug:
            print("round", i, "\tvalue", val, "\twidth", rig - lef)
        if np.sign(val) != np.sign(F(lef)):
            rig = mid
        else:
            lef = mid
    return (lef + rig) / 2

# 不动点迭代模板
def fpi(F, init_, eps = 1e-8, round = 20, debug = False):
    '''不动点迭代(FPI)模板
    '''
    x = F(init_)
    pre = init_
    diff = np.abs(x - pre)
    prediff = diff
    if prediff == 0:
        S = 0
    else:S = diff/prediff
    if debug:
        print("case\tvalue\t\t\tdiff\t\t\tS")
        print("%d\t%.15f\t%.15f\t%.15f" % (0, x, diff, S))
    for i in range(1, round):
        pre = x
        x = F(x)
        prediff = diff
        diff = np.abs(x - pre)
        if prediff == 0:
            S = 0
        else:S = diff/prediff
        if debug:
            print("%d\t%.15f\t%.15f\t%.15f" % (i, x, diff, S))
        if  diff < eps or math.isinf(x) or math.isnan(x):
            break
    return x, S

# 牛顿迭代模板
def newton(f, df, init_ = 0, rounds = 20, total = 0.5e-8, m = 1.0, debug = False, eps = 0.5e-16, trace = False):
    '''牛顿迭代模板
    '''
    theta = 1.0
    x = init_
    ei_1 = ei = 1.0

    if debug:
        print("round\t x\t\t\t ei+1/ei\t\t ei+1/ei^2")

    if trace:
        T = np.array([x])

    for i in range(rounds):
        fdx = df(x)
        delta = m * f(x) / fdx
        x = x - delta
        if trace:
            T = np.append(T, [x])
        ei_1 = np.abs(-delta)
        if debug:
            print(i+1, "\t", x, "\t", ei_1/ei, "\t", ei_1/(ei**2))
        if total != 0.0:
            if (np.abs(delta)/max(np.abs(x), theta) < total) or (np.abs(f(x)) < eps):
                break
        ei = ei_1
    if trace:
        return x, T
    return x

# 割线法求根
def secant(f, init_ = [1, 2], rounds = 20, total = 0.5e-8, m = 1.0, eps = 0.5e-16, debug = False):
    '''割线法求根
    '''
    x0, x1 = init_[0], init_[1]
    theta = 1.0

    if debug:
        print('round\t','x')

    for i in range(rounds):
        f0, f1 = f(x0), f(x1)
        dx = - f1/(f1-f0)*(x1 - x0)
        x0 = x1 + dx
        abserr = np.abs(dx)
        relerr = abserr/max(abs(x0), theta)

        if debug:
            print(i+1, '\t', x0)

        if np.abs(f(x0)) < eps or relerr < total:
            break
        x1, x0 = x0, x1
    return x0

# 试位法模板
def FP(f, x = [1.0, 2.0], rounds = 50, total = 0.5e-8, m = 1.0, eps = 0.5e-16, debug = False):
    '''试位法模板
    '''
    fx = [f(x[0]), f(x[1])]
    mid = 0
    if np.sign(fx[0]) == np.sign(fx[1]):
        return

    if debug:
        print('rounds \t\t a\t\t b \t\t c')

    for i in range(rounds):
        fx = [f(x[0]), f(x[1])]
        mid = x[0] - fx[0]/(fx[0]-fx[1])*(x[0]-x[1])
        fm = f(mid)
        
        if debug:
            print(i+1,'\t', x[0], '\t', x[1], '\t', mid)
        
        if np.abs(fm) < eps or np.abs(x[1]-x[0]) < total:
            break
        if np.sign(fm) == np.sign(fx[0]):
            x[0] = mid
        else:
            x[1] = mid
    return mid

# 反二插值法模板
def IOI(f, x = [1, 2, 3], rounds = 50, total = 0.5e-8, m = 1.0, eps = 0.5e-16, debug = False):
    '''反二插值法模板
    '''
    theta = 1
    x = np.array(x)

    if debug:
        print('rounds \t x')

    for i in range(rounds):
        fx = f(x)

        if (fx == np.zeros_like(fx)).any():
            return fx[fx == np.zeros_like(fx)]

        q, r, s = fx[0]/fx[1], fx[2]/fx[1], fx[2]/fx[0]
        delta = - (r*(r-q)*(x[2]-x[1])+(1-r)*s*(x[2]-x[0]))/((q-1)*(r-1)*(s-1))

        x = np.array([x[1], x[2], x[2] + delta])
        abserr = np.abs(delta)
        relerr = abserr/max(x[2], theta)

        if debug:
            print(i+1, '\t', x[2])

        if np.abs(f(x[2])) < eps or relerr < total:
            break
    return x[2]

# 高斯消元模板
def Gauss(A, B):
    '''高斯消元模板
    '''
    A = np.mat(A)
    B = np.mat(B)
    if A.shape[0] != A.shape[1]:
        print('A must be n*n matrix')
        return
    n = A.shape[0]
    x = np.ndarray((n, 1))
    for k in range(n - 1):
        mk = A[k+1:n, k]/A[k, k]
        A[k+1:n, k:n] = A[k+1:n, k:n] - mk@(A[k, k:n].reshape(1,-1))
        B[k+1:n] = B[k+1:n] - mk@B[k]
    for i in range(n-1, -1, -1):
        x[i] = B[i]/A[i, i]
        if(i != 0):
            B[:i, 0] = B[:i, 0] - float(x[i]) * A[:i, i]
    return x

# 希尔伯特矩阵生成模板
def Hilbert(n):
    A = np.ndarray((n, n))
    for i in range(n):
        for j in range(n):
            A[i, j] = 1/(i + j +1)
    return A

# jacobi方法求方程组
def jacobi(A, b, N=25, x=None, eps=1e-6):
    if x is None:
        pre = x = np.zeros(len(A[0]))
    D = np.diag(A)
    R = A - np.diagflat(D)

    for i in range(N):
        x = (b - np.dot(R,x)) / D
        ferr = np.linalg.norm((x - pre), ord = np.inf)
        pre = x
        if ferr < eps :
            print("After %d steps" % (i))
            berr = np.linalg.norm(A@x - b, ord = np.inf)
            print("backward_error =", berr)
            print("forward_error =", ferr)
            break
        
    return x

# SOR/Gauss-Seidel方法求线性方程组
def SOR(A, b, w=1.0, N=25, x=None, eps=1e-6):
    if x is None:
        pre = x = np.zeros(len(A[0]))

    L, D, U = np.tril(A, -1), np.diagflat(np.diag(A)), np.triu(A, 1)

    for i in range(N):
        x = np.linalg.inv(w*L + D)@((1.0 - w)*D@x - w*U@x) + w*np.linalg.inv(D + w*L)@b
        ferr = np.linalg.norm((x - pre), ord = np.inf)
        pre = x

        if ferr < eps:
            print("After %d steps" % (i))
            berr = np.linalg.norm(A@x - b, ord = np.inf)
            print("backward_error =", berr)
            print("forward_error =", ferr)
            break

    return x

# 共轭梯度法求线性方程组
def conj_grad(A, b, N=25, x=None, eps=1e-6, debug=False):
    if x is None:
        pre = x = np.zeros(len(A[0]))
    b0 = b
    d = r = b
    for i in range(N):
        temp = r.T@r
        a = temp/(d.T@A@d)
        x = x + a*d
        r = r - a*A@d
        b = r.T@r/temp
        d = r + b*d

        ferr = np.linalg.norm(x - pre, ord = np.inf)
        berr = np.linalg.norm(r, ord = np.inf)
        pre = x
        if debug:
            print(i, '\t', ferr, '\t', berr)
        
        if ferr < eps:
            print("After %d steps" % (i))
            print("backward_error/norm(r, inf) =", berr)
            print("forward_error =", ferr)
            break
    
    return x

# 牛顿法求非线性方程组
def n_newton(f, df, x, N=25, eps=1e-6):
    for i in range(N):
        A, b = df(x), -f(x)
        b = b.reshape(-1, 1)
        s = Comp.Gauss(A, b)
        s = s.ravel()
        x = x + s

        ferr = np.linalg.norm(s, ord = np.inf)
        berr = np.linalg.norm(f(x), ord = np.inf)
        if ferr < eps:
            print("After %d steps" % (i))
            print("backward_error =", berr)
            print("forward_error =", ferr)
            break
    
    return x

# broyden法非线性方程组
def broyden2(f, x0, x1, B=None, N=25, eps=1e-6):
    if B is None:
        B = np.mat(np.eye(len(x0)))

    for i in range(N):
        delta = np.mat(f(x1) - f(x0)).reshape((-1, 1))
        sigma = np.mat(x1 - x0).reshape((-1, 1))
        B = B + ((sigma - B @ delta) @ sigma.T @ B)/(sigma.T @ B @ delta)
        x0 = x1 - (B @ np.mat(f(x1)).reshape((-1, 1))).ravel()
        x0 = np.array(x0).ravel()

        ferr = np.linalg.norm(x0 - x1, ord = np.inf)
        berr = np.linalg.norm(f(x0), ord = np.inf)
        x1, x0 = x0, x1
        if ferr < eps:
            print("After %d steps" % (i))
            print("backward_error =", berr)
            print("forward_error =", ferr)
            break

    return x1

# LU分解
def lu(A, debug = False):
    A = np.mat(np.array(A).copy(), dtype=np.float_)
    if A.ndim > 2:
        print("Dims of A must be 2")
        return
    row, col = A.shape[0], A.shape[1]
    min_row_col = min(row, col)

    if debug:
        print("A", A, "\nrows = %d cols = %d min_row_col = %d" % (row, col, min_row_col))

    for i in range(row - 1):
        A[i+1:, i] = A[i+1:, i]/A[i, i]
        A[i+1:, i+1:] = A[i+1:, i+1:] - A[i+1:, i] @ A[i, i+1:]

        if debug:
            print(A)

    L = np.tril(A, -1).copy() + np.eye(A.shape[0], A.shape[1])
    U = np.triu(A).copy()
    return L, U

# PA = LU 分解
def plu(A):
    A = np.mat(np.array(A).copy(), dtype=np.float_)
    if A.ndim > 2:
        print("Dims of A must be 2")
        return

    row, col = A.shape[0], A.shape[1]
    min_row_col = min(row, col)
    P = np.eye(A.shape[0], A.shape[1])

    for i in range(row - 1):
        mx_el = np.max(np.abs(A[i:, i]))
        fr, _ = np.where(mx_el == np.abs(A[i:, i]))
        fr = int(fr[0]) + i
        if fr != i:
            A[(i, fr), 0:] = A[(fr, i), 0:]
            P[(i, fr), 0:] = P[(fr, i), 0:]
        A[i+1:, i] = A[i+1:, i]/A[i, i]
        A[i+1:, i+1:] = A[i+1:, i+1:] - A[i+1:, i] * A[i, i+1:]

    U = np.triu(A)
    L = np.tril(A, -1) + np.eye(A.shape[0], A.shape[1])
    return P, L, U