数值问题

总述

涉及连续性数学问题

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

特点

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

    • 泰勒展开求解$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,则牛顿法失效
    • 导数绝对值越大,牛顿法越有效
  • 牛顿法不会把根括起来