1. 从“猜数字”到“找交点”方程求根的直观理解在数学建模和科学计算的无数场景里我们常常会遇到一个核心问题如何找到一个方程的解这个问题听起来简单比如求解x^2 - 2 0我们立刻能想到答案是x √2或x -√2。但在实际工程和科研中方程往往复杂得多比如描述化学反应速率的非线性方程、描述经济模型的隐式方程或者像x - cos(x) 0这样看似简单却无法用初等函数写出解析解的方程。这时候我们就需要依赖数值方法通过计算机“算”出满足一定精度的近似解。这个过程就是“方程求根”。你可以把它想象成一个“猜数字”游戏我们不知道方程f(x) 0的根x*具体是多少但我们可以不断尝试不同的x值计算对应的f(x)并根据f(x)的符号和大小智能地调整下一次猜测的方向和步长最终逼近那个使函数值为零的点。从几何上看方程f(x) 0的根就是函数y f(x)图像与x轴的交点。数值求根法本质上就是设计一套高效的策略让计算机自动、快速地找到这个或这些交点。今天我们就深入探讨两种最基础、也最强大的数值求根方法一般迭代法和牛顿法。我会结合 MATLAB 和 Python 这两种在数学建模和工程领域最常用的工具不仅告诉你算法公式怎么写更会剖析它们背后的数学思想、收敛的奥秘、编程实现的细节以及在实际建模中如何选择、如何避坑。无论你是正在备战数学建模竞赛的学生还是需要处理实际工程问题的工程师理解并掌握这两种方法都将是你的核心技能之一。2. 不动点迭代将求根问题转化为“寻找收敛点”一般迭代法更学术的名称是“不动点迭代法”。它的核心思想非常巧妙不是直接去解f(x) 0而是将原方程等价地改写为x g(x)的形式。这个新函数g(x)称为迭代函数。方程x g(x)的解x*满足x* g(x*)我们称之为函数g(x)的一个“不动点”。原方程f(x)0的根就是这个不动点。为什么这么做因为x g(x)的形式天然地暗示了一种迭代策略从一个初始猜测值x0开始我们通过公式x_{k1} g(x_k)不断地产生一个新的序列x1, x2, x3, ...。如果这个序列收敛到某个极限x*并且函数g(x)连续那么两边取极限就有x* g(x*)即x*就是不动点也就是原方程的根。2.1 构造迭代函数艺术与科学的结合将f(x)0改写为x g(x)的方式有无数种但并非每一种都能保证迭代收敛。这是迭代法的第一个关键点也是体现经验的地方。举个例子求解方程f(x) x^3 - x - 1 0。构造方式Ax g1(x) x^3 - 1。直接从f(x)0移项得到x^3 x 1再开立方不这里我们直接写成x (x1)^(1/3)是更好的选择但先看这个g1。从x01.5开始迭代x1 1.5^3 - 1 3.375 - 1 2.375x2 2.375^3 - 1 ≈ 13.396 - 1 12.396数值迅速发散这不是我们想要的。构造方式Bx g2(x) (x 1)^(1/3)。这是从x^3 x 1两边开立方根得到的。从同样的x01.5开始x1 (1.51)^(1/3) ≈ 1.357x2 (1.3571)^(1/3) ≈ 1.331x3 ≈ 1.326x4 ≈ 1.324... 序列稳定地向真实根约1.3247靠近。为什么g1发散而g2收敛秘密在于迭代函数g(x)在根x*附近的导数绝对值|g(x*)|。这引出了不动点迭代收敛的局部收敛定理设x*是x g(x)的一个不动点如果迭代函数g(x)在包含x*的某个开区间内连续可微且满足|g(x*)| 1那么存在一个以x*为中心的区间从该区间内任意一点x0出发的迭代序列{x_k}都收敛到x*。并且当0 |g(x*)| 1时收敛是线性的误差大致按等比数列减少当g(x*) 0时收敛速度可能更快如平方收敛。对于g1(x)x^3-1g1(x)3x^2在根x*≈1.3247处|g1(x*)| ≈ 3*(1.3247)^2 ≈ 5.26 1不满足收敛条件。 对于g2(x)(x1)^(1/3)g2(x) (1/3)*(x1)^(-2/3)在x*≈1.3247处|g2(x*)| ≈ 0.2 1满足收敛条件。实操心得在构造迭代函数时我们的目标是让g(x)在根附近的导数绝对值尽可能小最好小于1。常见技巧包括将x单独解出放在等式一边利用代数变形如分子有理化有时甚至需要引入一个松弛因子。没有绝对最好的公式需要结合函数图像和初步判断进行尝试。2.2 MATLAB与Python实现代码与收敛控制理解了原理实现起来就清晰了。我们需要一个循环反复计算x_new g(x_old)并设置合理的停止条件。停止条件通常有两个误差条件|x_new - x_old| tol其中tol是预设的误差容限如1e-6。残差条件|f(x_new)| tol直接看函数值是否足够接近零。 通常两者结合使用并设置最大迭代次数以防不收敛陷入死循环。MATLAB实现示例function [root, iter] fixed_point_iteration(g, x0, tol, max_iter) % FIXED_POINT_ITERATION 不动点迭代法求根 % 输入: g - 迭代函数句柄, x0 - 初始猜测值 % tol - 容差, max_iter - 最大迭代次数 % 输出: root - 求得的根, iter - 实际迭代次数 x_old x0; iter 0; fprintf(迭代过程:\n); fprintf(k\t\t x_k\t\t\t |x_k - x_{k-1}|\n); fprintf(-----------------------------------\n); fprintf(0\t\t %.10f\n, x_old); for iter 1:max_iter x_new g(x_old); % 核心迭代步骤 diff abs(x_new - x_old); fprintf(%d\t\t %.10f\t\t %.10e\n, iter, x_new, diff); % 检查收敛条件 if diff tol root x_new; fprintf(在容差 %.1e 下收敛迭代次数: %d\n, tol, iter); return; end x_old x_new; % 为下一次迭代更新值 end warning(达到最大迭代次数 %d 仍未收敛当前值: %.10f, max_iter, x_new); root x_new; end % 使用示例求解 x - cos(x) 0 构造 g(x) cos(x) g (x) cos(x); x0 0.5; % 初始猜测从图像看根在0.7附近0.5也可以 tol 1e-8; max_iter 100; [root, iter] fixed_point_iteration(g, x0, tol, max_iter); fprintf(求得近似根: %.10f\n, root);Python实现示例import numpy as np def fixed_point_iteration(g, x0, tol1e-8, max_iter100): 不动点迭代法求根 参数: g: 迭代函数 callable x0: 初始猜测值 float tol: 容差 float max_iter: 最大迭代次数 int 返回: root: 求得的根 float iter_count: 实际迭代次数 int history: 迭代历史记录 list of tuples (k, x_k, diff) x_old x0 history [(0, x_old, None)] # 记录迭代历史 (迭代次数, x值, 变化量) print(迭代过程:) print(f{k:4} {x_k:20} {|x_k - x_{k-1}|:20}) print(- * 50) print(f{0:4} {x_old:20.10f}) for k in range(1, max_iter 1): x_new g(x_old) # 核心迭代步骤 diff abs(x_new - x_old) history.append((k, x_new, diff)) print(f{k:4} {x_new:20.10f} {diff:20.10e}) # 检查收敛条件 if diff tol: print(f在容差 {tol:.1e} 下收敛迭代次数: {k}) return x_new, k, history x_old x_new # 为下一次迭代更新值 print(f警告: 达到最大迭代次数 {max_iter} 仍未收敛最后值: {x_new:.10f}) return x_new, max_iter, history # 使用示例求解 x - cos(x) 0 构造 g(x) cos(x) import math g_func lambda x: math.cos(x) x0 0.5 root, iters, hist fixed_point_iteration(g_func, x0) print(f求得近似根: {root:.10f}) # 可选绘制收敛过程 import matplotlib.pyplot as plt iters_list, vals_list, _ zip(*hist) plt.figure(figsize(10, 5)) plt.subplot(1, 2, 1) plt.plot(iters_list, vals_list, bo-, linewidth2, markersize4) plt.xlabel(迭代次数 k) plt.ylabel(迭代值 x_k) plt.title(迭代值随迭代次数的变化) plt.grid(True, alpha0.3) plt.subplot(1, 2, 2) # 绘制函数和yx直线直观展示不动点 x_range np.linspace(0, 1, 400) y_g [g_func(xi) for xi in x_range] plt.plot(x_range, x_range, k--, labely x) plt.plot(x_range, y_g, r-, labely g(x) cos(x)) # 绘制迭代轨迹蛛网图 for i in range(min(10, len(hist)-1)): x_k, x_k_next hist[i][1], hist[i1][1] plt.plot([x_k, x_k], [x_k, x_k_next], b:) plt.plot([x_k, x_k_next], [x_k_next, x_k_next], b:) plt.scatter([root], [root], colorgreen, s100, zorder5, labelf不动点 (~{root:.4f})) plt.xlabel(x) plt.ylabel(y) plt.title(不动点迭代的几何图示蛛网图) plt.legend() plt.grid(True, alpha0.3) plt.tight_layout() plt.show()代码解读与避坑指南初始值选择迭代法收敛是“局部”的意味着初始猜测x0必须足够靠近真实的根x*。如果离得太远即使|g(x*)|1序列也可能发散。通常需要结合函数图像或通过简单试算如代入几个点看f(x)符号变化来确定根的粗略区间。停止准则示例中使用了相邻迭代值的绝对差。在实际高精度计算中有时也使用相对误差|x_new - x_old| / |x_new|特别是当根的数量级未知或可能接近0时。同时检查|f(x_new)|是个好习惯确保我们找到的点确实使函数值接近零。收敛诊断输出每一步的迭代值x_k和变化量diff非常有用。你可以观察diff是否在单调递减线性收敛的典型特征如果出现震荡或发散应立即停止并检查迭代函数g(x)的构造或初始值。发散处理代码中设置了最大迭代次数max_iter和警告信息这是防止程序无限循环的必要措施。如果迭代不收敛你需要重新审视问题是不是g(x*)的绝对值大于等于1是不是初始值离根太远是不是应该换一种迭代函数构造方法3. 牛顿法利用导数信息的“超级加速”如果说一般迭代法是“稳步推进”那么牛顿法就是“精准制导”。它利用了函数更多的信息——不仅知道函数值f(x)还知道其导数f(x)——从而实现了更快的收敛速度。牛顿法的迭代公式来源于函数在当前点x_k处的线性近似泰勒展开一阶项。我们想找f(x)0的根。在当前点x_k处将f(x)线性近似f(x) ≈ f(x_k) f(x_k)(x - x_k)。令这个线性近似等于零解出x就得到了下一个迭代点x_{k1}f(x_k) f(x_k)(x_{k1} - x_k) 0x_{k1} x_k - f(x_k) / f(x_k)这个公式有着极其直观的几何解释在曲线上点(x_k, f(x_k))处作切线这条切线与x轴的交点的横坐标就是x_{k1}。因此牛顿法也常被称为“切线法”。3.1 牛顿法的收敛性与优势牛顿法最吸引人的地方在于其局部平方收敛性。在根x*附近如果f(x*) ≠ 0且f(x)连续那么牛顿法产生的误差e_k |x_k - x*|满足e_{k1} ≈ C * (e_k)^2。这意味着每迭代一次有效位数大约会翻倍相比之下线性收敛的一般迭代法e_{k1} ≈ |g(x*)| * e_k要慢得多。但是牛顿法的“美好”是有严格前提的初始值必须足够好平方收敛是“局部”性质如果初始猜测x0离根太远牛顿法完全可能发散。导数不能为零公式中需要除以f(x_k)。如果在迭代过程中某点的导数为零或非常接近零计算将溢出或产生巨大步长导致失败。根x*处的导数f(x*)也不能为零否则是重根收敛速度会降级。需要计算导数你必须能提供导数函数f(x)的表达式或计算方法。对于复杂函数手动求导可能困难。3.2 动手实现牛顿法细节决定成败理解了风险和优势我们来看如何稳健地实现牛顿法。核心迭代公式x_new x_old - f(x_old)/f_prime(x_old)很简单但实现细节关乎成败。MATLAB实现示例function [root, iter, history] newton_method(f, f_prime, x0, tol, max_iter) % NEWTON_METHOD 牛顿法求根 % 输入: f - 原函数句柄, f_prime - 导函数句柄 % x0 - 初始猜测值, tol - 容差, max_iter - 最大迭代次数 % 输出: root - 求得的根, iter - 实际迭代次数, history - 迭代历史 x_old x0; iter 0; history [iter, x_old, nan, nan]; % [k, x_k, f(x_k), |diff|] fprintf(牛顿法迭代过程:\n); fprintf(k\t\t x_k\t\t\t f(x_k)\t\t\t |x_k - x_{k-1}|\n); fprintf(----------------------------------------------------------------\n); fprintf(%d\t\t %.10f\t %.10e\t\t N/A\n, iter, x_old, f(x_old)); for iter 1:max_iter f_val f(x_old); f_prime_val f_prime(x_old); % **关键保护检查导数是否为零或极小** if abs(f_prime_val) eps * 10 % eps是MATLAB的浮点数精度 warning(在 x %.10f 处导数值过小 (|f| %.2e) 迭代终止。, x_old, f_prime_val); root x_old; return; end % 核心牛顿迭代步骤 x_new x_old - f_val / f_prime_val; diff abs(x_new - x_old); history [history; [iter, x_new, f(x_new), diff]]; fprintf(%d\t\t %.10f\t %.10e\t %.10e\n, iter, x_new, f(x_new), diff); % 收敛判断通常结合差值和残差 if diff tol abs(f(x_new)) tol root x_new; fprintf(在容差 %.1e 下收敛迭代次数: %d\n, tol, iter); return; end x_old x_new; end warning(达到最大迭代次数 %d 仍未收敛最后值: %.10f, max_iter, x_new); root x_new; end % 使用示例求解 f(x) x^2 - 2 0 (求根号2) f (x) x^2 - 2; f_prime (x) 2*x; x0 1.5; % 初始猜测真实根约为1.4142 tol 1e-12; max_iter 20; [root, iter, hist] newton_method(f, f_prime, x0, tol, max_iter); fprintf(sqrt(2)的近似值: %.15f\n, root); fprintf(与真实值的绝对误差: %.2e\n, abs(root - sqrt(2)));Python实现示例def newton_method(f, f_prime, x0, tol1e-12, max_iter20): 牛顿法求根 参数: f: 原函数 callable f_prime: 导函数 callable x0: 初始猜测值 float tol: 容差 float max_iter: 最大迭代次数 int 返回: root: 求得的根 float iter_count: 实际迭代次数 int history: 迭代历史 list of [k, x_k, f(x_k), |diff|] x_old x0 history [[0, x_old, f(x_old), None]] print(牛顿法迭代过程:) print(f{k:4} {x_k:20} {f(x_k):25} {|x_k - x_{k-1}|:20}) print(- * 80) print(f{0:4} {x_old:20.15f} {f(x_old):25.15e} {N/A:20}) for k in range(1, max_iter 1): f_val f(x_old) f_prime_val f_prime(x_old) # **关键保护检查导数是否为零或极小** if abs(f_prime_val) 1e-15: # 一个很小的阈值 print(f警告: 在 x {x_old:.10f} 处导数值过小 (|f| {f_prime_val:.2e}) 迭代终止。) return x_old, k-1, history # 核心牛顿迭代步骤 x_new x_old - f_val / f_prime_val diff abs(x_new - x_old) f_new f(x_new) history.append([k, x_new, f_new, diff]) print(f{k:4} {x_new:20.15f} {f_new:25.15e} {diff:20.15e}) # 收敛判断 if diff tol and abs(f_new) tol: print(f在容差 {tol:.1e} 下收敛迭代次数: {k}) return x_new, k, history x_old x_new print(f警告: 达到最大迭代次数 {max_iter} 仍未收敛最后值: {x_new:.10f}) return x_new, max_iter, history # 使用示例求解 f(x) x^2 - 2 0 (求根号2) import math f_func lambda x: x**2 - 2 f_prime_func lambda x: 2*x x0 1.5 root, iters, hist newton_method(f_func, f_prime_func, x0) print(fsqrt(2)的近似值: {root:.15f}) print(f与math.sqrt(2)的绝对误差: {abs(root - math.sqrt(2)):.2e}) # 可视化牛顿法的几何过程 import numpy as np import matplotlib.pyplot as plt x_vals np.linspace(0.5, 2.5, 400) y_vals f_func(x_vals) plt.figure(figsize(10, 6)) plt.plot(x_vals, y_vals, b-, linewidth2, labelf(x) x^2 - 2) plt.axhline(y0, colork, linestyle-, alpha0.3) # x轴 # 绘制前几次迭代的切线 for i in range(min(4, iters1)): x_k, f_x_k hist[i][1], hist[i][2] if i iters: # 有下一个点才能画切线 slope f_prime_func(x_k) # 切线方程: y f(x_k) f(x_k)*(x - x_k) tangent_line lambda x: f_x_k slope * (x - x_k) # 绘制切线线段范围在x_k附近一小段 x_tangent np.linspace(x_k - 0.5, x_k 0.5, 50) y_tangent tangent_line(x_tangent) plt.plot(x_tangent, y_tangent, r--, linewidth1, alpha0.7) # 标记迭代点 plt.scatter(x_k, f_x_k, colorred, s50, zorder5) plt.text(x_k, f_x_k0.1, f$x_{i}$, fontsize12, hacenter) # 标记切线与x轴交点 (x_{k1}) if i1 len(hist): x_next hist[i1][1] plt.scatter(x_next, 0, colorgreen, s50, zorder5) plt.plot([x_k, x_next], [f_x_k, 0], g:, alpha0.5) # 垂线示意 plt.xlabel(x) plt.ylabel(f(x)) plt.title(牛顿法切线法的几何图示) plt.legend() plt.grid(True, alpha0.3) plt.tight_layout() plt.show()关键实现细节与避坑经验导数安全检测代码中if abs(f_prime_val) eps * 10MATLAB或if abs(f_prime_val) 1e-15Python是生死线。没有这个检查当迭代点接近函数极值点或拐点导数为零时程序会因除以零或极小数而产生Inf、NaN或巨大步长导致崩溃或发散。这个阈值需要根据函数尺度灵活调整。收敛条件牛顿法通常收敛很快所以容差tol可以设得很小如1e-12。停止条件应同时检查迭代步长diff和函数值f(x_new)。有时diff很小但f(x_new)还不够小例如在平台区域所以双重检查更保险。初始值选择牛顿法对初始值敏感。一个实用的策略是先用对分法或试位法等稳健但较慢的方法将根隔离到一个较小区间并得到一个较好的初始近似然后再用牛顿法加速。也可以从不同初始点尝试几次观察收敛情况。处理重根如果根x*是m重根即f(x*)f(x*)...f^{(m-1)}(x*)0, f^{(m)}(x*)≠0标准牛顿法会降级为线性收敛。改进方法是使用修正的牛顿公式x_{k1} x_k - m * f(x_k) / f(x_k)但这需要知道重数m。4. 当导数不可得割线法与拟牛顿法的思路牛顿法需要显式的导数f(x)这有时是个障碍。对于复杂函数、黑箱函数或只能通过实验测量得到函数值的情况我们无法获得解析的导数。这时候割线法提供了一种巧妙的替代方案。割线法的思想是用差商来近似导数。它需要两个初始点x_{k-1}和x_k然后用这两点间的割线斜率(f(x_k) - f(x_{k-1})) / (x_k - x_{k-1})来代替牛顿法中的导数f(x_k)。于是迭代公式变为x_{k1} x_k - f(x_k) * (x_k - x_{k-1}) / (f(x_k) - f(x_{k-1}))从几何上看牛顿法是用切线找下一个点割线法是用过前两点的割线与x轴的交点来找下一个点。割线法的特点优点不需要计算导数只需要函数值。对于计算f(x)成本高昂的问题割线法可能比牛顿法需要额外计算f(x)更高效。收敛速度超线性收敛收敛阶约为(1√5)/2 ≈ 1.618黄金分割率比一般迭代法快但比牛顿法慢。缺点需要两个初始点。和牛顿法一样初始点选择不当可能导致发散。分母f(x_k) - f(x_{k-1})可能接近零需要类似除零保护。MATLAB割线法简单实现function [root, iter] secant_method(f, x0, x1, tol, max_iter) iter 0; fprintf(割线法迭代:\n); fprintf(k\t\t x_k\t\t\t f(x_k)\n); fprintf(-----------------------------------\n); f0 f(x0); f1 f(x1); fprintf(%d\t\t %.10f\t %.10e\n, 0, x0, f0); fprintf(%d\t\t %.10f\t %.10e\n, 1, x1, f1); for iter 2:max_iter if abs(f1 - f0) eps % 防止除零 warning(函数值差过小可能无法继续迭代。); root x1; return; end % 割线法核心公式 x_new x1 - f1 * (x1 - x0) / (f1 - f0); f_new f(x_new); fprintf(%d\t\t %.10f\t %.10e\n, iter, x_new, f_new); if abs(x_new - x1) tol abs(f_new) tol root x_new; fprintf(收敛于迭代次数 %d\n, iter); return; end % 更新点准备下一次迭代 x0 x1; f0 f1; x1 x_new; f1 f_new; end warning(未在最大迭代次数内收敛。); root x1; end在优化和更高维的问题中还有更复杂的拟牛顿法如BFGS、DFP算法它们通过维护一个不断更新的矩阵来近似海森矩阵二阶导数矩阵的逆同样避免了直接计算导数。这是数值优化领域的核心内容但在单变量方程求根中割线法通常已足够好用。5. 实战建模场景如何选择与联用这些方法在数学建模竞赛或实际工程问题中你很少会孤立地使用某一种方法。一个稳健、高效的求根策略往往是多种方法的组合。下面我结合几个典型场景分享我的经验。场景一方程形式简单导数易得且对初始值有粗略估计首选牛顿法。它的平方收敛速度能让你用极少的迭代次数通常5-10次就达到机器精度。例如在物理仿真中求解一个非线性方程来确定某个参数或者在校正模型中求解一个标量方程。操作要点一定要实现上一节提到的导数保护。如果迭代中出现f(x)过小应触发回退策略比如暂时切换到对分法几步或者提示用户检查初始值和函数性质。场景二函数是“黑箱”或求导极其困难选择割线法或不需要导数的迭代法。例如函数f(x)是另一个复杂仿真程序或实验测量的输出你只能得到f(x)的值不知道其解析形式。操作要点割线法需要两个初始点。这两个点应位于根的两侧即f(x0)*f(x1)0这样能保证早期迭代的稳定性。如果无法确定根的两侧可以用两个相近的点开始但需密切监控收敛情况。场景三对稳健性要求极高不求最快但求最稳采用混合策略。这是最推荐的做法尤其是在自动化脚本或通用求解器中。第一步隔离根。使用对分法。这是一个绝对稳健的方法只要找到区间[a, b]满足f(a)*f(b)0它保证收敛。虽然速度慢线性收敛每次区间减半但它不依赖函数性质总能将根的范围缩小到一个很小的区间。第二步精细化求解。将对分法得到的最终区间中点或一个端点作为牛顿法或割线法的初始值。由于此时初始值已经非常接近真根牛顿法能安全、快速地达到高精度。MATLAB中fzero函数的策略MATLAB内置的fzero函数就采用了类似的混合策略。它首先尝试用割线法或逆二次插值法快速逼近如果进展不顺或符号变化会自动切换到对分法以确保稳健性。在Python中scipy.optimize.root_scalar函数提供了多种方法bisect,newton,secant,brentq等其中brentq是结合了对分法、割线法和逆二次插值的混合方法通常是最佳的单变量求根选择。场景四求解多项式方程的全部根对于多项式有更专用的方法如伴随矩阵法MATLAB的roots函数或拉盖尔方法。牛顿法也可用于逐个寻找实根但需要配合多项式降阶每找到一个根就用多项式除法除掉对应的因子(x - root)以避免收敛到已找到的根。一个综合案例求解超越方程f(x) e^{-x} - sin(x) 0这个方程在[0, 2]区间内有根。我们设计一个混合求解流程绘图定位先用MATLAB或Python画出f(x)在[0, 2]的图像肉眼观察根的大概位置约在0.5到0.6之间。import numpy as np; import matplotlib.pyplot as plt x np.linspace(0, 2, 200); y np.exp(-x) - np.sin(x) plt.plot(x, y); plt.axhline(y0, colork); plt.grid(); plt.show()稳健起步由于对函数性质不确定先用对分法在[0.5, 0.6]区间内迭代几次将区间长度缩小到1e-2以下。加速求精取对分法最后区间的中点x0 ≈ 0.55作为牛顿法的初始值。提供导数f(x) -e^{-x} - cos(x)。实现与验证import math def f(x): return math.exp(-x) - math.sin(x) def f_prime(x): return -math.exp(-x) - math.cos(x) # 步骤2: 对分法 (简易实现) def bisection(f, a, b, tol): fa, fb f(a), f(b) if fa * fb 0: raise ValueError(区间两端函数值同号。) while (b - a) / 2 tol: c (a b) / 2 fc f(c) if fc 0: return c if fa * fc 0: b, fb c, fc else: a, fa c, fc return (a b) / 2 # 先用对分法得到一个好起点 root_rough bisection(f, 0.5, 0.6, 1e-3) print(f对分法得到的粗解: {root_rough}) # 步骤3: 用牛顿法精细化 root_refined, iters, _ newton_method(f, f_prime, root_rough, tol1e-12) print(f牛顿法精细化后的解: {root_refined:.15f}) print(f函数值 f(root) {f(root_refined):.2e})这种“稳健方法先行快速方法收尾”的策略在数学建模中非常实用既能避免因初始值差导致的牛顿法失败又能利用牛顿法的高效率获得高精度解。6. 从单变量到多变量思想的延伸虽然本文聚焦于单变量方程求根但无论是迭代法还是牛顿法其思想都可以直接推广到多变量非线性方程组的求解这在数学建模中更为常见例如求解一个包含多个未知数的动力学系统平衡点。对于方程组F(x) 0其中x是向量F是向量值函数不动点迭代的推广形式是x_{k1} G(x_k)其中G是一个向量函数。牛顿法的推广形式成为牛顿-拉夫森方法x_{k1} x_k - J_F(x_k)^{-1} F(x_k)。这里J_F(x_k)是雅可比矩阵一阶偏导数矩阵。核心步骤从标量的除法变成了求解一个线性方程组J_F(x_k) * s -F(x_k)得到步长s。多变量情况下的挑战急剧增加初始猜测更难、雅可比矩阵的计算或近似更复杂、收敛性分析更困难。但万变不离其宗理解单变量情况下的收敛原理、实现细节和 pitfalls是迈向高维问题坚实的基础。最后再分享一个我个人的小技巧在编写任何求根算法时一定要加入详尽的迭代过程输出和可视化。就像我上面代码中做的那样打印出每一步的x_k、f(x_k)和步长。这不仅仅是调试的需要更是你理解算法行为、诊断收敛问题是震荡、发散还是缓慢爬行的最直接工具。一幅像上面那样的几何迭代图能让你对方法的动态过程有直觉上的把握这是任何文字描述都无法替代的。