数值分析上机：秦九韶算法与二分法求根的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也就是最高次数 n10000系数从 2 到 10001 按顺序递增再加上常数项 1。直接法的思路就是循环中一项一项加每次用幂运算x ** (n1)得到对应的次幂再乘上系数。原始代码很短但“省事”并不等于“省运算”。先看第一版import time x 2 sum_val 1 start time.perf_counter() for n in range(10000): sum_val sum_val (n2) * (x ** (n1)) length len(str(sum_val)) end time.perf_counter() print(result, sum_val) print(有效数字位数, length) print(time%.2g 秒\n % (end-start))这段代码能跑但有两个明显问题一是x ** (n1)每次都把 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 (n2) * (x ** (n1)) 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 次连乘那么从 n0 到 n9999所有幂运算叠加的乘法总次数约等于sum_{k1}^{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))循环从 i10000 递减到 1初始 S 是 10001。第一次循环 S x * 10001 10000得到 10001x 10000第二次 S x*(10001x10000) 9999得到 10001x^2 10000x 9999。经过 10000 次迭代后S 恰好等于 10001x^10000 10000x^9999 ... 2x 1和直接法的结果完全一致。参数说明pow_n是多项式最高次数这里等于 10000a_n是最高次系数等于pow_n 1range(pow_n, 0, -1)生成从 10000 到 1 的整数序列结束值 0 不会被包含所以最后一个 i 是 1。循环内每次只做一次乘法和一次加法运算次数被严格控制在 n 量级。这就是“减少运算次数”的关键秦九韶完全不计算任何 x^k而是让 x 在每轮乘法中自然累积幂次。n 次多项式只需要 n 次乘法和 n 次加法复杂度从 O(n^2) 级别降到了 O(n)。2.3 运行结果对比x 越大秦九韶优势越明显实验原文测了 x0.1、1、2 三个点统计结果整理成下面的表x算法函数结果乘法次数加法次数用时(秒)0.1直接法1.234567901234568约 50005000约 300000.00710.1秦九韶1.23456790123456810000100000.00431直接法与理论一致约 50005000约 300000.01311秦九韶与理论一致10000100000.00652直接法大整数约 3015 位约 50005000约 300000.1422秦九韶大整数约 3015 位10000100000.033注意 x1 时多项式变成 1 2 3 ... 10001结果应该就是 50015001但表格里的结果在报告中没写全这里不纠结精确值。重点是两种算法算出的函数值完全一致但运算次数和耗时差异巨大。直接法比秦九韶多出来的时间几乎都耗在幂运算上。x0.1 的时候小数幂运算量不大两者差距只有 0.0028 秒到了 x2大整数幂运算让直接法耗时要 0.142 秒秦九韶只用了 0.033 秒快了约 4.3 倍。这个现象很典型输入数值越大幂运算的中间结果膨胀越厉害暴力循环的开销会非线性上升而秦九韶始终维持线性的乘加次数。实验报告的结论是“算法 2 越有益于较大数值运算”这句话放到生产环境同样成立。当你处理的不是一万次多项式而是几十万次或带大系数时直接法会慢到无法接受秦九韶的线性复杂度反而是唯一的可行解。3. 二分法求 sqrt(2)从迭代表看懂区间收敛顺便解决“数据往 Word 贴”的麻烦3.1 二分法原理零点存在性定理与迭代终止条件二分法依赖的是闭区间上连续函数的零点存在性定理如果 f(a) 和 f(b) 异号那么 (a,b) 内至少有一个根。每次取中点 m (ab)/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 0f(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, encodingutf-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判断根在左半区间所以把上限更新为 middleelif f(xup)*f(middle) 0判断根在右半区间把下限更新为 middle如果两者都不满足说明 f(middle) 等于 0直接跳出循环。注意这里用了elif不是两个独立的if否则 f(middle) 恰好为 0 时会什么都不做循环卡死。这个代码还会把每次迭代的记录写入bisection_log.tsv我用的是w模式每次运行都会清空旧文件避免上次的日志混进新报告里。输出格式用了{:.15g}这样对浮点数自动保留 15 位有效数字不会把 1.4142135623730951 这种长尾巴完全打出来难看。运行这一段前几步输出大概是迭代次数下限上限中点f(中点)符号11.31.51.4负21.41.51.45正31.41.451.425正41.41.4251.4125负51.41251.4251.41875正随着迭代继续上下限会慢慢向 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, encodingutf-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可以一行代码把数据变成 Excelimport pandas as pd df pd.DataFrame(rows, columns[iter, xlow, xup, mid, sign]) df.to_excel(bisection_log.xlsx, indexFalse)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忘记加elifif f(xlow) * f(middle) 0: xup middle if f(xup) * f(middle) 0: xlow middle当f(middle)恰好等于 0 时两个条件都为 Falsexlow和xup都不更新循环变成无限循环。就算f(middle)不为 0两个条件也只有一个成立所以两个if不会同时触发但万一浮点数误差导致两个都成立上下限会同时赋值区间反而变大。解决方法是改成if ... elif ... else: break并且在初始时校验f(xlow)*f(xup) 0是否成立。我一般会在代码里直接抛出ValueError否则区间外求解的结果看着像根实际完全不对。4.3 现象用a模式写记录文件第二次运行后文件里全是旧数据原实验代码里用了open(1.txt, modea)这是追加模式。第一次运行没问题第二次运行会把新日志接在旧日志后面文件越来越大而实验报告要求只贴本次迭代的数据结果我差点把 17 次迭代的历史记录和上一次的结果混在一起贴上去。解决方法是改用w模式每次运行覆盖写入。如果确实想保留历史可以手动加时间戳import time filename fbisection_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))当时以为能算出有效数字位数结果发现对大整数来说它返回的是十进制位数不是有效数字位数。比如结果等于 50015001len(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.04f(1.5) 0.25异号正常。但如果输入 [1, 1.5]f(1) -1f(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, markero) plt.yscale(log) plt.xlabel(迭代次数) plt.ylabel(区间长度) plt.title(二分法区间收敛过程) plt.savefig(bisection_convergence.png, dpi150) plt.show()这段代码里lengths需要你在二分法循环里追加xup - xlow。画出来之后曲线会是一条快速下坠的折线到 10 次迭代左右就接近 1e-4非常直观。第二张图是直接法和秦九韶算法在三个 x 值下的耗时对比。因为直接法在 x2 时耗时明显偏长普通柱状图会让小数值对应的柱子几乎看不见所以同样用对数刻度import matplotlib.pyplot as plt x_vals [x0.1, x1, x2] 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, dpi150) plt.show()这里用plt.yscale(log)之后x0.1 和 x1 的时间差距就不会被 x2 压扁。如果你用的是中文需要在代码开头加plt.rcParams[font.sans-serif] [SimHei]否则标题会乱码。我后来每次都先把字体设置写上已经成为习惯。画完图再填进实验报告整个文档的观感立刻提升一个档次。相比表格里的原始数字评阅老师更愿意看这张“区间长度下降”的图因为它一眼就证明二分法在收敛而不是你手动编了一堆数据。从那以后我每次数值分析实验都会强制自己跑完数据后顺手画图顺便把迭代日志存成 tsv两件事加起来不超过五分钟却能让报告少挨很多质疑。希望帮到你。本文还有配套的精品资源点击获取