Python常用算法——迭代算法
深入浅出迭代算法：从数学思想到Python实战
引言
在计算机科学和数值分析领域，迭代算法是一种最基本、最重要的算法思想。无论是求解复杂的数学方程，还是处理大规模线性方程组，迭代法都扮演着不可替代的角色。本文将带领读者深入理解迭代算法的核心思想，并通过多个Python实战案例，展示迭代算法在不同场景下的应用。
一、迭代算法思想基础
1.1 什么是迭代算法
迭代法（Iterative Method），也称辗转法，是一种不断用变量的旧值递推新值的过程。简单来说，就是让计算机重复执行一组指令，每次执行都从变量的原值推出一个新值。
与迭代法相对应的是直接法（一次解法），例如通过开方直接求解方程 x³ = 4。直接法能够一次性解决问题，但在面对复杂方程（如五次及以上代数方程）或大规模线性方程组时，往往难以找到直接解法，此时迭代法就成为重要的替代方案。
1.2 迭代法的数学定义
从数学角度看，对于给定的线性方程组，通常可以变换成如下形式：
x = Bx + f
其中x、B、f同为矩阵。通过如下公式逐步代入求近似解的方法称为迭代法：
xₖ₊₁ = Bxₖ + f
这里xₖ代表迭代k次得到的x。如果当k→∞时，xₖ趋近于某个x*，则称此法收敛，x*就是方程组的解；否则称为迭代发散。
1.3 迭代算法的三个关键要素
要成功运用迭代算法解决问题，需要把握好以下三个方面：
（1）确定迭代变量 在要解决的问题中，至少存在一个可以直接或间接地从旧值递推出新值的变量，这个变量就是迭代变量。
（2）建立迭代关系表达式 这是解决迭代问题的关键。迭代关系式描述了如何从前一个值推出下一个值的公式或关系。在实际应用中，通常通过递推或倒推等方法建立迭代关系式。
（3）控制迭代过程 必须确定迭代何时结束，不能让程序无限循环下去。通常有两种情况：
- 迭代次数确定：构建固定次数的循环来控制
- 迭代次数不确定：需要分析出结束迭代的条件（如误差小于某个阈值）
1.4 递归与迭代的区别
很多初学者容易混淆递归和迭代。它们的主要区别在于：
递归：自顶向下逐步拓展需求，最后自底向上逐项运算。即从f(n)拓展到f(1)，再由f(1)逐步算回f(n)。递归是在函数内部调用自身。
迭代：直接自底向上逐项运算，由f(1)算到f(n)。迭代是循环求值的过程。
虽然递归的效率通常低于迭代，但随着计算机性能的提升，且递归代码更易于理解、可读性强，对于初学者来说，递归往往是更好的入门选择。
二、经典迭代算法实战案例
下面通过多个Python实例，展示迭代算法在不同场景下的应用。
2.1 斐波那契数列（fib.py）
斐波那契数列是最经典的迭代问题之一：每个数等于前两个数之和。下面是用迭代算法实现的斐波那契数列计算：
def fab(n):
n1 = 1
n2 = 1
n3 = 1
if n < 1:
print('对不起，输入有误！')
return -1
else:
while (n - 2) > 0:
n3 = n2 + n1
n1 = n2
n2 = n3
n -= 1
return n3
month = int(input('请输入月数：'))
result = fab(month)
print("%d 个月后的兔子数量为 %d" % (month, result))
这个程序通过while循环不断迭代更新n1、n2、n3的值，直到计算出第n个月的兔子数量。每次迭代中，n3存储最新的计算结果，n1和n2分别向前移动一位。
2.2 牛顿迭代法求方程根（niudun01.py）
牛顿迭代法（Newton's Method）是求解方程根的经典方法。下面是用牛顿法求解方程 f(x)=x³-x-1 的根的代码：
def f(xi):
return xi*xi*xi - xi - 1
def f1(xi):
return 3*xi*xi - 1
x = []
x.append(0.5)
eps = 1e-14
error = abs(f(x[-1]))
number_iteration = 0
while error > eps:
x.append(x[-1] - f(x[-1]) / f1(x[-1]))
error = abs(f(x[-1]))
number_iteration += 1
print('牛顿法迭代次数为%d次' % number_iteration)
print('方程的根x*为%f' % x[-1])
print('f(x*)的值为%f' % f(x[-1]))
牛顿迭代法的核心思想是：在当前位置用切线近似替代原函数，将切线与x轴的交点作为下一个近似点。迭代公式为：xₖ₊₁ = xₖ - f(xₖ)/f'(xₖ)。当|f(xₖ)|小于设定的误差限时，迭代结束。
2.3 通用的牛顿迭代框架（niudun05.py）
下面的代码实现了一个更通用的牛顿迭代法，用户可以输入任意表达式：
from sympy import *
import math
import sys
x = symbols("x")
f = input("请输入方程的表达式：")
f = integrate(f, x)
f = diff(f)
list_in = input("请输入初始值x,极小值sigma以及想要计算的次数: ").split(" ")
x0 = float(list_in[0])
sigma = float(list_in[1])
N = int(list_in[2])
if diff(f).subs(x, x0) == 0:
print(x0)
sys.exit()
account = 1
x_new = x0
while account <= N:
x_new = x0 - f.subs(x, x0) / diff(f).subs(x, x0)
if math.fabs(x_new - x0) < sigma:
print(x_new)
sys.exit()
x0 = x_new
account += 1
print("无法解出，不满足条件")
这个程序利用了SymPy库进行符号运算，用户可以输入任意数学表达式，程序会先对表达式进行积分再求导，然后应用牛顿迭代法求解。这种灵活性使得代码可以适应各种不同的方程。
2.4 牛顿法求平方根（niudun03.py）
牛顿迭代法也可以用来求平方根。下面是比较二分法和牛顿法求平方根的完整实现：
import math
from math import sqrt
def check_precision(l, h, p, len1):
"""检查是否达到了精确位"""
l = str(l)
h = str(h)
if len(l) <= len1 + p or len(h) <= len1 + p:
return False
for i in range(len1, p + len1):
if l[i] != h[i]:
return False
return True
def print_result(x, len1, p):
"""按精度要求打印结果"""
x = str(x)
if len(x) - len1 < p:
s = x[:len1] + "." + x[len1:] + '0' * (p - len(x) + len1)
else:
s = x[:len1] + "." + x[len1:len1 + p]
print(s)
def newton_sqrt(x, p):
"""牛顿迭代法求平方根"""
x0 = int(sqrt(x))
if x0 * x0 == x:
print_result(x0, len(str(x0)), p)
return
len1 = len(str(x0))
g = 1
q = x // g
g = (g + q) // 2
while not check_precision(g, q, p, len1):
x = x * 100
g = g * 10
q = x // g
g = (g + q) // 2
return print_result(g, len1, p)
# 主程序循环
while True:
x = int(input("请输入待开方数："))
p = int(input("请输入精度："))
print("newton_sqrt:", end="")
newton_sqrt(x, p)
这个程序的核心思想是：求x的平方根等价于求方程f(g)=g²-x=0的根。牛顿迭代公式为：gₖ₊₁ = (gₖ + x/gₖ)/2。每次迭代都取猜测值和商的平均值，这样可以快速逼近真实的平方根。
2.5 牛顿法求解多元函数极值（niudun04.py）
牛顿法不仅可以求解方程，还可以用于多元函数的优化问题。下面是用牛顿法求解Rosenbrock函数的极小值点：
import numpy as np
import matplotlib.pyplot as plt
def jacobian(x):
"""梯度（一阶导数）"""
return np.array([-400*x[0]*(x[1]-x[0]**2)-2*(1-x[0]),
200*(x[1]-x[0]**2)])
def hessian(x):
"""海森矩阵（二阶导数）"""
return np.array([[-400*(x[1]-3*x[0]**2)+2, -400*x[0]],
[-400*x[0], 200]])
def newton(x0):
"""牛顿法优化"""
W = np.zeros((2, 10**3))
i = 1
imax = 1000
W[:, 0] = x0
x = x0
delta = 1
alpha = 1
while i < imax and delta > 10**(-5):
p = -np.dot(np.linalg.inv(hessian(x)), jacobian(x))
x0 = x
x = x + alpha * p
W[:, i] = x
delta = sum((x - x0)**2)
print('第', i, '次迭代结果:', x)
i += 1
return W[:, 0:i]
# 可视化
X1 = np.arange(-1.5, 1.5+0.05, 0.05)
X2 = np.arange(-3.5, 2+0.05, 0.05)
[x1, x2] = np.meshgrid(X1, X2)
f = 100*(x2-x1**2)**2 + (1-x1)**2
plt.contour(x1, x2, f, 20)
x0 = np.array([-1.2, 1])
W = newton(x0)
plt.plot(W[0, :], W[1, :], 'g*', W[0, :], W[1, :])
plt.show()
Rosenbrock函数是一个经典的非凸测试函数，其全局最小值位于(1,1)处。牛顿法利用目标函数的二阶导数（海森矩阵）信息，能够快速收敛到极值点。程序最后还通过等高线图直观展示了迭代点的收敛轨迹。
2.6 雅可比迭代法（yakebi.py）
雅可比迭代法（Jacobi Iteration）是求解线性方程组的基本迭代方法。下面是用Python实现的雅可比迭代法：
import numpy as np
max_iter = 100
Delta = 0.01
n = 3 # 矩阵维数
# 读取数据
f = open("123.txt")
# 文件内容：
# 2,-1,-1
# 1,5,-1
# 1,1,10
# -5,11,8
def read_tensor(f, shape):
data = []
for i in range(n):
line = f.readline()
data.append(list(map(eval, line.split(","))))
return np.array(data).reshape(shape)
def read_vector(f):
line = f.readline()
line = line.replace("\n", "")
line = list(map(eval, line.split(",")))
return np.array(line)
shape = (n, n)
A = read_tensor(f, shape)
b = read_vector(f)
f.close()
print('A:', A)
print('b:', b)
# 求LU = L + U（非对角线元素）
LU = np.copy(A)
for i in range(n):
LU[i, i] = 0
LU = 0 - LU
# 求D（对角线矩阵）
D = A - LU
# 迭代求解
x = np.ones(n)
DLU = np.dot(np.linalg.inv(D), LU)
Db = np.dot(np.linalg.inv(D), b)
print('初始x:', x)
for iteration in range(max_iter):
y = np.dot(DLU, x) + Db
if np.max(np.fabs(x - y)) < Delta:
print('迭代次数:', iteration)
break
x = np.copy(y)
print('x:', x)
雅可比迭代法的基本思想是：将方程组Ax=b改写为x=D⁻¹(b-(L+U)x)，然后进行迭代求解。每次迭代都使用上一轮的所有分量值来计算新一轮的各个分量。
2.7 高斯-赛德尔迭代法（Gauss-Seidel.py）
高斯-赛德尔迭代法（Gauss-Seidel Method）是雅可比迭代法的改进版本：
import numpy as np
max_iter = 100
Delta = 0.01
n = 3
# 读取数据
f = open("123.txt")
# 文件内容同上
# ...
# 求U（上三角矩阵）
U = np.copy(A)
for i in range(n):
for j in range(n):
if j <= i:
U[i, j] = 0
U = 0 - U
# 求DL（下三角矩阵）
DL = A + U
# 迭代求解
x = np.ones(n)
DLU = np.dot(np.linalg.inv(DL), U)
DLb = np.dot(np.linalg.inv(DL), b)
print('初始x:', x)
for iteration in range(max_iter):
y = np.dot(DLU, x) + DLb
if np.max(np.fabs(x - y)) < Delta:
print('迭代次数:', iteration)
break
x = np.copy(y)
print('x:', x)
与雅可比迭代法不同，高斯-赛德尔迭代法在计算第i个分量时，会立即使用已经更新的前i-1个分量值，因此通常收敛速度更快。
2.8 更紧凑的高斯-赛德尔实现（Gauss-Seidel03.py）
下面是一个更紧凑的实现方式，直接在迭代过程中更新变量：
import numpy as np
A = np.array([[5.0, 2, 1], [-1, 4, 2], [2, -3, 10]])
B = np.array([-12.0, 20, 3])
x0 = np.array([1.0, 1, 1])
x = np.array([0.0, 0, 0])
times = 0
while True:
tempx = x0.copy()
for i in range(3):
temp = 0
for j in range(3):
if i != j:
temp += x0[j] * A[i][j]
x[i] = (B[i] - temp) / A[i][i]
x0[i] = x[i].copy()
calTemp = max(abs(x - tempx))
times += 1
if calTemp < 1e-4:
break
print('迭代次数:', times)
print('解:', x)
这种实现方式更加直观，在每次内循环中计算完x[i]后立即更新x0[i]，用于后续分量的计算。
2.9 简单迭代法求解非线性方程组（Iteration.py）
迭代法也可以用于求解非线性方程组。下面是用简单迭代法求解一个三维非线性方程组的例子：
import math
def takeStep(Xcur):
"""迭代函数：从当前值计算下一个近似值"""
Xnex = [0, 0, 0]
Xnex[0] = math.cos(Xcur[1] * Xcur[2]) / 3.0 + 1.0 / 6
Xnex[1] = math.sqrt(Xcur[0] * Xcur[0] + math.sin(Xcur[2]) + 1.06) / 9.0 - 0.1
Xnex[2] = -1 * math.exp(-1 * Xcur[0] * Xcur[1]) / 20.0 - (10 * math.pi - 3) / 60
return Xnex
def ColculateDistance(Xcur, Xnew):
"""计算两次迭代之间的距离"""
temp = [Xcur[0]-Xnew[0], Xcur[1]-Xnew[1], Xcur[2]-Xnew[2]]
dis = math.sqrt(temp[0]*temp[0] + temp[1]*temp[1] + temp[2]*temp[2])
return dis
def iteration(eps, maxIter):
"""主迭代函数"""
cur_eps = 10000
Xcur = [0.1, 0.1, -0.1]
Xnew = [0, 0, 0]
iterNum = 1
print("--------------------------开始迭代--------------------------")
print(" 迭代次数 | Xk1 | Xk2 | Xk3 | eps ")
while cur_eps > eps and iterNum < maxIter:
Xnew = takeStep(Xcur)
cur_eps = ColculateDistance(Xcur, Xnew)
print(" %d %.8f %.8f %.8f %.8f" %
(iterNum, Xcur[0], Xcur[1], Xcur[2], cur_eps))
iterNum += 1
Xcur = Xnew
return 0
iteration(10**-10, 200)
这个例子展示了如何将非线性方程组改写成x = g(x)的形式，然后通过不断迭代逼近方程组的解。程序会输出每次迭代的中间结果和误差，直观展示收敛过程。
2.10 角谷猜想（jiaogu.py）
最后，我们来看一个有趣的数学问题——角谷猜想（也称冰雹猜想）：
n = int(input("请输入一个正整数："))
while n != 1:
if n % 2 == 0:
k = n / 2
print("%d/2=%d" % (n, k))
n = k
else:
l = n * 3 + 1
print("%d*3+1=%d" % (n, l))
n = l
这个猜想的内容是：对于任意正整数，如果是偶数则除以2，如果是奇数则乘以3再加1，最终总会进入4→2→1的循环。虽然这个猜想至今未被证明，但我们可以用迭代程序验证对于给定的输入，这个迭代过程确实会收敛到1。
三、总结与思考
通过上述多个实例，我们可以看到迭代算法的广泛应用和强大功能。从简单的斐波那契数列，到复杂的牛顿法求根，再到线性方程组的雅可比和高斯-赛德尔迭代，以及非线性方程组的求解，迭代算法贯穿了数值计算的方方面面。
3.1 迭代算法的优点
- 简单易懂：迭代算法的逻辑通常比较直观，容易理解和实现
- 适用范围广：从简单数列到复杂方程组，迭代法都能派上用场
- 可控性好：可以通过调整误差限来控制计算精度
- 适合计算机：充分发挥了计算机擅长重复运算的特点
3.2 使用迭代算法的注意事项
- 收敛性：不是所有迭代都收敛，需要分析迭代格式的收敛条件
- 初始值选择：很多迭代方法（如牛顿法）对初始值敏感，需要合理选择
- 终止条件：需要合理设置迭代终止的条件，避免无限循环或过早停止
- 计算效率：对于大规模问题，需要考虑迭代的效率，选择合适的迭代方法
3.3 结语
迭代算法作为计算机科学和数值分析的基础思想，值得我们深入学习和掌握。希望通过本文的讲解和代码示例，能够帮助读者更好地理解迭代算法的本质，并在实际问题中灵活运用。
正如文章开头所说，迭代法的核心是"用旧值递推新值"。这种思想不仅体现在编程中，在我们的日常生活和学习中也处处可见。每一次的尝试和修正，都是在向目标不断靠近——这或许就是迭代算法带给我们的人生启示。
