Posted Updated Algorithm / Problem19 minutes read (About 2882 words)0 visits

数值问题

总述

涉及连续性数学问题

  • 解方程、方程组，计算定积分，函数求值
  • 和离散数学中：图、树、排序、组合相对

特点

  • 大部分此类问题只能近似求解

    • 泰勒展开求解$e^x$
    • Composite Trapezoidal rule：组合梯形法则，计算 定积分

  • 此类问题大部分要操作实数，而实数在计算机内部只能近似表示 ，大量对近近似数的算术操作可能会叠加误差，输出错误结果

整数乘法

俄式乘法

两个正整数n、m相乘的非主流算法

算法

  • 反复应用以下公式，简化每步的计算
  • 以$1 * m$作为算法终止条件

特点

  • n为奇数步骤中的m，可以最后累加即可

  • 算法中只有折半、加倍、相加操作

    • 手动计算非常简便
    • 计算机硬件对折半、加倍只需要移位就可

  • 减常因子法

大整数乘法

算法

考虑a、b两个n位整数，n为偶数

  • 从中间把数字分段，得到$a_1, a_0, b_1, b_0$

  • 则有

    • $c_2 = a_1 * b_1$
    • $c_0 = a_0 * b_0$
    • $c_1 = (a_1 + a_0) * (b_1 + b_0) - (c_2 + c_0)

  • 若n/2也是偶数，可以使用相同的方法计算得到三个乘法表达式

    • 若n为2的幂次，就可以得到计算n位数积的递归乘法
    • n迭代到足够小时，递归就可以停止

特点

  • 算法效率

    • 乘法次数递推式：$M(n)=3M(n/2), M(1)=1$，则 $M(n) = n^(log_2 3) \approx n^{1.585}$
    • 加法次数递推式：$A(n)=3A(n/2) + cn, A(1)=1$，则 $A(n) \in \Theta(n^{log_2 3})$

  • 算法有渐进效率优势，实际性能依赖于计算机系统、算法实现 质量，在某些情况下

    • 计算8位乘法时，分治算法速度快于传统方法
    • 计算超过100位时，速度是传统算法2倍

  • 分治法

欧几里得算法

计算最大公约数、最大公倍数

最大公约数

  • n为0，返回m作为结果结束
  • 将m处以n的余数赋给r
  • 将n付给m，r赋给n，返回第一步
1
2
3
4
5
6
Euclid(m, n)
	while n != 0 do
		r = m mod n
		m = n
		n = r
	return m

最大公倍数

  • 利用最大公约数计算最小公倍数

特点

  • 变治法（+减可变规模）

特定点求值

霍纳法则（计算多项式）

霍纳法则：不断将x作为公因子提取出来，合并降次后的项，然后 计算多项式在特定点的值

算法

1
2
3
4
5
6
7
8
Horner(P[0..n], x)
	// 用霍纳法则求多项式在给定点的值
	// 输入：多项式系数数组P[0..n]、数字x
	// 输出：多项式在x点的值
	p = P[n]
	for i = n-1 downto 0 do
		p = x*p + P[i]
	return p

特点

  • 算法效率

    • 效率始终为n，只相当于直接计算中$a_n x^n$的乘法数量

  • 变治法

二进制（计算）幂

将幂次转换为二进制位串，利用二进制位串简化计算

从左至右二进制幂

  • 对位串应用霍纳法则
1
2
3
4
5
6
7
8
9
10
LeftRightBinaryExponentiation(a, B[0..n-1])
	// 从左至右二进制幂算法计算a^n
	// 输入：数字a、表示幂次的二级制位串B[0..n-1]
	// 输出：a^n的值
	product = a
	for i = n-1 downto 0 do
		product = product * product
		if B[i] == 1:
			prduct = product * a
	return product

从右至左二进制幂

  • 累次计算二进制位串中为1部分值，将其累乘
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
RightLeftBinaryExponentiation(a, B[0..n-1])
	// 从右至左二进制幂算法
	// 输入：数字a、表示幂次的二级制位串B[0..n-1]
	// 输出：a^n的值
	term = a
	if B[0] == 1
		product = a
	else
		product = 1
		// 保存累乘值
	for i = i to n do
		term *= 2
		// 保存二进制位为1部分值
		if B[i] = 1
			product = product * term
	return product

特点

  • 算法效率

    • 两个算法效率取决于位串长度，是对数级的

  • 变治法

应用

  • 在密码技术中，需要对超过100位十进制整数进行乘法运算，而 计算机往往不能直接运算

矩阵乘法

Strassen矩阵乘法

算法

若A、B是两个n阶方阵（若n不是2幂次，考虑填充0）

  • 将A、B、C均分块为4个n/2子矩阵
  • 递归使用Strassen方程中定义的矩阵M进行计算计算C各个子阵

算法特点

  • 对2 * 2分块计算，Strassen算法执行了7次乘法、18次加减法， 蛮力算法需要执行8次乘法、4次加法

  • 算法效率

    • 乘法次数递推式：$M(n) = 7M(n/2), M(1) = 1$，则 $M(n) = 7^{log2 n} = n^{log_2 7} \approx n{2.807}$
    • 加法次数递推式：$A(n) = 7A(n/2) + 18(n/2)^2, A(1)=0$ ，则$A(n) \in \Theta(n^{log_2 7})$
    • 矩阵趋于无穷大时，算法表现出的渐进效率卓越

  • 还有一些算法能运行时间$\in \Theta(n^\alpha)$，最小能达到 2.376，但是这些算法乘法常量很大、算法复杂，没有实用价值

  • 矩阵乘法效率下界为$n^2$，目前得到的最优效率和其还有很大 距离

  • 分治法

线性方程组

  • 假设方程组系数矩阵为n阶方阵，且解唯一
  • 主要思想都是高斯消元法（变治法），只是出于效率、误差有 不同实现方式

前向消去法

算法

1
2
3
4
5
6
7
8
9
10
11
ForwardElimination(A[1..n, 1..n], b[1..n])
	// 对方程组扩展矩阵[A|b]使用高斯消元法
	// 输入：矩阵A[1..n, 1..n]，向量b[1..n]
	// 输出：扩展的上三角矩阵
	for i = 1 to n do
		A[i, n+1] = b[i]
		// 得到扩展矩阵
	for i = 1 to n-1 do
		for j = i+1 to n do
			for k = n+1 downto i do
				A[j, k] = A[j, k] - A[i, k]*A[j, i] / A[i, i]

算法特点

  • 前向消去法不一定正确

    • 如果A[i, i]==0，不能以其作为除数，此时需要交换行 （解唯一时总是存在非0行）
    • A[i, i]非常小，导致比例因子A[j, i] / A[i, i]非常大， 产生大的舍入误差

  • 最内层循环效率低

部分选主元法

算法

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
BetterForwardElimination(A[1..n, 1..n], b[1..n])
	// 用部分选主元法实现高斯消去
	// 输入：矩阵A[1..n, 1..n]，向量b[1..n]
	// 输出：扩展的上三角矩阵
	for i = 1 to n do
		A[i, n+1] = b[i]
	for i = 1 to n-1 do
		pivotrow = i
		for j = i+1 to n do
			if |A[j, i]| > A[pivot, i]
				pivotrow = j
				// 选择第i列系数最大的行作为第i次迭代基点
				// 保证比例因子绝对值不会大于1
			for k = i to n+1 do
				swap(A[i, k], A[pivot, k])
			for j = j+1 to n do
				temp = A[j, i] / A[i, i]
				// 这样只需要计算依次比例因子
				for k = i to n+1 do
					A[j, k] = A[j, k] - A[i, k] * temp

特点

  • 部分选主元法克服了前向消去法弊端
    • 最内层乘法（加法）执行次数为 $\frac {n(n-1)(2n+5) 6 \approx \frac n^3 3 \in \Theta(n^3)$
    • 始终能保证比例因子绝对值不大于1

反向替换法

在得到上三角系数矩阵中

  • 从最后一个方程中可以立刻求出$x_n$
  • 将$xn$带入倒数第二个方程求出$x{n-1}$
  • 逐次递推得到所以解

特点

  • 算法时间效率$\in \Theta(n^2)$

高斯消去法应用

  • 矩阵（可逆矩阵）中应用
    • LU分解（Doolittle分解）
    • Cholesky分解（要求矩阵正定）
    • 求逆
    • 求行列式
  • 高斯消元法整个算法效率取决于消去部分，是立方级
    • 事实上此方法在计算机上求解大规模方程组很难，因为舍入 误差在计算过程中会不断累积

非线性方程求解

  • 5次及以上多项式没有只包含多项式系数、算术操作、开根号 的通用求根公式

  • 方程的代数解并不具有很大的意义，充其量只是为方程的根设置 一个符号，然后再说方程有一个根等于这个符号（高斯）

平分法

基于连续函数界值定理，类似于连续版折半查找

算法

在区间[a, b]的端点上，$f(x)$符号取反

  • 计算$f(x{mid}), x{mid}= \frac {a+b} 2$的值
  • 若$f(x_{mid})=0$，则求得一个
  • 否则选择使得$f(x)$能在端点上取得相反值区间$[a, x{mid}$ 、$[x{mid}, b]$
  • 当包含根得区间小于预定义得$\epsilon > 0$时，就可以停止 算法
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
Bisection(f(x), a, b, eps, N)
	// 评分法求f(x) = 0的一个根
	// 输入：f(a)f(b) < 0，eps绝对误差上界，N迭代次数上界
	// 输出：(a, b)上的一个根近似（精确）值，或包含根的区间
	n = 1
		// 迭代计数
	while n <= N do
		x = (a + b)/2
		if x - a < eps
			 return x
		fval = f(x)
		if fval = 0
			return x
		if fval*f(a) < 0:
			b = x
		else
			a = x
		n += 1
	return a, b
		// 达到迭代限制，返回包含根的区间

特点

  • 求解精度
    • 理论上迭代次数足够$x_n$可以任意接近真实根$x^{*}$
    • 实际上，机器使用0表示非常小的值，$\epsilon$小于特定 机器阈值时，算法不会停止，也无法得到满足条件的解
    • d对目标函数求值时可能会发生舍入误差
  • 缺点
    • 相较于其他已经算法，收敛速度较慢
    • 并且无法扩展到更加一般的方程、方程组领域
  • 优点
    • 区间特性容易检验

Method of False Position

试位法：类似于连续版差值查找

算法

  • 类似平分法每次也使用某个区间$[a_n, b_n]$括住连续函数的 根，函数在端点取值符号相反

  • 使用穿过$(a_n, f(a_n)), (b_n, f(b_n))$的直线在x轴截距 $x=\frac {a_nf(b_n) - b_nf(a_n)} {f(b_n) - f(a_n)}$作为 分割点

特点

  • 对许多实例，试位法收敛速度较平分法更快

牛顿法

Newton-Raphon Method

算法

  • 方法产生近似解序列：函数切线在x轴截距 $x_{n+1} = x_n - \frac {f(x_n)} {f^{‘}(x_n)},n=0,1,\cdots$

特点

  • 大多数情况下，若初值$x_0$足够接近根，牛顿法能够保证序列 收敛与根；对远离根的初值，无法保证一定会收敛

  • 优点

    • 牛顿法相较于平分法、试位法收敛速度更快，选定合适的 初值能够快速收敛
    • 能够应用于更一般类型的方程、方程组

  • 方法每次都迭代需要重新求函数、导数值

    • 导数值等于0，则牛顿法失效
    • 导数绝对值越大，牛顿法越有效

  • 牛顿法不会把根括起来

数值问题

https://xyy15926.github.io/Algorithm/Problem/numeric_analysis.html

Author

UBeaRLy

Posted on

2019-03-21

Updated on

2021-08-02

Licensed under

#

Like this article? Support the author with

Comments

Categories

Follow

Recents

Advertisement

follow.it

×