平方根求解-趋近法步长保守策略数学推导

平方根求解 — 趋近法步长保守策略数学推导

预计阅读时间：4 分钟

目录

• 问题描述
• 保守步长选择推导
• 实现代码
• 注

问题描述

循环退出要求的是函数值误差满足：

$$|x^2-c|\le \epsilon$$

令真实解为 $r=\sqrt{c}$，设当前估计值为 $x=r+\Delta$，其中 $\Delta$ 是自变量偏差（即 $|x-r|=|\Delta |$）。

则有：

$$|x^2-c|=|(r+\Delta {)}^2-r^2|=|2r\Delta +{\Delta }^2|\le 2r|\Delta |+{\Delta }^2$$

如果误差 $\le \epsilon$，并且假设步长很小（即 $|\Delta |\ll 1$），可忽略 ${\Delta }^2$ 项，得到近似条件：

$$2r|\Delta |\lesssim \epsilon \Rightarrow |\Delta |\lesssim \frac{\epsilon }{2r}$$

由于趋近法的最大可能偏差就是步长 $\delta$（因为你只能停在某个 $x=k\cdot \delta$ 上，而真实解可能在两个步长之间），因此为保证 $|\Delta |<\delta \le \frac{\epsilon }{2r}$，必须取：

$$\delta \le \frac{\epsilon }{2\sqrt{c}}$$

但这里有个问题：不能在不知道 $\sqrt{c}$ 的前提下使用 $\sqrt{c}$ —— 这是要计算的目标。

保守步长选择推导

因此，为了构造一个可计算的上界，用一个已知且大于 $\sqrt{c}$ 的量来代替分母中的 $\sqrt{c}$。一个简单且安全的（但略松的）选择是：

$$\sqrt{c}\le \max (1,c)$$

更合理的做法是直接定义一个无需开方就能计算的上界 $U\ge \sqrt{c}$，例如：

• 若 $c\le 1$，则 $\sqrt{c}\le 1$；
• 若 $c>1$，则 $\sqrt{c}<c$。

因此，可以取

$$U=\max (1,c)$$

于是，保守的步长选择为：

$$\delta =\frac{\epsilon }{2U}=\frac{\epsilon }{2\cdot \max (1,c)}$$

注

上界 $U$ 的选择可以根据具体问题更精细地调整（例如可以基于初值估计或区间界来取更紧的上界），这里给出的是一个简单、通用且安全的保守估计。

实现代码

如下为具体实现。

import time
import csv
import math
import mpmathmpmath.mp.dps = 50  # 设置高精度def linear_approximation_sqrt(target, precision=0.0001):"""使用线性逼近法计算平方根。从初始猜测开始，每次增加步长直到足够接近。"""start_time = time.perf_counter_ns()step = precision / (2 * max(1, target)) # 保守策略iterations = 0approx = 0for candidate in range(0, int(target) + 1):if candidate * candidate > target and approx == 0:approx = candidate - 1records = []power = 0max_iterations = 1_000_000_000  # 增加最大迭代限制while abs(approx * approx - target) > precision and iterations < max_iterations:approx += stepiterations += 1if iterations == 2**power:records.append((iterations, approx, time.perf_counter_ns() - start_time))power += 1end_time = time.perf_counter_ns()time_ns = end_time - start_timereturn approx, iterations, time_ns, recordsdef binary_search_sqrt(target, precision=1e-13):"""使用二分查找法计算平方根。不断缩小范围直到达到精度。"""start_time = time.perf_counter_ns()iterations = 0lower_bound = 0.0upper_bound = max(1.0, float(target)) # max(1.0, target) 确保上界至少为 1approx = (lower_bound + upper_bound) / 2.0records = []# 为防止由于数值精度或错误的区间导致无限循环，增加一个保守的停止条件max_iterations = 10_000_000while abs(approx * approx - target) > precision and iterations < max_iterations:records.append((iterations, approx, time.perf_counter_ns() - start_time))if approx * approx < target:lower_bound = approxelse:upper_bound = approxprev_approx = approxapprox = (lower_bound + upper_bound) / 2.0iterations += 1# 如果中点不再改变（浮点原因），退出避免无限循环if approx == prev_approx:breakend_time = time.perf_counter_ns()time_ns = end_time - start_timereturn approx, iterations, time_ns, recordsdef newton_raphson_sqrt(target, precision=1e-13):"""使用牛顿-拉夫逊法计算平方根。使用公式迭代改进近似值。"""start_time = time.perf_counter_ns()approx = target / 2iterations = 0records = []while abs(approx * approx - target) > precision:records.append((iterations, approx, time.perf_counter_ns() - start_time))approx = (approx + target / approx) / 2iterations += 1end_time = time.perf_counter_ns()time_ns = end_time - start_timereturn approx, iterations, time_ns, records# 实验矩阵
experiments = [("1-A", 0.1234, 1e-3),("1-B", 0.1234, 1e-6),("1-C", 0.1234, 1e-8),("2-A", 10.0, 1e-3),("2-B", 10.0, 1e-6),("2-C", 10.0, 1e-8),("3-A", 1234.0, 1e-3),("3-B", 1234.0, 1e-6),("3-C", 1234.0, 1e-8),
]# 方法列表
methods = [("线性逼近法", linear_approximation_sqrt),("二分查找法", binary_search_sqrt),("牛顿-拉夫逊法", newton_raphson_sqrt),
]# 主程序执行
with open("平方根实验结果.txt", "w", encoding="utf-8") as f, open("平方根实验结果.csv", "w", newline='', encoding="utf-8") as csv_f, open("records.csv", "w", newline='', encoding="utf-8") as rec_csv:writer = csv.writer(csv_f)rec_writer = csv.writer(rec_csv)writer.writerow(["实验ID", "方法", "近似值", "迭代次数", "运行时间(ns)", "记录点数", "初始误差", "最终误差", "平均误差变化"])rec_writer.writerow(["实验ID", "方法", "迭代次数", "近似值", "误差", "相对时间(ns)"])for exp_id, target, precision in experiments:print(f"Running experiment {exp_id} with target {target}, precision {precision}")true_sqrt = mpmath.sqrt(target)f.write(f"\n实验 {exp_id} (Target={target}, Precision={precision}):\n")for method_name, method_func in methods:print(f"  Running {method_name}...")approx, iter_count, time_ns, records = method_func(target, precision)print(f"  Finished {method_name}: approx {approx:.6f}, iterations {iter_count}, time {time_ns} ns")final_error = abs(float(approx) - float(true_sqrt))f.write(f"  {method_name}: 近似值={approx:.13f}, 迭代次数={iter_count}, 运行时间={time_ns} ns, 最终误差={final_error:.2e}\n")if records:for iter_num, app, rel_time in records:err = abs(float(app) - float(true_sqrt))rec_writer.writerow([exp_id, method_name, iter_num, f"{app:.13f}", f"{err:.2e}", rel_time])f.write(f"    记录已写入 records.csv\n")initial_error = abs(float(records[0][1]) - float(true_sqrt))f.write(f"    记录点数: {len(records)}, 初始误差: {initial_error:.2e}, 最终误差: {final_error:.2e}\n")if len(records) > 1:error_changes = [abs(float(records[i][1]) - float(true_sqrt)) - abs(float(records[i-1][1]) - float(true_sqrt)) for i in range(1, len(records))]avg_error_change = sum(error_changes) / len(error_changes)f.write(f"    平均误差变化: {avg_error_change:.2e}\n")else:avg_error_change = 0writer.writerow([exp_id, method_name, f"{approx:.13f}", iter_count, time_ns, len(records), f"{initial_error:.2e}", f"{final_error:.2e}", f"{avg_error_change:.2e}"])else:f.write("    无记录\n")writer.writerow([exp_id, method_name, f"{approx:.13f}", iter_count, time_ns, 0, "", "", ""])print("结果已写入平方根实验结果.txt、平方根实验结果.csv 和 records.csv")