chapter6_iteration_solve_linear_system_of_equations/iteration_methods.py

@@ -0,0 +1,164 @@
+# !/usr/bin/env python
+# -*- coding: utf-8 -*-
+"""
+@Time        : 2020/12/16 10:33
+@Author      : Albert Darren
+@Contact     : 2563491540@qq.com
+@File        : iteration_methods.py
+@Version     : Version 1.0.0
+@Description : TODO 自己实现雅各比迭代法，高斯-赛德尔迭代法,SOR迭代法，谱半径函数，希尔伯特矩阵函数
+@Created By  : PyCharm
+"""
+
+import numpy as np
+from scipy.linalg import eigvals, inv
+
+import \
+    AI.NumericalAnalysis.chapter5_directive_solve_linear_system_of_equations.DirectiveSolveLinearSystemEquations as dsl
+
+
+def vector_norm(vector: np.ndarray, p=None):
+    """
+    计算向量的p-范数
+    :param vector: 实向量或者复向量
+    :param p: 指定类型的范数，默认是oo范数
+    :return: 指定向量范数和p值
+    """
+    if p is None:
+        return abs(vector).max(), p
+    elif p >= 1:
+        return np.power(np.sum(np.power(abs(vector), p)), 1 / p), p
+    else:
+        raise Exception("error,p must be an integer , greater than  or equal to 1")
+
+
+def jacobi(coefficient_matrix: np.ndarray,
+           right_hand_side_vector: np.ndarray,
+           initial_vector: np.ndarray, epsilon=1e-4):
+    d = np.diag(np.diag(coefficient_matrix))
+    l = -np.tril(coefficient_matrix, k=-1)
+    u = -np.triu(coefficient_matrix, k=1)
+    jacobi_matrix = inv(d) @ (l + u)
+    if spectral_radius(jacobi_matrix) < 1:
+        f = inv(d) @ right_hand_side_vector
+        x_1 = initial_vector
+        x_2 = jacobi_matrix @ x_1 + f
+        iteration_count = 1
+        while 1:
+            norm, _ = vector_norm(x_2 - x_1)
+            if norm < epsilon:
+                return x_2, iteration_count
+            x_1 = x_2
+            x_2 = jacobi_matrix @ x_1 + f
+            iteration_count += 1
+    raise Exception("jacobi iteration is not convergent")
+
+
+def gauss_seidel(coefficient_matrix: np.ndarray,
+                 right_hand_side_vector: np.ndarray,
+                 initial_vector: np.ndarray, epsilon=1e-4):
+    d = np.diag(np.diag(coefficient_matrix))
+    l = -np.tril(coefficient_matrix, k=-1)
+    u = -np.triu(coefficient_matrix, k=1)
+    gauss_seidel_matrix = inv(d - l) @ u
+    if spectral_radius(gauss_seidel_matrix) < 1:
+        f = inv(d - l) @ right_hand_side_vector
+        x_1 = initial_vector
+        x_2 = gauss_seidel_matrix @ x_1 + f
+        iteration_count = 1
+        while 1:
+            norm, _ = vector_norm(x_2 - x_1)
+            if norm < epsilon:
+                return x_2, iteration_count
+            x_1 = x_2
+            x_2 = gauss_seidel_matrix @ x_1 + f
+            iteration_count += 1
+    raise Exception("jacobi iteration is not convergent")
+
+
+def successive_over_relaxation(coefficient_matrix: np.ndarray,
+                               right_hand_side_vector: np.ndarray,
+                               initial_vector: np.ndarray,
+                               true_solution: np.ndarray, omega=1.0, epsilon=5e-6):
+    d = np.diag(np.diag(coefficient_matrix))
+    l = -np.tril(coefficient_matrix, k=-1)
+    u = -np.triu(coefficient_matrix, k=1)
+    l_omega = inv(d - omega * l) @ ((1 - omega) * d + omega * u)
+    if spectral_radius(l_omega) < 1:
+        f = omega * inv(d - omega * l) @ right_hand_side_vector
+        x = l_omega @ initial_vector + f
+        iteration_count = 1
+        while 1:
+            norm, _ = vector_norm(true_solution - x)
+            if norm < epsilon:
+                return x, iteration_count
+            x = l_omega @ x + f
+            iteration_count += 1
+    raise Exception("jacobi iteration is not convergent")
+
+
+def hilbert_matrix(order: int):
+    hilbert = np.zeros((order, order))
+    for i in range(order):
+        for j in range(order):
+            hilbert[i, j] = 1 / (i + j + 1)
+    return hilbert
+
+
+def spectral_radius(square_matrix: np.ndarray):
+    if square_matrix.shape[0] == square_matrix.shape[1]:
+        return abs(eigvals(square_matrix)).max()
+    raise Exception("\n{} is not a square matrix".format(square_matrix))
+
+
+if __name__ == '__main__':
+    # 第一组测试方程组，来源详见李庆扬数值分析第5版P180-181
+    # A = np.array([[8, -3, 2], [4, 11, -1], [6, 3, 12]], dtype=np.float64)
+    # b = np.array([20, 33, 36], dtype=np.float64).reshape((3, 1))
+    # true_solution_vector = dsl.gaussian_elimination(A, b)
+    # b0 = np.array([[0, 3 / 8, -2 / 8], [-4 / 11, 0, 1 / 11], [-6 / 12, -3 / 12, 0]])
+    # print(spectral_radius(b0))
+
+    # 第二组测试方程组，来源详见李庆扬数值分析第5版P209 e.g.1
+    """
+    A = np.array([[5, 2, 1], [-1, 4, 2], [2, -3, 10]], dtype=np.float64)
+    b = np.array([-12, 20, 3], dtype=np.float64).reshape((3, 1))
+    x0 = np.array([1, 1, 1], dtype=np.float64).reshape((3, 1))
+    root1, count1 = jacobi(A, b, x0)
+    print("jacobi公式迭代{}次达到精度要求，\n近似解x=\n{}".format(count1, root1))
+    root2, count2 = gauss_seidel(A, b, x0)
+    print("Gauss_seidel公式迭代{}次达到精度要求，\n近似解x=\n{}".format(count2, root2))
+    """
+    # 第三组测试方程组，来源详见李庆扬数值分析第5版P210 e.g.7
+    """
+    A = np.array([[4, -1, 0], [-1, 4, -1], [0, -1, 4]], dtype=np.float64)
+    b = np.array([1, 4, -3], dtype=np.float64).reshape((3, 1))
+    true_roots = np.array([1 / 2, 1, -1 / 2]).reshape((3, 1))
+    x0 = np.array([0, 0, 0]).reshape((3, 1))
+    # omega=1.03
+    root1, count1 = successive_over_relaxation(A, b, x0, true_roots, 1.03)
+    print("SOR方法,当取Omega=1.03,初值x0=\n{}\n迭代{}次可达到精度要求，此时近似解x=\n{}".format(x0, count1, root1))
+    # omega=1
+    root2, count2 = successive_over_relaxation(A, b, x0, true_roots)
+    print("SOR方法,当取Omega=1,初值x0=\n{}\n迭代{}次可达到精度要求，此时近似解x=\n{}".format(x0, count2, root2))
+    # omega=1.1
+    root3, count3 = successive_over_relaxation(A, b, x0, true_roots, 1.1)
+    print("SOR方法,当取Omega=1.1,初值x0=\n{}\n迭代{}次可达到精度要求，此时近似解x=\n{}".format(x0, count3, root3))
+    """
+    # 第四组测试方程组，来源详见李庆扬数值分析第5版P211 e.g.1
+    A = hilbert_matrix(10)
+    b = A @ np.ones((10, 1))
+    x0 = np.zeros((10, 1))
+    true_roots = dsl.gaussian_elimination(A, b)
+    print("列主元高斯消去法精确解x=\n{}".format(true_roots))
+    # root1, count1 = jacobi(A, b, x0, 1e-6)
+    # print("jacobi公式迭代{}次达到精度要求，\n近似解x=\n{}".format(count1, root1))
+    # omega=1
+    root2, count2 = successive_over_relaxation(A, b, x0, true_roots)
+    print("SOR方法,当取Omega=1,初值x0=\n{}\n迭代{}次可达到精度要求，此时近似解x=\n{}".format(x0, count2, root2))
+    # omega=1.25
+    root3, count3 = successive_over_relaxation(A, b, x0, true_roots, 1.25)
+    print("SOR方法,当取Omega=1.25,初值x0=\n{}\n迭代{}次可达到精度要求，此时近似解x=\n{}".format(x0, count3, root3))
+    # omega=1.5
+    root4, count4 = successive_over_relaxation(A, b, x0, true_roots, 1.5)
+    print("SOR方法,当取Omega=1.5,初值x0=\n{}\n迭代{}次可达到精度要求，此时近似解x=\n{}".format(x0, count4, root4))

chapter7_Numerical_methods_for_solving_nonlinear_equations_and_systems/single_variable_nonlinear_equation.py

@@ -0,0 +1,143 @@
+# !/usr/bin/env python
+# -*- coding: utf-8 -*-
+"""
+@Time        : 2020/12/15 15:58
+@Author      : Albert Darren
+@Contact     : 2563491540@qq.com
+@File        : single_variable_nonlinear_equation.py
+@Version     : Version 1.0.0
+@Description : TODO 自己实现单变量非线性方程组f(x)=0常见数值解法
+@Created By  : PyCharm
+"""
+from sympy.abc import x
+from sympy import init_printing, evalf, exp, solve, Eq
+
+init_printing(use_unicode=True)
+
+
+def bisection(f_x, interval_start: (float, int), interval_end: (float, int), epsilon=5e-3):
+    """
+    二分法求解单变量非线性方程组f(x)=0在区间[a,b]上的准确到小数点后指定精度的近似解
+    :param f_x: 单变量非线性函数f(x)
+    :param interval_start: 求根区间始点
+    :param interval_end: 求根区间终点
+    :param epsilon: 小数点精度
+    :return: 近似解x,二分迭代次数k
+    """
+    if f_x.subs(x, interval_start) * f_x.subs(x, interval_end) < 0:
+        # 初始化二分次数，区间长度
+        bisection_count = 0
+        interval_length = abs(interval_start - interval_end)
+        while 1:
+            middle_point = (interval_start + interval_end) / 2
+            bisection_count += 1
+            value = f_x.subs(x, middle_point)
+            if interval_length <= epsilon or value == 0:  # 达到预定精度epsilon或者等于精确解时结束循环
+                return middle_point, bisection_count
+            if value * f_x.subs(x, interval_start) < 0:
+                interval_end = middle_point
+            else:
+                interval_start = middle_point
+            interval_length = abs(interval_start - interval_end)
+    else:
+        raise Exception("function has no solution in the specified interval")
+
+
+def newtown_tangent(f_x, x0, true_root, epsilon=5e-4):
+    """
+    牛顿法求解单变量非线性方程组f(x)=0在区间[a,b]上的准确到小数点后指定精度的近似解
+    :param f_x: 单变量非线性函数f(x)
+    :param x0: 迭代初始近似值
+    :param true_root: 精确解
+    :param epsilon: 小数点精度
+    :return: 近似解x,二分迭代次数k
+    """
+    f = f_x.subs(x, x0)
+    diff_f_x = f_x.diff()
+    diff_f = diff_f_x.subs(x, x0)
+    iteration_count = 0
+    while 1:
+        x_k = x0 - f / diff_f
+        iteration_count += 1
+        if abs(x_k - true_root) < epsilon:
+            return x_k, iteration_count
+        x0 = x_k
+        f = f_x.subs(x, x0)
+        diff_f = diff_f_x.subs(x, x0)
+
+
+def reduce_newtown():
+    pass
+
+
+def newtown_downhill():
+    pass
+
+
+def secant(f_x, x0, x1, true_root, epsilon=5e-4):
+    """
+    弦截法求解单变量非线性方程组f(x)=0在区间[a,b]上的准确到小数点后指定精度的近似解
+    :param f_x: 单变量非线性函数f(x)
+    :param x0: 第一个迭代初始近似值
+    :param x1: 第二个迭代初始近似值
+    :param true_root: 精确解
+    :param epsilon: 小数点精度
+    :return: 近似解x,二分迭代次数k
+    """
+    f_x0 = f_x.subs(x, x0)
+    f_x1 = f_x.subs(x, x1)
+    difference_quotient = (f_x1 - f_x0) / (x1 - x0)
+    iteration_count = 0
+    while 1:
+        x_k = x1 - f_x1 / difference_quotient
+        iteration_count += 1
+        if abs(x_k - true_root) < epsilon:
+            return x_k, iteration_count
+        x0 = x1
+        x1 = x_k
+        f_x0 = f_x.subs(x, x0)
+        f_x1 = f_x.subs(x, x1)
+        difference_quotient = (f_x1 - f_x0) / (x1 - x0)
+
+
+def parabola():
+    pass
+
+
+if __name__ == '__main__':
+    # 第一组测试函数，来源详见李庆扬数值分析第5版P214，e.g.2
+
+    # f = x ** 3 - x - 1
+    # appropriate_root1, iteration_count1 = bisection(f, 1.0, 1.5)
+    # print("二分法求方程近似解x={},迭代次数={}".format(appropriate_root1,iteration_count1))
+
+    # 测试区间中点恰好是精确解时的情况
+    # f = x - 1
+    # appropriate_root2, iteration_count2 = bisection(f, 0.0, 2.0)
+    # print("测试区间中点恰好是精确解时,近似解x={},迭代次数={}".format(appropriate_root2, iteration_count2))
+    # 第二组测试函数，来源详见李庆扬数值分析第5版P238，e.g.7
+
+    # f = x ** 3 - 3 * x - 1
+    # initial_value = 2
+    # true_solution = 1.87938524
+    # root1, count1 = newtown_tangent(f, initial_value, true_solution)
+    # print("牛顿切线法近似解x={},迭代次数={}".format(root1.evalf(), count1))
+    # initial_value0, initial_value1 = 2, 1.9
+    # root2, count2 = secant(f, initial_value0, initial_value1, true_solution)
+    # print("弦截法近似解x={},迭代次数={}".format(root2.evalf(), count2))
+    # 第三组测试函数，来源详见李庆扬数值分析第5版P224-225
+
+    # c = 115
+    # f = x ** 2 - c
+    # true_solution = pow(c, 0.5)
+    # initial_value = 11
+    # solution, count = newtown_tangent(f, initial_value, true_solution, epsilon=1e-6)
+    # print("牛顿切线法近似解x={},迭代次数={}".format(solution.evalf(), count))
+    # 第四组测试函数，来源详见李庆扬数值分析第5版P229 e.g.10
+
+    f = x * exp(x) - 1
+    root = solve(Eq(f, 0), x)[0]
+    print("调用sympy函数库解x={}".format(root))
+    initial_value0, initial_value1 = 0.5, 0.6
+    root, count = secant(f, initial_value0, initial_value1, root)
+    print("牛顿切线法近似解x={},迭代次数={}".format(root.evalf(), count))