From 59f41ee600fcd01fedd246edabb61b891259b2fe Mon Sep 17 00:00:00 2001 From: "2563491540@qq.com" <2563491540@qq.com> Date: Fri, 24 Feb 2023 09:12:17 +0800 Subject: [PATCH] added two folder containing chapter 6 and 7 --- .../iteration_methods.py | 164 ++++++++++++++++++ .../single_variable_nonlinear_equation.py | 143 +++++++++++++++ 2 files changed, 307 insertions(+) create mode 100644 chapter6_iteration_solve_linear_system_of_equations/iteration_methods.py create mode 100644 chapter7_Numerical_methods_for_solving_nonlinear_equations_and_systems/single_variable_nonlinear_equation.py diff --git a/chapter6_iteration_solve_linear_system_of_equations/iteration_methods.py b/chapter6_iteration_solve_linear_system_of_equations/iteration_methods.py new file mode 100644 index 0000000..35d56ce --- /dev/null +++ b/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)) diff --git a/chapter7_Numerical_methods_for_solving_nonlinear_equations_and_systems/single_variable_nonlinear_equation.py b/chapter7_Numerical_methods_for_solving_nonlinear_equations_and_systems/single_variable_nonlinear_equation.py new file mode 100644 index 0000000..f968356 --- /dev/null +++ b/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)) -- 2.34.1