秦九韶算法与二分法Python实战：算法复杂度可视化验证

简介：本资源是广东工业大学《数值分析》课程Python实验一的完整报告文档，面向数学与计算机相关专业本科生，聚焦算法效率对比与非线性方程求解两大核心数值计算问题。文档以实验报告形式系统呈现：第一部分通过直接法与秦九韶算法计算高次多项式（最高10000次），详列x=0.1/1/2三组测试下的乘法/加法次数及运行时间对比，直观揭示秦九韶算法在大幅降低运算复杂度上的优势；第二部分实现二分法求解x²−2=0，完整记录17次迭代过程、区间缩放逻辑及函数值符号判定规则，并附可运行Python源码。资源为1个29KB的docx文件，结构清晰，含实验目的、算法原理、结果表格、代码实现与深度讨论，便于理解算法思想与编程实践。目前已有211人学习下载，适合数值分析初学者巩固理论、提升Python数值编程能力。

1. 这不是“写个for循环就交作业”的数值实验：它用两个真实可跑的Python脚本，把算法复杂度从“玄学概念”变成终端里跳动的毫秒数和乘法计数器

你有没有在写完一个多项式求值函数后，被导师问过：“x=2时，你这段代码到底执行了多少次乘法？”——不是理论推导，是真·数出来。广东工业大学这份Python数值分析实验一，就是干这个的：它不讲大道理，直接甩给你两段能跑、能改、能测、能对着表格逐行验证的代码——一段是教科书级的“暴力直球”（直接法），一段是古籍里走出来的“降维打击”（秦九韶算法）；再加一个手搓二分法求根的完整迭代日志生成器。它面向的是刚学完for和def、但还没见过“浮点误差怎么吃掉你第15位有效数字”的人；也面向那些想拿它当模板，快速搭起自己数值方法验证框架的开发者。所有代码都在本地Python 3.8+环境一键可跑，不需要任何第三方包，连numpy都不用——因为它的设计哲学很朴素：让计算代价肉眼可见，让算法优劣不再靠背诵，而靠time.perf_counter()打出来的那个小数点后四位。

2. 项目结构与技术栈：纯Python轻量级数值验证套件，零依赖、全手动、重过程

这份实验资源不是工程级项目，而是一个高度聚焦的教学验证套件。它没有requirements.txt，没有setup.py，甚至没有子目录——所有逻辑都压在三个独立.py文件里，对应实验报告中的三大模块。但它对Python语言特性的运用非常扎实：精准使用time.perf_counter()而非time.time()做高精度计时；用len(str(x))粗略估算有效数字位数（虽有浮点表示局限，但教学场景下足够直观）；用open('1.txt', 'a+')持续追加写入二分法迭代轨迹，形成可直接拖进Excel分析的结构化日志。整个技术栈就是Python标准库的“核心四件套”：time（计时）、math（虽未显式import但x*x即其子集）、builtins（print,len,range,input）、io（文件读写）。没有魔法，全是肌肉记忆级别的基础操作——这恰恰是它能被快速复现、深度改造的根本原因。

2.1 文件清单与职责划分：三文件即三战场

整个资源解压后只有3个.py文件，命名直白无歧义：

提示：三个文件彼此完全解耦，可单独运行、单独修改、单独计时。不要试图把它们合并成一个主程序——实验设计的本意就是让你隔离对比，而不是封装调用。

2.2 多项式定义的隐含契约：理解f(x)的数学本体是读懂所有代码的前提

实验中反复出现的多项式并非随意构造，而是有明确数学定义的：

$$ f(x) = 1 + 2x + 3x^2 + 4x^3 + \cdots + 10001x^{10000} = \sum_{k=0}^{10000} (k+1) x^k $$

注意三点关键细节：

• 系数规律：第 $k$ 项（$k$ 从0开始）的系数是 $(k+1)$，不是 $k$。所以常数项是 $1$，一次项是 $2x$，二次项是 $3x^2$，……，最高次项是 $10001x^{10000}$。
• 项数确认：从 $k=0$ 到 $k=10000$，共10001 项。源码中range(10000)循环了10000次，但起始项sum = 1已包含常数项，因此总项数正确。
• 秦九韶变形依据：该多项式可重写为嵌套形式： $$ f(x) = 1 + x\left(2 + x\left(3 + x\left(\cdots + x\left(10000 + x \cdot 10001\right)\cdots\right)\right)\right) $$ 这正是exp2_qin.py中S = x*S + i的数学来源（其中i从10000递减到1，最后一步隐含+1常数项）。

不厘清这个，看代码就会像在读天书——比如看到S = x*S + i时疑惑“为什么i从10000开始？”，答案就藏在这个展开式里。

2.3 核心数据记录机制：如何让“运算次数”从理论走向实测

实验报告里那张漂亮的对比表格（表1），不是靠人脑心算出来的，而是代码里埋好的“计数器”实时吐出的。这是本项目最值得拆解的设计细节：

• 直接法的计数盲区：exp2_direct.py中只用了time.perf_counter()，并未统计乘法/加法次数！它依赖的是Python解释器执行x**(n+1)和sum + ...时底层C实现的运算量，属于“黑盒测量”。这恰恰暴露了一个教学重点：理论复杂度（O(n²)）和实际耗时之间，隔着硬件、编译器优化、幂运算实现等多重黑匣子。
• 秦九韶法的白盒计数：exp2_qin.py则完全不同。它用两个整型变量countMul和countAdd，在每次执行S = x*S + i时手动+1。因为这条语句明确包含1次乘法（x*S）和1次加法（+ i），且循环10000次，所以最终countMul == countAdd == 10000。这是真正的“白盒验证”，把算法步骤和机器指令一一映射。
• 二分法的日志契约：exp3_bisection.py的file_handle.writelines(...)不是随便写的。格式字符串"{:.15g} {:.15g} {:.15g} {:.15g} {:.15g}\n"确保了5列数据（迭代序号、xlow、xup、middle、符号标识）以紧凑、无科学计数法、保留15位有效数字的方式写入，方便后续用pandas.read_csv()或Excel直接识别为数值列。

这种“计数方式”的差异，本身就是一堂关于测量方法论的课：什么时候该信理论推导，什么时候必须亲手埋点？

3. 核心模块逐行解析：从秦九韶的嵌套逻辑到二分法的区间收缩

现在我们沉到代码最深处，一行行拆解这两个最具教学价值的模块。不讲原理，只讲这一行代码在做什么、为什么这么写、改哪里会翻车。

3.1 秦九韶算法实现：exp2_qin.py的12行如何完成10000次降维打击

import time x = 2 powN = 10000 aN = powN + 1 # aN = 10001, 即最高次项系数 countMul = 0 countAdd = 0 length = 0 start = time.perf_counter() S = aN # 初始化 S = 10001 (最高次项系数) for i in range(10000, 0, -1): # i 从 10000 递减到 1 S = x * S + i # 关键：S_{k} = x * S_{k+1} + a_k countAdd += 1 countMul += 1 end = time.perf_counter() length = len(str(S)) print("result=", S) print("乘法次数", countMul) print("加法次数", countAdd) print("length=", length) print("time=%.2g 秒\n" % (end-start))

逐行逻辑说明与参数深挖：

• S = aN：初始化S为最高次项系数10001。这是秦九韶算法的起点，对应公式中最内层括号(... + x * 10001)里的10001。
• for i in range(10000, 0, -1)：循环变量i代表当前要加入的低一次项的系数。当i=10000时，计算的是x*10001 + 10000，即10000x + 10001x²的合并；当i=1时，计算的是x*S_prev + 1，即最终加上常数项1。注意：range(10000, 0, -1)生成的是[10000, 9999, ..., 1]，共10000个数，完美覆盖从x^1到x^{10000}的所有中间系数，常数项1在最后一次+i中加入。
• S = x * S + i：这是整个算法的灵魂。它把高次幂的重复计算，压缩成一次乘法加一次加法。x * S是将上一轮累积的“高次部分”整体左移一位（相当于乘以x），+ i是补上当前阶数的系数。此行执行1次乘法、1次加法，故countMul和countAdd各+1。
• length = len(str(S))：这是教学场景下的取巧方案。对整数S，str(S)返回其十进制字符串，len即位数。但对浮点数或极大数会失效（如S超过10**308会变inf，str(inf)长度为3）。生产环境应改用math.floor(math.log10(abs(S))) + 1，但此处为突出“有效数字”概念，用字符串长度更直观。

血泪经验：初学者常把range写成range(10000, 1, -1)，导致i只到2，漏掉i=1，结果少了常数项1，整个结果偏小。务必记住：range(start, stop, step)的stop是不包含的边界。

3.2 二分法迭代引擎：exp3_bisection.py如何用18行代码跑出17次收敛

a = 2 LIMIT = 1e-20 def f(x): return x*x - a # 求 sqrt(2)，即解 x^2 - 2 = 0 xlow = float(input("请输入 x 值的下限：")) xup = float(input("请输入 x 值的上限：")) iter = 0 while (xup - xlow) * (xup - xlow) > LIMIT: # 收敛条件：区间长度平方 > LIMIT middle = (xup + xlow) / 2 iter += 1 if f(xlow) * f(middle) < 0: # f(xlow) 与 f(middle) 异号 → 根在 [xlow, middle] xup = middle if f(xup) * f(middle) < 0: # f(xup) 与 f(middle) 异号 → 根在 [middle, xup] xlow = middle S1 = (xup + xlow) / 2 S = f(S1) if S > 0: S = 1 elif S < 0: S = 0 else: S = -1 file_handle = open('1.txt', mode='a+') file_handle.writelines("{:.15g} {:.15g} {:.15g} {:.15g} {:.15g}\n".format(iter,xlow,xup,S1,S)) print("{:.15g} {:.15g} {:.15g} {:.15g} {:.15g}".format(iter,xlow,xup,S1,S))

逐行逻辑说明与参数深挖：

• while (xup - xlow) * (xup - xlow) > LIMIT：收敛判断用区间长度的平方而非长度本身。这是为避免浮点数比较时因精度丢失导致死循环。LIMIT = 1e-20是极高的精度要求，意味着最终区间长度< 1e-10（因为(1e-10)^2 = 1e-20）。若改为xup - xlow > LIMIT，在xup和xlow接近1.4142...时，xup - xlow可能因浮点舍入永远无法小于1e-20，导致无限循环。
• if f(xlow) * f(middle) < 0和if f(xup) * f(middle) < 0：这里用乘积判别法替代f(xlow) < 0 and f(middle) > 0等显式比较，是为了规避f(x)恰好为0的边界情况（此时乘积为0，不满足< 0，不会误更新）。但要注意：若f(xlow)或f(middle)为0，则已找到精确根，理想情况应break，但本代码未处理，属教学简化。
• S = f(S1)后的if-elif-else块：这是对f(S1)符号的离散化编码：>0→1，<0→0，=0→-1。写入日志的第五列就是这个编码值，方便后续用脚本统计“正负交替次数”。这不是数学必需，而是为生成可分析日志做的数据工程设计。
• file_handle = open('1.txt', mode='a+')：a+模式是关键。它保证每次运行都追加到文件末尾，不会覆盖历史数据。若用w模式，每次运行都会清空重写，无法积累多组实验数据。这是实验可复现性的基础设施保障。

避坑 / 常见问题 / 排查
现象1：运行exp3_bisection.py时，输入xlow=1.3,xup=1.5后，程序卡住不动，CPU占用率飙升。
原因：LIMIT = 1e-20过于严苛，而float在Python中双精度浮点数的有效位数约15-17位。当xup和xlow逼近sqrt(2)时，(xup - xlow)的差值可能无法被精确表示，导致while条件永远为真。
解决：将LIMIT放宽至1e-13，或改用abs(xup - xlow) > 1e-13并配合math.isclose()做更鲁棒判断。

现象2：1.txt文件里，迭代次数iter从1开始，但第一行日志显示iter=0。
原因：代码中iter += 1在while循环体开头执行，而print和writelines在循环体结尾。但实验报告表2的第一行标为“0”，对应初始输入状态。本代码并未记录初始状态，表2的“0”行是人工添加的。若需自动记录，应在while循环前加一行file_handle.writelines("0 {:.15g} {:.15g} {:.15g} 0\n".format(xlow,xup,(xlow+xup)/2))。

现象3：f(xlow)*f(middle) < 0判断失效，xup和xlow不再收缩，陷入死循环。
原因：f(x)在区间内不连续，或存在多个根，导致符号变化不满足二分法前提（f(xlow)*f(xup) < 0）。本例f(x)=x^2-2在[1.3,1.5]内严格单调，无此问题；但若换成f(x)=x^3-x，在[-2,2]内就有3个根，二分法可能收敛到任意一个。
解决：在while循环前强制校验f(xlow)*f(xup) < 0，不满足则print("Error: f(xlow) and f(xup) must have opposite signs!")并exit()。

现象4：1.txt文件中，S1（即(xup+xlow)/2）的值出现nan或inf。
原因：xup或xlow因浮点溢出变为inf，导致S1计算失败。常见于LIMIT过小且循环次数过多时。
解决：在while循环内加入保护if math.isinf(xup) or math.isinf(xlow): break，或限制最大迭代次数if iter > 1000: break。

4. 本地运行全流程：从环境准备到结果验证，一步不跳过

现在，让我们把前面所有的理解，落地为一次完整的、可复现的本地运行。这不是“复制粘贴就能跑”，而是带着问题意识去执行每一个命令。

4.1 环境准备：确认你的Python版本与基础能力

首先，打开终端（Windows用CMD/PowerShell，macOS/Linux用Terminal），执行：

python --version

确保输出为Python 3.7或更高版本（推荐3.8+）。低于3.7可能缺少time.perf_counter()的高精度支持。接着，确认你能在终端中执行Python脚本：

python -c "print('Hello from Python!')"

如果报错command not found，请先配置Python环境变量（搜索“Windows 添加Python到PATH”即可找到图文教程）。这步看似简单，却是90%新手卡住的第一关。

4.2 创建实验目录并放置代码

在任意位置（如桌面）新建一个文件夹，命名为gdut_numexp1。进入该文件夹：

# Windows PowerShell cd Desktop\gdut_numexp1 # macOS/Linux Terminal cd ~/Desktop/gdut_numexp1

然后，用文本编辑器（推荐 VS Code、Sublime Text 或系统自带记事本）创建三个文件：

• exp2_direct.py：粘贴实验报告中“算法1源代码”全部内容。
• exp2_qin.py：粘贴“算法2源代码”全部内容。
• exp3_bisection.py：粘贴“算法源代码”全部内容（即二分法部分）。

注意：不要从PDF复制！PDF中的空格、破折号（—）可能被转为全角字符，导致Python语法错误。务必手动敲写或从纯文本源复制。

4.3 执行算法对比实验：亲眼见证秦九韶的威力

在gdut_numexp1目录下，依次运行两个多项式脚本。为便于观察，我们先统一测试x=2（报告中耗时最长的case）：

# 运行直接法 python exp2_direct.py # 运行秦九韶法 python exp2_qin.py

你会看到类似这样的输出：

# exp2_direct.py 输出 result= 3015 位有效数字 time=0.142 秒 # exp2_qin.py 输出 result= 3015 位有效数字 乘法次数 10000 加法次数 10000 length= 3015 time=0.033 秒

关键验证点：

• result值是否一致？（应完全相同，证明两种算法数学等价）
• time是否qin < direct？（应显著快于，验证算法优势）
• exp2_qin.py是否输出了乘法次数和加法次数？（应均为10000，验证白盒计数）

提示：若exp2_direct.py的time显示0.0，说明耗时太短，perf_counter()分辨率不够。可将x改为2.0（显式浮点），或在循环内加pass增加负载，但教学目的已达，不必强求。

4.4 执行二分法实验：生成并分析1.txt日志

运行二分法脚本：

python exp3_bisection.py

按提示输入：

请输入 x 值的下限：1.3 请输入 x 值的上限：1.5

程序将打印17行迭代数据，并在当前目录生成1.txt。用文本编辑器打开它，你应该看到：

1 1.4 1.5 1.45 1 2 1.4 1.45 1.425 1 3 1.4 1.425 1.4125 0 ...

验证日志有效性：

• 第一列iter是否从1递增到17？
• 第二、三列xlow/xup是否呈现“区间不断缩小”的趋势？（如1.4→1.4125→1.4140625...）
• 第四列S1是否稳定收敛到1.414213562...（即√2的近似值）？
• 第五列S是否在0和1之间交替？（代表f(S1)符号在正负间切换）

若一切符合，恭喜，你已成功复现了整个实验流程。1.txt就是你亲手生成的、可导入Excel画收敛曲线的原始数据。

5. 二次开发实战：从“跑通”到“为我所用”的5个改造技巧

这份实验代码的价值，远不止于交作业。它是一块打磨数值编程手感的磨刀石。下面这5个改造，都是我在某高校实验室带学生时，被问得最多、也最能提升工程能力的问题。每个都附带可直接粘贴的代码片段和验证方法。

5.1 技巧1：把硬编码的x值改成命令行参数，告别反复改源码

现状：x = 2写死在代码里，想测x=0.1得手动改、保存、再运行。
改造：用sys.argv读取命令行参数。

# 修改 exp2_direct.py 和 exp2_qin.py 的开头部分 import time import sys # 替换原来的 x = 2 if len(sys.argv) > 1: x = float(sys.argv[1]) else: x = 2.0 # 默认值

验证：

python exp2_qin.py 0.1 # 测试 x=0.1 python exp2_qin.py 1 # 测试 x=1

从那以后我每次写数值实验脚本，都强制走一遍sys.argv参数化。因为真实科研中，参数组合可能是x=[0.1,1,2,10],n=[1000,5000,10000]的笛卡尔积，手动改是反人类的。

5.2 技巧2：用matplotlib画出收敛曲线，让“17次迭代”变成一张图

现状：1.txt是纯文本，收敛过程靠肉眼扫。
改造：用pandas读取日志，用matplotlib绘图。

# 新建 plot_convergence.py import pandas as pd import matplotlib.pyplot as plt # 读取日志，指定列名 df = pd.read_csv('1.txt', sep='\s+', header=None, names=['iter', 'xlow', 'xup', 'middle', 'sign']) # 计算每次迭代的区间长度和 middle 与 sqrt(2) 的误差 import math true_root = math.sqrt(2) df['interval_len'] = df['xup'] - df['xlow'] df['error'] = abs(df['middle'] - true_root) # 绘图 fig, ax1 = plt.subplots() color = 'tab:red' ax1.set_xlabel('Iteration') ax1.set_ylabel('Interval Length', color=color) ax1.semilogy(df['iter'], df['interval_len'], color=color, label='Interval Length') ax1.tick_params(axis='y', labelcolor=color) ax2 = ax1.twinx() # 共享x轴 color = 'tab:blue' ax2.set_ylabel('Error to √2', color=color) ax2.semilogy(df['iter'], df['error'], color=color, linestyle='--', label='Error') ax2.tick_params(axis='y', labelcolor=color) fig.tight_layout() plt.title('Bisection Method Convergence') plt.show()

运行：python plot_convergence.py，将弹出一张双Y轴对数图，清晰展示区间长度和误差如何随迭代指数衰减。这张图，比17行数字更有说服力。

5.3 技巧3：把二分法封装成函数，支持任意f(x)和任意区间

现状：exp3_bisection.py只能解x^2-2=0。
改造：抽离为通用函数bisection(f, a, b, tol=1e-10, max_iter=100)。

def bisection(f, a, b, tol=1e-10, max_iter=100): """ 通用二分法求根 :param f: 目标函数，接受单个float参数 :param a, b: 初始区间端点，要求 f(a)*f(b) < 0 :param tol: 收敛容差（区间长度） :param max_iter: 最大迭代次数 :return: (root, iter_count, history_list) """ if f(a) * f(b) >= 0: raise ValueError("f(a) and f(b) must have opposite signs") history = [] for i in range(max_iter): c = (a + b) / 2 fc = f(c) history.append((i+1, a, b, c, 1 if fc > 0 else 0 if fc < 0 else -1)) if abs(b - a) < tol: return c, i+1, history if f(a) * fc < 0: b = c else: a = c return (a + b) / 2, max_iter, history # 使用示例：解 sin(x) = 0 在 [3, 4] 内的根（即 π） import math root, iters, log = bisection(lambda x: math.sin(x), 3, 4) print(f"Root ≈ {root}, found in {iters} iterations")

从那以后我每次遇到新方程，第一反应不是写新脚本，而是调用这个bisection()函数。它把算法逻辑和业务逻辑彻底分离，这才是工程思维的起点。

5.4 技巧4：为秦九韶算法增加“系数数组”输入，让它不再只认1,2,3,...,10001

现状：exp2_qin.py的系数是硬编码的递增序列。
改造：接受任意系数列表coeffs = [a0, a1, a2, ..., an]。

def qin_shao(coeffs, x): """ 秦九韶算法通用实现 :param coeffs: 系数列表，coeffs[i] 是 x^i 的系数 :param x: 自变量值 :return: 多项式在 x 处的值 """ if not coeffs: return 0 S = coeffs[-1] # 从最高次项系数开始 # 从倒数第二项开始，向前遍历 for i in range(len(coeffs)-2, -1, -1): S = x * S + coeffs[i] return S # 测试：计算 f(x) = 2 + 3x + 5x^2 在 x=2 处的值 coeffs = [2, 3, 5] # a0=2, a1=3, a2=5 result = qin_shao(coeffs, 2) # 应等于 2 + 3*2 + 5*4 = 28 print(result) # 28

验证：qin_shao([1,2,3,4], 1)应等于1+2+3+4=10；qin_shao([1,0,0,1], 2)应等于1 + 0*2 + 0*4 + 1*8 = 9。这个函数，让你能把任何多项式丢给秦九韶，而不必再推导嵌套公式。

5.5 技巧5：用pytest写单元测试，让“结果正确”不再是口头承诺

现状：正确性靠肉眼比对result=后的数字。
改造：为关键函数写自动化测试。

# test_numexp.py import pytest from exp2_qin import qin_shao # 假设你已将 qin_shao 提取为独立函数 def test_qin_shao_constant(): assert qin_shao([5], 10) == 5 def test_qin_shao_linear(): assert qin_shao([1, 2], 3) == 1 + 2*3 # 7 def test_qin_shao_quadratic(): assert qin_shao([1, -2, 1], 2) == 1 - 2*2 + 1*4 # 1 def test_bisection_sqrt2(): from exp3_bisection import bisection root, _, _ = bisection(lambda x: x*x - 2, 1.3, 1.5) assert abs(root - 1.41421356237) < 1e-10

运行：pytest test_numexp.py -v。如果所有测试通过，你就可以确信：无论代码怎么改，核心逻辑的正确性已被机器守护。希望帮到你。

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