1 计算方法概论

Posted 2021-03-16 philolif

1 概论

计算方法又称数值分析，主要研究科研和工程中数学问题的算法设计及相关理论。这一部分主要介绍计算方法当中一些基本的概念。

1.1 引言

1.1.1 计算方法的意义

由于计算机性能的快速发展，极大拓展了数学的应用范围和能力。现今人们可以通过高性能计算机来模拟各种自然现象的产生（如计算电磁学，计算化学），这些模拟实验离不开高性能的硬件和高效的计算方法。使用时间、空间复杂度低的算法，可以在同等硬件条件下，更加高效地完成任务。

1.1.2 研究对象和特点

数值计算方法：研究适合利用计算机进行科学计算的方法，计算机的特点是只能计算离散，有限的问题。

基本步骤是：

graph LR 实际问题-->建立数学模型; 建立数学模型-->算法设计; 算法设计-->程序设计; 程序设计-->上机计算; 上机计算-->求得结果;

其中计算方法重点研究算法设计。

1.2 算法与效率

1.2.1 算法

解决某类数学问题的数值方法成为数值算法，为了能够使算法在计算机上实现，必须将待解决问题分解成有限次的四则运算问题。还需要注意的是不同的算法常常具有不同的时间、空间复杂度，算法选择得不好有可能导致计算结果迟迟无法获得。

算法常具有的特征：

1. 表现为一个无穷过程的截断

   求 (sin(x),xin[0,{pi}/{4}]) 的值，其中 (sin(x)) 的泰勒展开：

   [sin(x)=x-frac{x^3}{3!}+frac{x^5}{5!}-frac{x^7}{7!}+cdots+(-1)^nfrac{x^{2n+1}}{(2n+1)!}+cdots ]

   为了能够使用计算机来计算，并且精度又能满足要求，我们只能取该无穷级数的前几项进行计算，作为 (sin(x)) 的近似值。

   例如这里我们选择前n项进行近似计算，根据拉格朗日余项便可以估计误差的大小：

   [|R_n|=|(-1)^{2n+3}frac{xi^{2n+3}}{(2n+3)!}|lefrac{1}{(2n+3)!}(frac{pi}{4})^{2n+3},xiin[0,pi/4] ]

   通过无穷级数进行计算的时候，需要注意级数的收敛性和误差。

2. 表现为连续过程的离散化

   计算积分：

   [I=int_0^1frac{1}{1+x}dx ]

   根据黎曼积分的定义进行计算，将原始区间均分为n份，即：

   [[0,1]Rightarrow[0,frac{1}n],cdots,[frac{n-1}n,1] ]

   并且令 (x_iin[(i-1)/n,i/n]) ，那么

   [Iapproxfrac{1}nsum_{i=1}^n{f(x_i)} ]

   当n越大，那么精度越高。

3. 表现为迭代的形式

   迭代是指简单算法的多次重复，后一次使用前一次的结果，在程序中表现为循环。一个简单的例子通过迭代求 (sqrt{a}) 的近似值：

   [x_{k+1}=frac{1}2(x_k+frac{a}{x_k}) ]

   可以证明上述数列收敛于 (sqrt{a}) ，因此可以给定一个初始值 (x_0) 便可以通过迭代求解。上述式子的一个直观理解是：因为 ((x_k)(a/x_k)=a) ，所以 (x_k) 和 (a/x_k) 必然有一个大于等于 (sqrt{a}) ，另一个小于等于 (sqrt{a}) 。因此两者的均值将会是 (sqrt{a}) 的一个更好近似值。

4. 采用构造性证明方法

   构造性证明是指若要证明一个问题有解，就具体地构造出一种求解该问题的具体步骤和过程，这样一来不仅证明了原问题有解，并给出了一种求解的方法。而非构造性证明则只是证明解的存在性，并未给出具体的算法。

1.2.2 算法的效率

同一个数学问题可能会有多种数值算法，如何比较一个算法的好坏，主要通过两个指标：计算结果和计算精度。算法的计算效率是由它的计算工作量体现的，由于在计算机中进行乘除运算（浮点数运算）的开销比加减运算的开销大得多，因此我们统计算法的计算量时可以仅关注算法进行浮点数运算的次数（一次浮点数运算记为1flop）。我们很自然地会选择精度满足要求，并且计算量尽可能少的算法。

这里将讨论计算多项式的两种不同解法：

[p(x)=a_nx^n+a_{n-1}x^{n-1}+cdots+a_1x+a_0 ]

算法1：

初始化

[left{egin{array}{} t_0=1u_0=a_0 end{array} ight. ]

迭代过程：

[left{egin{array}{} t_k=xt_{k-1}u_k=u_{k-1}+a_kt_k end{array} ight.quad k=1,2,cdots,n ]

迭代n次即可求得所需的结果，计算量为2n次乘法，以及n次加法，因此计算量为 2nflop 。

算法2（秦九韶算法）：

[egin{split} p(x)&=a_nx^n+a_{n-1}x^{n-1}+cdots+a_1x+a_0&=(a_nx^{n-1}+a_{n-1}x^{n-2}+cdots+a_1)x+a_0&=(cdots(a_nx+a_{n-1})x+cdots+a_1)x+a_0 end{split} ]

迭代过程：

[left{egin{array}{} v_0=a_nv_k=v_{k-1}x+a_{n-k} end{array} ight.quad k=1,2,cdots,n ]

此时就仅需n次加法和n次乘法即可，因此计算量为 (n flop)

1.3 误差

进行测量的时候不可避免地会出现误差，如何将误差降低到可以接受的范围，便是我们需要研究的。

1.3.1 误差的来源

1. 模型误差：由于实际问题往往都是高度复杂的，在建立数学模型的时候，会不可避免地忽略某些次要的因素，进行必要的简化，这就带来了与实际问题的误差。
2. 测量误差：测量已知参数时，数据带来的误差。
3. 截断误差：设计算法时，进行必要的近似处理，例如收敛的无穷级数计算只取前几项。
4. 舍入误差：由于计算机所能表示的数范围是有限的（模运算），因此每一步都可能出现四舍五入的现象，这时就出现了舍入误差。

数值分析主要考虑截断误差，以及舍入误差。其余因素带来的误差一般不考虑。

1.3.2 误差的概念

误差的定义： 设 (x) 为准确值，(p) 是 (x) 的一个近似值，则 (e(p)=x-p) 成为近似值 (p) 的绝对误差，简称误差。当 (x eq0) 时，将 (e_r(p)=e(p)/x) 成为近似值 (p) 的相对误差。

因为准确值 (x) 往往是未知的（不然也就没有必要使用近似值了），因此绝对误差和相对误差都是未知的，无法知道确切值。为此我们引入绝对误差限（误差限）和相对误差限。

误差限： 如果近似值 (p) 满足 (|e(p)|=|x-p|ledelta{p}) ，则称 (delta p) 为近似值 (p) 的绝对误差限。将 (delta_rp=delta{p}/|x|) 称为相对误差限，但 (x) 未知，一般使用 (delta_rp=delta{p}/|p|) 作为相对误差限。

需要注意的是，绝对误差是有单位的，而相对误差的量纲为单位1。在多数情况下，相对误差更加有用。

接下来我们通过计算四舍五入方法的绝对误差限，引入有效数字的概念，假设准确值为 (x) ，近似值为 (p)：

[x=0.p_1p_2cdots p_ncdots imes10^mp=left{egin{split} &pm0.p_1p_2cdots p_n imes10^m,&p_{n+1}le4 ext{(四舍)}&pm0.p_1p_2cdots(p_n+1) imes10^m,;;&p_{n+1}ge5 ext{(五入)} end{split} ight. ]

四舍时的误差限：

[egin{split}|x-p| &=|0.p_1p_2cdots p_ncdots-p| imes10^m&le|0.p_1p_2cdots p_n499cdots-0.p_1p_2cdots p_n| imes10^m&=0.underbrace{00cdots0}_{n}499cdots imes10^m&=0.5 imes10^{m-n} end{split} ]

这里 ( 0.499cdots=0.5) 。

五入时的误差限：

[egin{split}|x-p| &=|0.p_1p_2cdots(p_n+1)-0.p_1p_2cdots p_ncdots| imes10^m&=|0.underbrace{0cdots0}_{n-1}1-0.0cdots p_{n+1}p_{n+2}cdots| imes10^m&le{(1-0.p_{n+1})} imes10^{m-n}&le0.5 imes10^{m-n} end{split} ]

综上可知四舍五入的绝对误差限是 ( 0.5 imes{10^{m-n}}) ，也即误差限是被保留最后数位上的半个单位。如果一个数是经过四舍五入得到的，那么我们将绝对误差限中的n定义为该近似值的有效数字。 注意这里的n值应该是将近似值改写成 ( 0.p_1p_2cdots) 这种形式，并且 (p_1 e0) 时候的值。

还有几点需要注意的是，小数点的左右移动并不影响有效数字，并且如果一个数是经过四舍五入得到的近似值，那么该数的每一位都是有效数字。

更严格来说，如果判断一个近似值有多少个有效数字应该按照上述定义进行，也就是说找到一个最大的整数n使得下式成立：

[|x-p|le0.5 imes10^{m-n} ]

其中 (x) 是准确值， (p) 是近似值。

如果 (x) 和 (p) 满足以下条件：

[|x-p|le0.5 imes10^{-n} ]

其中 (n) 是大于零的正整数，此时我们称近似值 (p) 精确到小数点后第 (n) 位。

有效数字和相对误差具有以下定理：假设 (p=pm0.p_1p_2cdots {p_n} imes10^m,p_1 e0) 是 (x) 的近似值

1. (p) 具有n位有效数字，则相对误差满足

   [left|e_r(p) ight|=left|frac{e(p)}{p} ight|=left|frac{x-p}{p} ight| leleft|frac{0.5 imes10^{m-n}}{0.p_1 imes10^{m}} ight|=frac{5}{p_1} imes10^{-n} ]

2. 如果相对误差满足以下不等式，则近似值 (p) 至少有 (n) 位有效数字

   [left|e_r(p) ight|lefrac{5}{p_1+1} imes10^{-n} ]

   证明如下：

   [egin{split} |x-p|=|p|left|frac{x-p}{p} ight|=|p||e_r(p)|&le|p|frac{5}{p_1+1} imes10^{-n}&le|0.(p_1+1) imes10^m|frac{5}{p_1+1} imes10^{-n}\&=0.5 imes10^{m-n} end{split} ]

   根据有效数字的定义可知，近似值 (p) 至少有n位有效数字

1.3.3 误差的传播

假设函数 (f(x)) 二阶可微， (p) 是 (x) 的近似值，利用带拉格朗日余项的泰勒展开：

[f(x)-f(p)=f^prime(p)(x-p)+frac{f^{primeprime}(xi)}{2}(x-p)^2 ]

其中 (xi) 位于 (x) 和 (p) 两个数之间。

忽略二次项，并取绝对值可以得到：

[|f(x)-f(p)|le|f^prime(p)|delta{p} ]

也就是：

[delta{f(p)}le|f^prime(p)|delta{p}qquad ext{or}qquaddelta_rf(p)lefrac{|f^prime(p)|delta{p}}{|f(p)|} ]

事实上只有当近似值足够靠近真实值时，忽略高阶项的操作才是可取的。

同样对于多元函数，利用多元函数的泰勒展开得到一阶近似，或者说一阶微分可以得到误差估计式：

[e(u)=|f(x,y)-f(p,q)|leleft|frac{partial{f(p,q)}}{partial x} ight|delta{p}+ left|frac{partial{f(p,q)}}{partial y} ight|delta{q} ]

也即是说：

[delta{u}leleft|frac{partial{f(p,q)}}{partial x} ight|delta{p}+ left|frac{partial{f(p,q)}}{partial y} ight|delta{q}\delta_r{u}lefrac{left|frac{partial{f(p,q)}}{partial x} ight|delta{p}+ left|frac{partial{f(p,q)}}{partial y} ight|delta{q}}{|f(p,q)|} ]

利用上述公式，当给定具体的函数之后，便可以估计误差：

[delta{(ppm q)}ledelta(p)+delta{(q)}\delta_r{(ppm q)}lefrac{delta(p)+delta{(q)}}{|ppm q|}\delta(pq)le{|q|delta{(p)}}+|p|delta(q) ]

1.4 数值计算中的几个原则

1.避免两个相近数相减 ：两个相近的数相减，有效数字会大大损失

[delta_r(x-y)=frac{delta(x-y)}{|x-y|}lefrac{delta(x)+delta(y)}{|x-y|} ]

显然，当x和y两数相近时， (|x-y|) 会很小，导致相对误差限很大。

例如：

[sqrt{170}-13=0.0384048cdots ]

如果选择4位有效数字进行计算，那么结果只有一位有效数字：

[sqrt{170}-13approx13.04-13=0.04 ]

将其修改为：

[sqrt{170}-13=frac{1}{sqrt{170}+13}=frac{1}{13.04+13}approx0.03840 ]

就可以得到4位有效数字的结果。

2.避免除数的绝对值远小于被除数的绝对值

[delta(frac{x}{y})lefrac{|x|delta(y)+|y|delta(x)}{y^2} ]

如果除数 (|y|) 相对 (|x|) 来说，很小的话会导致误差限变得很大。

3.防止大数吃小数

计算机在进行运算的时候，首先要将参加运算的数都表示成绝对值小于1，而阶码相同的数。

例如，为了计算:( 10^{9}+1) ，就需要将其改写为:( 0.1 imes10^{10}+0.0000000001 imes10^{10})。

如果计算机只能表示8位小数，那么最终的计算结果就是 ( 0.1 imes10^{10}) 。这相当于是大数“吃掉了”小数，这种情况有时允许，而有时候应该避免。

4.简化计算步骤，减少运算次数，避免误差积累

这一点也很好理解，在计算机中浮点数本身所能表示的数就是有限的（因此计算过程中不可避免地出现一定的误差），所以计算过程中使用步骤少，计算简单的方法一般都能够有效地减小最终结果的误差。

1.5 病态问题、数值稳定性与避免误差危害

1.5.1 病态问题与条件数

计算函数值 (f(x)) 时使用近似值 (xprime) 代替 (x) ，则其相对误差： ((x-x^prime)/x) ，函数值的相对误差： ((f(x)-f(x^prime))/f(x)) ，定义条件数 (C) ：

[left|frac{(f(x)-f(x^prime))/f(x)}{(x-x^prime)/x} ight|approx left|frac{xf^prime(x)}{f(x)} ight|=C ]

如果条件数 (C) 的值很大，就说明当 (x) 发生一个微小的相对变化时，也会使得函数值发生很大的变化。这也说明即使我们选择的近似值 (x‘) 和实际值 (x) 很接近，也会使得最终的计算结果 (f(x‘)) 与实际值 (f(x)) 差距很大。出现这种现象的问题称为病态问题，非病态问题则称为良态问题。病态和良态是相对的，没有严格的界定。一般认为当 (Cgg1) 时，就是一个病态问题。

1.5.2 数值方法的稳定性

一个算法，如果初始数据的微小变化，引起最终结果的改变也是微小的，那么该算法就是数值稳定的算法（否则称为数值不稳定），这也意味着该问题是非病态问题。如果对于某些特定的数据算法是稳定的，而对于其他数据是不稳定的，则称算法具有条件稳定性。若对任何数据，算法都是稳定的，则称算法具有无条件稳定性。

假设一个算法有初始误差 (varepsilon_0>0) ，由它引起的算法执行第n步之后的误差是 (varepsilon_n) ，有两种情况会发生：

1. 如果 (varepsilon_napprox Cnvarepsilon_0) ，其中 (C) 是与 (n) 无关的常数，则称误差的增长是线性的。
2. 如果 (varepsilon_napprox C^nvarepsilon_0) ，则称误差的增长是指数型的。

线性增长的误差往往不可避免，如果 (C) 和 (varepsilon_0) 均较小，结果一般可接受。对于指数型的误差增长， (C>1) 的情况应该尽量避免，否则当计算的步骤n较大时，会导致误差非常大。