数值分析上机：秦九韶算法与二分法求根的Python实现及性能对比

简介：这是广东工业大学数值分析课程的Python上机实验报告，围绕“减少运算次数”与“求解非线性方程”两个主题展开。前一部分对比直接法与秦九韶算法在计算同一多项式时的表现，通过录制Python源代码、统计乘法/加法次数及耗时，验证了秦九韶算法在较大数值运算上的显著优势；后一部分实现二分法求方程根，记录了每次迭代的区间上下限、中间点函数值及正负性，展示逐步逼近根值的过程。资源以docx格式呈现，压缩包内仅含1个文档，整体约29KB，轻量便于使用。目前已有211人学习下载。报告除了给出完整可运行的算法代码、实验结果表与结论讨论外，还包含了程序注释和输出示例，适合计算机或数学相关专业学生对照复现，深入理解算法效率分析与数值逼近思想，也可作为实验报告撰写的结构参考。

1. 数值分析上机第一回：同是一万次多项式求值，秦九韶把乘法次数从五千万压到一万

数值分析第一课上机实验，拿到题目时很多人第一反应都是 for 循环直接累加 x 的高次幂，跑完回来看统计数字才意识到自己跟“暴力求解”绑定了多久。我这次拿到的模拟项目X要求用 Python 算一个一万次多项式的值，并分别记录直接法和秦九韶算法的乘法、加法次数以及耗时。结果很直观：直接法理论乘法次数接近五千万次，而秦九韶算法只用了 10000 次乘法和 10000 次加法，x 越大优势越明显。文章把这两套代码、二分法求根的完整实现，以及我调试时踩过的几个坑原样拆给你。适合正在写数值分析实验报告、想验证算法效率差异，或者需要导出大量迭代数据的同学。

2. 多项式求值两种算法：直接法与秦九韶，差距不在代码量而在乘法次数

2.1 直接法：暴力循环累加 x 的高次幂

实验要求计算的多项式是：

f(x) = 1 + 2x + 3x^2 + ... + 10001x^10000

也就是最高次数 n=10000，系数从 2 到 10001 按顺序递增，再加上常数项 1。直接法的思路就是循环中一项一项加，每次用幂运算x ** (n+1)得到对应的次幂，再乘上系数。

原始代码很短，但“省事”并不等于“省运算”。先看第一版：

import time x = 2 sum_val = 1 start = time.perf_counter() for n in range(10000): sum_val = sum_val + (n+2) * (x ** (n+1)) length = len(str(sum_val)) end = time.perf_counter() print("result=", sum_val) print("有效数字位数", length) print("time=%.2g 秒\n" % (end-start))

这段代码能跑，但有两个明显问题：一是x ** (n+1)每次都把 x 的幂从头算起，前面算过的中间结果完全没有复用；二是len(str(sum_val))统计出来的并不是有效数字位数，而是数字字符串总长度，对大整数来说这个值会随数位增长，但跟“有效数字”概念不是一回事。

为了跟秦九韶公平对比，我重写一版直接法，加上运算次数计数器：

import time x = 2 sum_val = 1 count_mul = 0 count_add = 0 start = time.perf_counter() for n in range(10000): term = (n+2) * (x ** (n+1)) count_mul += 1 count_add += 1 sum_val += term count_add += 1 end = time.perf_counter() print("result=", sum_val) print("乘法次数", count_mul) print("加法次数", count_add) print("time=%.2g 秒\n" % (end-start))

我这里的计数器只统计了显式的乘法和加法，没有拆开幂运算内部的乘法。如果按最朴素的理解，把x ** k看成 k 次连乘，那么从 n=0 到 n=9999，所有幂运算叠加的乘法总次数约等于：

sum_{k=1}^{10000} k = 50005000

再加上每轮循环里一次系数乘法，总乘法次数接近 5000 万。实验报告里记录的算法 1 乘法次数 50006000，也正是拿这个口径统计的。所以直接法的问题不在代码能不能跑，而在每一次幂运算都把之前的“增量”重新算了一遍。

2.2 秦九韶算法：把一万次乘方拆成一万次乘加

秦九韶算法的核心是改写多项式形式，把嵌套的乘法一层层剥出来。对 f(x) = a_n x^n + a_{n-1} x^{n-1} + ... + a_0，可以写成：

f(x) = (((a_n x + a_{n-1}) x + a_{n-2}) x + ... ) x + a_0

放到这个具体多项式里，因为 a_i 恰好是下标加 1，所以计算顺序就变成了从最高次向低次叠加系数。实现代码很简洁：

import time x = 2 pow_n = 10000 a_n = pow_n + 1 # 最高次项系数，即 10001 count_mul = 0 count_add = 0 start = time.perf_counter() S = a_n for i in range(pow_n, 0, -1): S = x * S + i count_mul += 1 count_add += 1 end = time.perf_counter() length = len(str(S)) print("result=", S) print("乘法次数", count_mul) print("加法次数", count_add) print("length=", length) print("time=%.2g 秒\n" % (end-start))

循环从 i=10000 递减到 1，初始 S 是 10001。第一次循环 S = x * 10001 + 10000，得到 10001x + 10000；第二次 S = x*(10001x+10000) + 9999，得到 10001x^2 + 10000x + 9999。经过 10000 次迭代后，S 恰好等于 10001x^10000 + 10000x^9999 + ... + 2x + 1，和直接法的结果完全一致。

参数说明：pow_n是多项式最高次数，这里等于 10000；a_n是最高次系数，等于pow_n + 1；range(pow_n, 0, -1)生成从 10000 到 1 的整数序列，结束值 0 不会被包含，所以最后一个 i 是 1。循环内每次只做一次乘法和一次加法，运算次数被严格控制在 n 量级。

这就是“减少运算次数”的关键：秦九韶完全不计算任何 x^k，而是让 x 在每轮乘法中自然累积幂次。n 次多项式只需要 n 次乘法和 n 次加法，复杂度从 O(n^2) 级别降到了 O(n)。

2.3 运行结果对比：x 越大，秦九韶优势越明显

实验原文测了 x=0.1、1、2 三个点，统计结果整理成下面的表：

注意 x=1 时，多项式变成 1 + 2 + 3 + ... + 10001，结果应该就是 50015001，但表格里的结果在报告中没写全，这里不纠结精确值。重点是两种算法算出的函数值完全一致，但运算次数和耗时差异巨大。

直接法比秦九韶多出来的时间几乎都耗在幂运算上。x=0.1 的时候，小数幂运算量不大，两者差距只有 0.0028 秒；到了 x=2，大整数幂运算让直接法耗时要 0.142 秒，秦九韶只用了 0.033 秒，快了约 4.3 倍。这个现象很典型：输入数值越大，幂运算的中间结果膨胀越厉害，暴力循环的开销会非线性上升，而秦九韶始终维持线性的乘加次数。

实验报告的结论是“算法 2 越有益于较大数值运算”，这句话放到生产环境同样成立。当你处理的不是一万次多项式，而是几十万次或带大系数时，直接法会慢到无法接受，秦九韶的线性复杂度反而是唯一的可行解。

3. 二分法求 sqrt(2)：从迭代表看懂区间收敛，顺便解决“数据往 Word 贴”的麻烦

3.1 二分法原理：零点存在性定理与迭代终止条件

二分法依赖的是闭区间上连续函数的零点存在性定理：如果 f(a) 和 f(b) 异号，那么 (a,b) 内至少有一个根。每次取中点 m = (a+b)/2，根据 f(m) 的符号决定把区间缩小到左半还是右半。每迭代一次，区间长度减半，n 次迭代后精度理论上能达到 (b-a)/2^n。

实验里要求解的是方程 x^2 - 2 = 0 的正根，也就是 sqrt(2)。初始区间取 [1.3, 1.5]，f(1.3) = 1.69 - 2 = -0.31 < 0，f(1.5) = 2.25 - 2 = 0.25 > 0，满足异号条件，所以可以放心迭代。

原报告里还打印了一张迭代表，记录每一步的下限、上限、中点和中点的函数值符号。其实这就是一个“观察算法在动”的好方法。我第一次跑的时候只看了最终结果，根本不知道区间是怎么收缩的；后来把日志打印出来才发现，中点的符号变化非常规律。

3.2 代码实现：带迭代日志的二分法

原文代码有两个地方不太稳妥，一是用了a+写入模式，二是没有处理 f(middle) == 0 的情况。我重新写了一个更安全的版本，同时把日志输出从控制台改成文件：

import time a = 2 LIMIT = 1e-20 def f(x): return x * x - a xlow = float(input("请输入 x 值的下限：")) xup = float(input("请输入 x 值的上限：")) if f(xlow) * f(xup) > 0: raise ValueError("f(xlow) 和 f(xup) 必须异号") iter_num = 0 start = time.perf_counter() with open('bisection_log.tsv', 'w', encoding='utf-8') as log: log.write("iter\txlow\txup\tmid\tsign\n") while (xup - xlow) * (xup - xlow) > LIMIT: middle = (xup + xlow) / 2 iter_num += 1 if f(xlow) * f(middle) < 0: xup = middle elif f(xup) * f(middle) < 0: xlow = middle else: # f(middle) == 0，正好命中根 break sign = 1 if f(middle) > 0 else (-1 if f(middle) < 0 else 0) log.write("{}\t{:.15g}\t{:.15g}\t{:.15g}\t{}\n".format( iter_num, xlow, xup, middle, sign)) print("{}\t{:.15g}\t{:.15g}\t{:.15g}\t{}".format( iter_num, xlow, xup, middle, sign)) result = (xlow + xup) / 2 print("最终根 = {:.15g}".format(result)) print("迭代次数 =", iter_num) print("耗时 = %.4f 秒" % (time.perf_counter() - start))

逻辑说明：while (xup-xlow)*(xup-xlow) > LIMIT比较的是区间长度的平方，所以当长度小于 1e-10 时就会停止。LIMIT设成 1e-20 意味着理论精度约 1e-10 位，已经足够展示收敛过程了。如果设得更小，比如 1e-30，在 Python 浮点数下其实也达不到，因为 double 的精度极限就在那里。

if f(xlow)*f(middle) < 0判断根在左半区间，所以把上限更新为 middle；elif f(xup)*f(middle) < 0判断根在右半区间，把下限更新为 middle；如果两者都不满足，说明 f(middle) 等于 0，直接跳出循环。注意这里用了elif，不是两个独立的if，否则 f(middle) 恰好为 0 时会什么都不做，循环卡死。

这个代码还会把每次迭代的记录写入bisection_log.tsv，我用的是'w'模式，每次运行都会清空旧文件，避免上次的日志混进新报告里。输出格式用了{:.15g}，这样对浮点数自动保留 15 位有效数字，不会把 1.4142135623730951 这种长尾巴完全打出来难看。

运行这一段，前几步输出大概是：

随着迭代继续，上下限会慢慢向 sqrt(2) 附近挤压。到第 17 次左右，区间上下限已经到 1.4142 附近，中点的前四位小数基本稳定。

3.3 附加题：把大量迭代数据写进 txt 和 Excel，解决“贴 Word 贴到手软”

实验报告最后有一道附加题，问“大量输出数据怎么方便地贴到 Word 或者直接写成矩阵到 txt”。原文答案是open()加a+模式。这个思路没问题，但实际写实验报告时，更好的做法是直接生成制表符分隔的 tsv 文件，或者用 pandas 输出 Excel。这样 Word 和 Excel 都能直接识别成表格。

如果不想装 pandas，标准库的 csv 模块就能写出一个规整的表格文件：

import csv rows = [] with open('bisection_log.tsv', 'w', encoding='utf-8', newline='') as f: writer = csv.writer(f, delimiter='\t') writer.writerow(['iter', 'xlow', 'xup', 'mid', 'sign']) # 循环时用 writer.writerow([iter, xlow, xup, middle, sign])

如果机器上有 pandas，可以一行代码把数据变成 Excel：

import pandas as pd df = pd.DataFrame(rows, columns=['iter', 'xlow', 'xup', 'mid', 'sign']) df.to_excel('bisection_log.xlsx', index=False)

rows可以在二分法循环里 append，每个元素是一个列表，代表一行数据。这样导出的 Excel 可以直接塞进实验报告的附录，不需要手动复制表格。

TSV 文件用\t分隔，复制进 Word 后选中文本，点击“插入表格”就能转成表格。这也是我后来写数值分析报告的习惯：任何超过 10 行的输出，一律写成 tsv 或 Excel，绝不手贴。

4. 避坑：数值分析上机实验里我踩过的五个典型问题

4.1 现象：把 LIMIT 设成 1e-20，程序跑半天不结束

有次为了追求“高精度”，我把二分法的LIMIT改成1e-20，然后看着迭代次数一直涨，输出结果却不再变化。原因是 Python 浮点数只有约 15 到 17 位十进制精度，当区间长度缩到 1e-14 左右时，xup - xlow已经无法通过浮点数继续缩小，(xup-xlow)*(xup-xlow)反而可能因为舍入误差停在一个大于 LIMIT 的值上，于是 while 循环永远为真，程序进入死循环。

解决方法是把终止条件设成区间长度和函数值双重判断，比如：

while (xup - xlow) > 1e-10 and abs(f(middle)) > 1e-10:

不要迷信 LIMIT 越小越好。对普通实验，1e-12就足够了，再小只会让迭代卡在浮点噪音里。

4.2 现象：二分法直接把根输出到区间之外，或者根本不对

有次我把区间更新写成两个独立的if，忘记加elif：

if f(xlow) * f(middle) < 0: xup = middle if f(xup) * f(middle) < 0: xlow = middle

当f(middle)恰好等于 0 时，两个条件都为 False，xlow和xup都不更新，循环变成无限循环。就算f(middle)不为 0，两个条件也只有一个成立，所以两个if不会同时触发，但万一浮点数误差导致两个都成立，上下限会同时赋值，区间反而变大。

解决方法是改成if ... elif ... else: break，并且在初始时校验f(xlow)*f(xup) < 0是否成立。我一般会在代码里直接抛出ValueError，否则区间外求解的结果看着像根，实际完全不对。

4.3 现象：用a+模式写记录文件，第二次运行后文件里全是旧数据

原实验代码里用了open('1.txt', mode='a+')，这是追加模式。第一次运行没问题，第二次运行会把新日志接在旧日志后面，文件越来越大，而实验报告要求只贴本次迭代的数据，结果我差点把 17 次迭代的历史记录和上一次的结果混在一起贴上去。

解决方法是改用'w'模式，每次运行覆盖写入。如果确实想保留历史，可以手动加时间戳：

import time filename = f"bisection_log_{time.strftime('%Y%m%d_%H%M%S')}.tsv"

另外建议用with open(...) as f而不是裸open()，这样不用记着f.close()，也不用担心中途报错导致文件句柄泄漏。

4.4 现象：len(str(sum_val))报出来的“有效数字位数”不对

实验代码里有一行：

length = len(str(sum_val))

当时以为能算出有效数字位数，结果发现对大整数来说，它返回的是十进制位数，不是有效数字位数。比如结果等于 50015001，len(str(...))得到 8，但这个数的有效数字就是 8 位吗？不一定。有效数字要求从第一个非零数字开始计算，末尾的 0 算有效，中间和前面的 0 不算，而str长度会把整个数的所有位都数进去。

正确做法是判断浮点数结果时用format(sum_val, '.15g')，或者用 Decimal 模块：

from decimal import Decimal d = Decimal(sum_val) print(d.adjusted()) # 有效数字位数

不过对实验报告而言，最简单的方式是直接把结果格式化保留 15 位有效数字，不要去纠结“位数”。我又一次意识到，实验代码里的注释写“有效数字位数”，其实根本不准确，这种小地方最容易在答辩时被老师问住。

4.5 现象：二分法输入上下限时，f(a) 和 f(b) 同号，程序照样跑，结果收敛到假根

二分法看着简单，但边界条件没想清楚就会翻车。有次我测试区间 [1, 2] 求 x^2 - 2 = 0，没问题。后来我故意把区间改成 [1.4, 1.5]，f(1.4) = -0.04，f(1.5) = 0.25，异号，正常。但如果输入 [1, 1.5]，f(1) = -1，f(1.5) = 0.25，还是异号。最怕的是输入 [1, 1.2]，两个函数值都小于 0，程序照跑，但迭代完全取决于初始的符号判断，最后可能收敛到区间边缘，输出一个 f(x) 不等于 0 的假根。

解决方法是进入循环前强制校验：

if f(xlow) * f(xup) > 0: raise ValueError("区间端点的函数值必须异号")

加上这行后，不需要肉眼去看结果是否正确，程序直接拒绝进入不合理的状态。从那以后我每次写数值实验，都会先检查终止条件、输入校验和文件写入模式，确保这些基础问题不再重演。

5. 进阶：用 matplotlib 把收敛过程和性能对比变成能放进报告的两张图

数值分析实验报告除了要求贴表格和代码，还希望展示“算法在动”。文字描述“区间不断缩小”很抽象，不如直接画图。下面两张图是我后来做同类实验时固定的套路，既能展示二分法的收敛速度，也能对比两种多项式求值算法的耗时差异。

第一张图是二分法区间长度随迭代次数的变化曲线。把每次迭代后的xup - xlow存下来，画出来就可以看到区间长度呈指数下降。注意纵轴要用对数刻度，否则从 0.2 到 1e-10 的跨度根本看不清：

import matplotlib.pyplot as plt # lengths 是每次迭代后的区间长度列表 plt.plot(range(len(lengths)), lengths, marker='o') plt.yscale('log') plt.xlabel('迭代次数') plt.ylabel('区间长度') plt.title('二分法区间收敛过程') plt.savefig('bisection_convergence.png', dpi=150) plt.show()

这段代码里lengths需要你在二分法循环里追加xup - xlow。画出来之后，曲线会是一条快速下坠的折线，到 10 次迭代左右就接近 1e-4，非常直观。

第二张图是直接法和秦九韶算法在三个 x 值下的耗时对比。因为直接法在 x=2 时耗时明显偏长，普通柱状图会让小数值对应的柱子几乎看不见，所以同样用对数刻度：

import matplotlib.pyplot as plt x_vals = ['x=0.1', 'x=1', 'x=2'] direct = [0.0071, 0.0131, 0.142] qin = [0.0043, 0.0065, 0.033] bar_width = 0.35 index = range(len(x_vals)) plt.bar(index, direct, bar_width, label='直接法') plt.bar([i + bar_width for i in index], qin, bar_width, label='秦九韶') plt.yscale('log') plt.xticks([i + bar_width / 2 for i in index], x_vals) plt.ylabel('耗时(秒)') plt.title('多项式求值耗时对比') plt.legend() plt.savefig('polynomial_time_compare.png', dpi=150) plt.show()

这里用plt.yscale('log')之后，x=0.1 和 x=1 的时间差距就不会被 x=2 压扁。如果你用的是中文，需要在代码开头加plt.rcParams['font.sans-serif'] = ['SimHei']，否则标题会乱码。我后来每次都先把字体设置写上，已经成为习惯。

画完图再填进实验报告，整个文档的观感立刻提升一个档次。相比表格里的原始数字，评阅老师更愿意看这张“区间长度下降”的图，因为它一眼就证明二分法在收敛，而不是你手动编了一堆数据。从那以后我每次数值分析实验都会强制自己跑完数据后顺手画图，顺便把迭代日志存成 tsv，两件事加起来不超过五分钟，却能让报告少挨很多质疑。希望帮到你。

本文还有配套的精品资源，点击获取