用Python复现数值分析经典算法从插值到方程组求解附Jupyter Notebook数值分析作为连接数学理论与工程实践的桥梁其核心算法在科学计算、金融建模、机器学习等领域有着广泛应用。对于编程学习者而言仅理解数学推导远远不够——将课本公式转化为可执行的代码才是真正掌握这些算法的关键。本文将带您用Python的NumPy和SciPy库在Jupyter Notebook中重现Lagrange插值、最小二乘拟合、Gauss消元法等经典算法通过可视化计算过程深化理解并验证教材中的经典例题结果。1. 环境准备与基础工具在开始算法实现前需要配置适合科学计算的Python环境。推荐使用Anaconda发行版它集成了Jupyter Notebook和主要科学计算库。以下是基础环境的搭建步骤conda create -n numerical_analysis python3.9 conda activate numerical_analysis conda install numpy scipy matplotlib pandas sympy核心库的功能定位NumPy提供高效的数组运算和线性代数操作SciPy包含优化、积分、插值等高级数学函数Matplotlib实现数据可视化直观展示计算结果SymPy用于符号计算可验证解析解提示在Jupyter Notebook中使用%matplotlib inline命令可使图表直接显示在单元格下方2. 插值法从理论到代码实现插值法是数值分析的基础工具用于通过已知数据点构建近似函数。我们以经典的Lagrange插值为起点逐步实现更复杂的算法。2.1 Lagrange插值实现Lagrange插值多项式定义为$$ L(x) \sum_{i0}^{n} y_i \prod_{\substack{j0 \ j \neq i}}^{n} \frac{x-x_j}{x_i-x_j} $$Python实现代码import numpy as np import matplotlib.pyplot as plt def lagrange_interp(x_points, y_points, x): n len(x_points) result 0.0 for i in range(n): term y_points[i] for j in range(n): if i ! j: term * (x - x_points[j]) / (x_points[i] - x_points[j]) result term return result # 课本例题验证 x_data np.array([0, 1, 2, 3]) y_data np.array([1, 2, 9, 28]) x_test np.linspace(-0.5, 3.5, 100) y_interp [lagrange_interp(x_data, y_data, x) for x in x_test] plt.figure(figsize(10,6)) plt.scatter(x_data, y_data, colorred, labelData points) plt.plot(x_test, y_interp, labelLagrange interpolation) plt.legend(); plt.grid(True) plt.title(Lagrange Interpolation Example)2.2 分段线性与三次样条插值当插值点较多时高次多项式会出现Runge现象。此时分段插值更为稳定from scipy.interpolate import interp1d # 分段线性插值 linear_interp interp1d(x_data, y_data, kindlinear) # 三次样条插值 cubic_interp interp1d(x_data, y_data, kindcubic) # 可视化对比 plt.figure(figsize(12,6)) plt.scatter(x_data, y_data, colorred, labelData points) plt.plot(x_test, linear_interp(x_test), labelPiecewise linear) plt.plot(x_test, cubic_interp(x_test), labelCubic spline) plt.legend(); plt.grid(True)3. 数据拟合最小二乘法实践当数据存在噪声时拟合比插值更为合适。我们以线性回归为例演示最小二乘法的实现。3.1 线性最小二乘拟合对于模型$y a bx$正规方程的解为$$ \begin{bmatrix} a \ b \end{bmatrix} (X^TX)^{-1}X^Ty \quad \text{其中} \quad X \begin{bmatrix} 1 x_1 \ 1 x_2 \ \vdots \vdots \ 1 x_n \end{bmatrix} $$Python实现def least_squares_fit(x, y): X np.vstack([np.ones_like(x), x]).T coeff np.linalg.inv(X.T X) X.T y return coeff # 生成带噪声的数据 np.random.seed(42) x_data np.linspace(0, 10, 20) y_true 2 1.5 * x_data y_noisy y_true np.random.normal(0, 2, sizelen(x_data)) # 拟合并可视化 a, b least_squares_fit(x_data, y_noisy) plt.figure(figsize(10,6)) plt.scatter(x_data, y_noisy, labelNoisy data) plt.plot(x_data, a b*x_data, r, labelfFit: y{a:.2f}{b:.2f}x) plt.plot(x_data, y_true, g--, labelTrue relationship) plt.legend(); plt.grid(True)3.2 非线性拟合示例对于非线性模型可使用scipy.optimize.curve_fitfrom scipy.optimize import curve_fit def exp_model(x, a, b, c): return a * np.exp(b * x) c popt, pcov curve_fit(exp_model, x_data, y_noisy) plt.figure(figsize(10,6)) plt.scatter(x_data, y_noisy) plt.plot(x_data, exp_model(x_data, *popt), r-)4. 线性方程组求解方法解线性方程组是数值分析的核心问题我们比较直接法和迭代法的实现。4.1 Gauss消元法实现基本Gauss消元法包含三个步骤前向消元主元归一化回代求解def gauss_elimination(A, b): n len(b) # 前向消元 for k in range(n-1): for i in range(k1, n): factor A[i,k] / A[k,k] A[i,k:] - factor * A[k,k:] b[i] - factor * b[k] # 回代求解 x np.zeros(n) for i in range(n-1, -1, -1): x[i] (b[i] - np.dot(A[i,i1:], x[i1:])) / A[i,i] return x # 课本例题验证 A np.array([[2, -1, 1], [3, 3, 9], [3, 3, 5]], dtypefloat) b np.array([2, 15, 10], dtypefloat) solution gauss_elimination(A.copy(), b.copy()) print(f解向量: {solution})4.2 迭代法Jacobi与Gauss-Seidel对于大型稀疏矩阵迭代法更为高效def jacobi(A, b, max_iter100, tol1e-8): n len(b) x np.zeros(n) for _ in range(max_iter): x_new np.zeros(n) for i in range(n): s np.dot(A[i,:i], x[:i]) np.dot(A[i,i1:], x[i1:]) x_new[i] (b[i] - s) / A[i,i] if np.linalg.norm(x_new - x) tol: break x x_new return x def gauss_seidel(A, b, max_iter100, tol1e-8): n len(b) x np.zeros(n) for _ in range(max_iter): x_old x.copy() for i in range(n): s np.dot(A[i,:i], x[:i]) np.dot(A[i,i1:], x[i1:]) x[i] (b[i] - s) / A[i,i] if np.linalg.norm(x - x_old) tol: break return x # 收敛矩阵测试 A_diag np.array([[4, -1, 0], [-1, 4, -1], [0, -1, 4]]) b_test np.array([1, 1, 1]) print(fJacobi解: {jacobi(A_diag, b_test)}) print(fGauss-Seidel解: {gauss_seidel(A_diag, b_test)})5. 数值积分方法比较数值积分在无法求得解析积分时尤为重要我们实现梯形法和Simpson法。5.1 复合梯形公式$$ \int_a^b f(x)dx \approx \frac{h}{2}\left[f(a) 2\sum_{i1}^{n-1}f(x_i) f(b)\right] $$def trapezoidal_rule(f, a, b, n100): h (b - a) / n x np.linspace(a, b, n1) y f(x) return h * (0.5*y[0] 0.5*y[-1] np.sum(y[1:-1])) # 测试积分 ∫_0^π sin(x)dx result trapezoidal_rule(np.sin, 0, np.pi, 100) print(f梯形法结果: {result} (误差: {abs(2-result):.2e}))5.2 Simpson法实现复合Simpson公式具有更高精度def simpson_rule(f, a, b, n100): h (b - a) / n x np.linspace(a, b, n1) y f(x) return h/3 * (y[0] y[-1] 4*np.sum(y[1:-1:2]) 2*np.sum(y[2:-2:2])) result simpson_rule(np.sin, 0, np.pi, 10) print(fSimpson法结果: {result} (误差: {abs(2-result):.2e}))6. 非线性方程求解非线性方程的数值解法包括二分法、Newton法等各有适用场景。6.1 二分法实现二分法简单但收敛速度较慢def bisection(f, a, b, tol1e-6, max_iter100): if f(a)*f(b) 0: raise ValueError(函数在区间端点同号) for _ in range(max_iter): c (a b) / 2 if abs(f(c)) tol or (b-a)/2 tol: return c if f(c)*f(a) 0: b c else: a c return (a b) / 2 # 求解x^3 - x - 2 0 f lambda x: x**3 - x - 2 root bisection(f, 1, 2) print(f二分法求根结果: {root:.6f})6.2 Newton迭代法Newton法具有二次收敛速度但需要导数信息def newton_method(f, df, x0, tol1e-6, max_iter100): x x0 for _ in range(max_iter): fx f(x) if abs(fx) tol: return x dfx df(x) if dfx 0: break x - fx / dfx return x f lambda x: x**3 - x - 2 df lambda x: 3*x**2 - 1 root newton_method(f, df, 1.5) print(fNewton法求根结果: {root:.6f})在实际项目中选择算法需要权衡收敛速度、计算复杂度和稳定性。例如对于病态矩阵方程组直接法可能产生较大误差此时预处理后的迭代法更为可靠。
用Python复现数值分析经典算法：从插值到方程组求解（附Jupyter Notebook）
用Python复现数值分析经典算法从插值到方程组求解附Jupyter Notebook数值分析作为连接数学理论与工程实践的桥梁其核心算法在科学计算、金融建模、机器学习等领域有着广泛应用。对于编程学习者而言仅理解数学推导远远不够——将课本公式转化为可执行的代码才是真正掌握这些算法的关键。本文将带您用Python的NumPy和SciPy库在Jupyter Notebook中重现Lagrange插值、最小二乘拟合、Gauss消元法等经典算法通过可视化计算过程深化理解并验证教材中的经典例题结果。1. 环境准备与基础工具在开始算法实现前需要配置适合科学计算的Python环境。推荐使用Anaconda发行版它集成了Jupyter Notebook和主要科学计算库。以下是基础环境的搭建步骤conda create -n numerical_analysis python3.9 conda activate numerical_analysis conda install numpy scipy matplotlib pandas sympy核心库的功能定位NumPy提供高效的数组运算和线性代数操作SciPy包含优化、积分、插值等高级数学函数Matplotlib实现数据可视化直观展示计算结果SymPy用于符号计算可验证解析解提示在Jupyter Notebook中使用%matplotlib inline命令可使图表直接显示在单元格下方2. 插值法从理论到代码实现插值法是数值分析的基础工具用于通过已知数据点构建近似函数。我们以经典的Lagrange插值为起点逐步实现更复杂的算法。2.1 Lagrange插值实现Lagrange插值多项式定义为$$ L(x) \sum_{i0}^{n} y_i \prod_{\substack{j0 \ j \neq i}}^{n} \frac{x-x_j}{x_i-x_j} $$Python实现代码import numpy as np import matplotlib.pyplot as plt def lagrange_interp(x_points, y_points, x): n len(x_points) result 0.0 for i in range(n): term y_points[i] for j in range(n): if i ! j: term * (x - x_points[j]) / (x_points[i] - x_points[j]) result term return result # 课本例题验证 x_data np.array([0, 1, 2, 3]) y_data np.array([1, 2, 9, 28]) x_test np.linspace(-0.5, 3.5, 100) y_interp [lagrange_interp(x_data, y_data, x) for x in x_test] plt.figure(figsize(10,6)) plt.scatter(x_data, y_data, colorred, labelData points) plt.plot(x_test, y_interp, labelLagrange interpolation) plt.legend(); plt.grid(True) plt.title(Lagrange Interpolation Example)2.2 分段线性与三次样条插值当插值点较多时高次多项式会出现Runge现象。此时分段插值更为稳定from scipy.interpolate import interp1d # 分段线性插值 linear_interp interp1d(x_data, y_data, kindlinear) # 三次样条插值 cubic_interp interp1d(x_data, y_data, kindcubic) # 可视化对比 plt.figure(figsize(12,6)) plt.scatter(x_data, y_data, colorred, labelData points) plt.plot(x_test, linear_interp(x_test), labelPiecewise linear) plt.plot(x_test, cubic_interp(x_test), labelCubic spline) plt.legend(); plt.grid(True)3. 数据拟合最小二乘法实践当数据存在噪声时拟合比插值更为合适。我们以线性回归为例演示最小二乘法的实现。3.1 线性最小二乘拟合对于模型$y a bx$正规方程的解为$$ \begin{bmatrix} a \ b \end{bmatrix} (X^TX)^{-1}X^Ty \quad \text{其中} \quad X \begin{bmatrix} 1 x_1 \ 1 x_2 \ \vdots \vdots \ 1 x_n \end{bmatrix} $$Python实现def least_squares_fit(x, y): X np.vstack([np.ones_like(x), x]).T coeff np.linalg.inv(X.T X) X.T y return coeff # 生成带噪声的数据 np.random.seed(42) x_data np.linspace(0, 10, 20) y_true 2 1.5 * x_data y_noisy y_true np.random.normal(0, 2, sizelen(x_data)) # 拟合并可视化 a, b least_squares_fit(x_data, y_noisy) plt.figure(figsize(10,6)) plt.scatter(x_data, y_noisy, labelNoisy data) plt.plot(x_data, a b*x_data, r, labelfFit: y{a:.2f}{b:.2f}x) plt.plot(x_data, y_true, g--, labelTrue relationship) plt.legend(); plt.grid(True)3.2 非线性拟合示例对于非线性模型可使用scipy.optimize.curve_fitfrom scipy.optimize import curve_fit def exp_model(x, a, b, c): return a * np.exp(b * x) c popt, pcov curve_fit(exp_model, x_data, y_noisy) plt.figure(figsize(10,6)) plt.scatter(x_data, y_noisy) plt.plot(x_data, exp_model(x_data, *popt), r-)4. 线性方程组求解方法解线性方程组是数值分析的核心问题我们比较直接法和迭代法的实现。4.1 Gauss消元法实现基本Gauss消元法包含三个步骤前向消元主元归一化回代求解def gauss_elimination(A, b): n len(b) # 前向消元 for k in range(n-1): for i in range(k1, n): factor A[i,k] / A[k,k] A[i,k:] - factor * A[k,k:] b[i] - factor * b[k] # 回代求解 x np.zeros(n) for i in range(n-1, -1, -1): x[i] (b[i] - np.dot(A[i,i1:], x[i1:])) / A[i,i] return x # 课本例题验证 A np.array([[2, -1, 1], [3, 3, 9], [3, 3, 5]], dtypefloat) b np.array([2, 15, 10], dtypefloat) solution gauss_elimination(A.copy(), b.copy()) print(f解向量: {solution})4.2 迭代法Jacobi与Gauss-Seidel对于大型稀疏矩阵迭代法更为高效def jacobi(A, b, max_iter100, tol1e-8): n len(b) x np.zeros(n) for _ in range(max_iter): x_new np.zeros(n) for i in range(n): s np.dot(A[i,:i], x[:i]) np.dot(A[i,i1:], x[i1:]) x_new[i] (b[i] - s) / A[i,i] if np.linalg.norm(x_new - x) tol: break x x_new return x def gauss_seidel(A, b, max_iter100, tol1e-8): n len(b) x np.zeros(n) for _ in range(max_iter): x_old x.copy() for i in range(n): s np.dot(A[i,:i], x[:i]) np.dot(A[i,i1:], x[i1:]) x[i] (b[i] - s) / A[i,i] if np.linalg.norm(x - x_old) tol: break return x # 收敛矩阵测试 A_diag np.array([[4, -1, 0], [-1, 4, -1], [0, -1, 4]]) b_test np.array([1, 1, 1]) print(fJacobi解: {jacobi(A_diag, b_test)}) print(fGauss-Seidel解: {gauss_seidel(A_diag, b_test)})5. 数值积分方法比较数值积分在无法求得解析积分时尤为重要我们实现梯形法和Simpson法。5.1 复合梯形公式$$ \int_a^b f(x)dx \approx \frac{h}{2}\left[f(a) 2\sum_{i1}^{n-1}f(x_i) f(b)\right] $$def trapezoidal_rule(f, a, b, n100): h (b - a) / n x np.linspace(a, b, n1) y f(x) return h * (0.5*y[0] 0.5*y[-1] np.sum(y[1:-1])) # 测试积分 ∫_0^π sin(x)dx result trapezoidal_rule(np.sin, 0, np.pi, 100) print(f梯形法结果: {result} (误差: {abs(2-result):.2e}))5.2 Simpson法实现复合Simpson公式具有更高精度def simpson_rule(f, a, b, n100): h (b - a) / n x np.linspace(a, b, n1) y f(x) return h/3 * (y[0] y[-1] 4*np.sum(y[1:-1:2]) 2*np.sum(y[2:-2:2])) result simpson_rule(np.sin, 0, np.pi, 10) print(fSimpson法结果: {result} (误差: {abs(2-result):.2e}))6. 非线性方程求解非线性方程的数值解法包括二分法、Newton法等各有适用场景。6.1 二分法实现二分法简单但收敛速度较慢def bisection(f, a, b, tol1e-6, max_iter100): if f(a)*f(b) 0: raise ValueError(函数在区间端点同号) for _ in range(max_iter): c (a b) / 2 if abs(f(c)) tol or (b-a)/2 tol: return c if f(c)*f(a) 0: b c else: a c return (a b) / 2 # 求解x^3 - x - 2 0 f lambda x: x**3 - x - 2 root bisection(f, 1, 2) print(f二分法求根结果: {root:.6f})6.2 Newton迭代法Newton法具有二次收敛速度但需要导数信息def newton_method(f, df, x0, tol1e-6, max_iter100): x x0 for _ in range(max_iter): fx f(x) if abs(fx) tol: return x dfx df(x) if dfx 0: break x - fx / dfx return x f lambda x: x**3 - x - 2 df lambda x: 3*x**2 - 1 root newton_method(f, df, 1.5) print(fNewton法求根结果: {root:.6f})在实际项目中选择算法需要权衡收敛速度、计算复杂度和稳定性。例如对于病态矩阵方程组直接法可能产生较大误差此时预处理后的迭代法更为可靠。
