SCAU华南农业大学数值计算期末复习笔记
|
457
|
0
|
SCAU

6745 字
|
33 分钟

数值计算问题概述

算法的描述

我们可用用下面一些方法来描述算法：

  1. 自然语言：重要会议的议程，一台晚会的节目单都可以看成使用自然语言描述的算法。用自然语言描述算法可用第一步，第二步，…，第n步的格式书写，形式灵活，容易理解，但不便于编程。
  2. 流程图：用带箭号的线表示执行流程，用不同的方框配合以文字表示不同的操作，最大的特点是结构清晰，含义准确，便于指导编程。算法语言：语法规定过于严格，不利于教学。
  3. 伪语言：用c语言的关键词来说明流程，用数学形式来说明运算，既容易编程，又容易理解。

关于算法的评价

  1. 程序的可读性好
  2. 节省计算时间
  3. 节省存储空间
  4. 数字稳定性好

计算量

  1. 概念：我们把作一次浮点数的乘法运算连同一次加法运算的计算量作为计算量的计量单位。提示：由于计算机执行加法运算的时间相对说来要少得多，所以我们也可以把作一次浮点数的乘法运算的计算量作为计算量的计量单位。
  2. 计算方式
    假设A=(ai j)l×m，B=(bi j)m×n分别为  k×m 和  m×n矩阵，记C=A·B,则 C 为k× n 矩阵，由于
    ci j=ai 1·b1 j+ ai 2·b2 j +...+ai m·bm j
    矩阵乘法的计算量可以通过以下公式计算： 

    k * m * n

其中，k 表示矩阵 A 的行数，m 表示矩阵 A 的列数和矩阵 B 的行数，n 表示矩阵 B 的列数。

例题

误差分析

误差来源

误差的来源是多方面的，但主要来源为：过失误差，描述误差，观测误差，截断误差和舍入误差。

  1. 过失误差
    • 过失误差通常是指由于设备的故障和人为的错误所导致的误差，比如物理设备在恶劣的工作环境下可能会是某些性能指标变坏，原始数据纪录错误，原始数据输入错误等。
    • 现在，由于计算机的普及，人人都可以自己变程序利用机器计算，所以产生过失误差的机会实际上大大增多了，所以“如何有效避免过失误差”的问题实际上也变得更加重要了。
    • 比如求一元二次方程$ ax2+bx+c=0
      x, =, \frac {-b\pm \sqrt {{b}^{2}-4ac}} {2a} $  编程计算，由于存在着同号两数相减的运算，所以在特殊情况下会产生较大误差而使得计算结果失去了意义。
  2. 描述误差
    • 为了便于数学分析和数值计算，人们对实际问题的数学描述通常只反映出主要因素之间的数量关系，而忽略次要因素的作用,由此产生的误差称为描述误差。
    • 举例：假设一门大炮的炮筒与水平面的夹角为θ，所发射的炮弹的出口速度为v0,那么t秒钟后炮弹的位置可表示为

      在这个模型中，我们忽略了空气的阻力。对实际问题进行数学描述通常称为是建立数学模型，所以描述误差也称为是模型误差。
  3. 观测误差
    描述实际问题或实际系统的数学模型中的某些参数往往是通过实验观测得到的。由试验得到的数据与实际数据之间的误差称为观测误差。
  4. 舍入误差
    几乎所有的计算工具，当然也包括电子计算机，都只能用一定数位的小数来近似地表示数位较多或无限的小数，由此产生的误差称为舍入误差。
  5. 截断误差
    假如真值 x_ 为近似值系列{xn}的极限,由于计算机只能执行有限步的计算过程，所以我们只能选取某个xN作为x _的近似值，由此产生的误差称为截断误差。

误差基本概念

  1. 误差分类
    1. 绝对误差
    2. 相对误差

两者定义：设 x 为某个量的真值，x 为 x 的近似值，称 |x - x|为近似值x的绝对误差；称|x - x|/|x*|为近似值x的相对误差

  1. 误差限
    既然在绝大多数情况下我们无法确定出真值x *，那么近似值x的误差、相对误差、以及绝对误差也都是无法确定的，但是我们总有办法估计出它们会落在那个范围之内。

    1. 定义  设x为真值 x 的近似值:
      • 若e>0满足条件|x*-x|≤e  , 则称e为x的绝对误差限（或误差限）
      • 若er>0满足条件|x-x|/|x|≤er,则称er为x的相对误差限.。
        提示：绝对误差限和相对误差限满足关系 e = er·|x|  所以只要知道其中的任意一个，即可求出另外一个。强调：我们平时口语中所讲的误差，实际所指的都是误差限。

有效数字

  • 重要性：有效数字的重要性在于，它大大简化了误差分析，使我们可以根据数值的表示形式即可简单地计算出绝对误差和相对误差，从而判断它的有效性。
  • 概念：如果真值  
    的近似值x的绝对误差限是它的某一个数位的半个单位，则称近似值x准确到这一位，且这一位一直到最左边第一个非零数字为止的所有数字都称为有效数字，有效数字的个数称为有效数字的位数。
    记号：若近似值x有p位有效数字，那么我们可以把x表示为
  • 计算
    若x *的近似值为  
     其中x1,x2,…,xp∈{0,1,…,9}

    • 绝对误差
      则有效位数为：
       (就是小数点最位置$0.0...1/2 ==0.5 * 10^{后面的n-原本的数小数点有多少个数字} $)
    • 相对误差
      若x*的近似值x有p位有效数字，那么我们有
  • 例子

    例：设

     分别求出x,y的绝对误差和相对误差（限）

解：

  • x,y的绝对误差分别为

    ![](https://www.yuque.com/api/services/graph/generate_redirect/latex?ex%3D0.00005%3D0.5%C3%9710%5E%7B-4%7D#card=math&code=ex%3D0.00005%3D0.5%C3%9710%5E%7B-4%7D&id=fXPMx)       $$ey=0.0000005×10-9=0.5×10^{-15}$$

  • x,y的相对误差分别为

利用微分估算误差

绝对误差

计算方法课程中，我们可以这样来理解上面的一阶泰勒展式：假设f是一个系统，x是我们希望输入的真值，从而得到希望输出的真值f(x)。但是实际输入的是x的近似值x=x+Δx,其误差为Δx，如果用f(x)作为f(x*)的近似值，则误差为      Δy =f'(x)·Δx + o(Δx)  略去比Δx高级的无穷小项，两边取绝对值

结论：在微分公式df(x)= f'(x)dx中如果把|dx|解释为x的绝对误差（限），那么|df(x)|就是所得到的计算结果f(x)的绝对误差限。

则有:

 ![](https://www.yuque.com/api/services/graph/generate_redirect/latex?%7C%5Cdelta%20y%7C%E2%89%88%7Cf'(x)%7C%C2%B7%7C%5Cdelta%20x%7C#card=math&code=%7C%5Cdelta%20y%7C%E2%89%88%7Cf%27%28x%29%7C%C2%B7%7C%5Cdelta%20x%7C&id=YCWrz)

其中：

是误差，
是求导带入x计算后的数据，
是要求的绝对误差限

结论：在微分公式df(x)= f'(x)dx中如果把|dx|解释为x的绝对误差（限），那么|df(x)|就是所得到的计算结果f(x)的绝对误差限。

例题：

例4.1：要做一个体积为1000

的正方体，要求体积的绝对误差不超过0.3
,试求边长的绝对误差限。

解： 设边长为x,体积为v，

        由 v=x3=1000 得x=10(cm)
        ![](https://www.yuque.com/api/services/graph/generate_redirect/latex?dv%20%3D%203x%5E2%20dx%20%3D%20300%C2%B7dx#card=math&code=dv%20%3D%203x%5E2%20dx%20%3D%20300%C2%B7dx&id=mlDDa)
        由300·dx=0.3 得 dx=1/1000

答：边长的绝对误差限为1/1000 (cm)

看完下面相对误差之后，可以求边长的相对误差限

有上面得到 Δx = 0.001

用公式可以的 得出  

 = 0.0001 cm$

相对误差

当我们把f(x)的绝对误差（限）近似地表示为df(x)时，我们可以进一步得到f(x)的相对误差（限）

结论：我们可以把|dln[f(x)]|作为f(x)的相对误差限的计算公式。

则有：

    ![](https://lsky-1316932265.cos.ap-guangzhou.myqcloud.com/uPic/image-20230530112249184.png#id=fjpFN&originHeight=208&originWidth=908&originalType=binary&ratio=1&rotation=0&showTitle=false&status=done&style=none&title=)

例题：

例4.2：设某机器上一个圆形的铁片零件的半径的设计要  求为 r =100±0.5 mm  试求这个圆形铁片的面积的绝对误差限和相对误差限。

解：

  1. 绝对误差：
    已知r =100，dr=0.5,由S=π·r2
    得： |ds| =|2π·r·dr|=2π·100·0.5=314 ( mm2)
  2. 相对误差

答：绝对误差限为$ 314 mm^2$，相对误差限为 1/100。

提示：绝对误差限是有量纲的，也就是要带单位，而相对误差限则是不带单位的数

选用计算方法应遵循的规则

  • 由于现在计算机的性能有大幅度提高，所以应当尽量采用双精度数计算，把提高精度放在更为重要的位置。
  • 应当尽量避免绝对值相近的两个正数相减，这样会吃掉许多有效数字位，这是计算方法中的头号杀手。
  • 应当尽量避免绝对值较大的数除以绝对值较小的数，如前所述这样会产生较大的绝对误差，甚至导致计算机数字溢出。
  • 应当防止“大数”吃“小数”绝对值相差很大的两个正数相加或相减，由于是对位相加，所以小数会吃亏。大量的这样类似的计算会有些问题。
  • 尽可能减少运算次数，防止误差扩散。

常用函数值计算方法

例子引入

某函数：

表示为  

A=(a0,a1,…,an)    程序设计时，可以把A说明为一个n+1维数组，此时约定A[0]存放a0，A[1]存放a1，…，A[n]存放an。

编程函数后： double PolyValue(double x,double*A,int n)

其中x为自变量；*A为A[0]的地址；n为多项式的次数。

逐项求和算法

求PolyValue(x,A,n)的数值解最容易想到的方法是按次数由低到高的顺序逐项求和。

为此可以在程序中说明一个临时变量power存放x的各次幂，并利用xk=xxk-1来简化计算。

double PolyValue(double x, double *A, int n)
{
  double power = 1.0, y;
  int k;
  y = A[0];
  for (k = 1; k <= n; k++)
  {
    power *= x;
    y += A[k] * power;
  }
  return y;
}

秦九韶算法

可以把多项式改写为便于递推的形式

A(x)=1+2x+3x2+4x3+5x4

A(x)=1+x(2+3x+4x2+5x3)=1+x(2+x(3+4x+5x2))

    =1+x(2+x(3+x(4+5x)))=1+x(2+x(3+x(4+x(5))))

即：yk=(ak+x(ak+1+x(…(an-1+x(an))…))). 约定yn=an，不难得到递推关系式，且y0就是所需要的结果。

y4=5 ；

yk=

，k=n-1,n-2,…,0

y3=$4+xy_4=4+x(5) $

y2=

#card=math&code=3%2Bxy_3%3D3%2Bx%284%2Bx%285%29%29&id=cavWs)

y1=

)#card=math&code=2%2Bxy_2%3D2%2Bx%283%2Bx%284%2Bx%285%29%29%29&id=Bfreu)

y0=

))#card=math&code=1%2Bxy_1%3D1%2Bx%282%2Bx%283%2Bx%284%2Bx%285%29%29%29%29&id=lO86V)

可以把多项式改写为便于递推的形式

%E2%80%A6))%E2%80%A6))#card=math&code=A%28x%29%3D%28%20a0%2B%20x%28%E2%80%A6%28ak%2Bx%28%E2%80%A6%28a%7Bn-1%7D%2Bx%28an%29%29%E2%80%A6%29%29%E2%80%A6%29%29&id=RRhOm)

即：

%E2%80%A6)))#card=math&code=y_k%3D%28a_k%2Bx%28a_k%2B1%2Bx%28%E2%80%A6%28an-1%2Bx%28an%29%29%E2%80%A6%29%29%29&id=zs6wq)

约定yn=an，不难得到递推关系式

double PolyValue(double x, double *A, int n)
{
  double y;
  int k;
  y = A[n];
  for (k = n - 1; k >= 0; k--)
    y = y * x + A[k];
  return y;
}

由于不需要保留yn，yn-1，…，y1等中间结果，所以程序中只用一个变量y来动态地表示它们。秦九韶算法的循环体内只有一次乘法运算，所以算法的计算量为O(n)，相当于程序3.01的一半。避免了大数吃小数的问题。

三角函数值计算方法

  1. 首先考虑正弦函数y=sin(x)在x∈[0,π]的函数值计算问题。利用y=sin(x)在x=0处的泰勒展式
  2. 转化得到

    对于x∈[0,π]而言，上面两式的收敛性都不成问题。对于求sin(x)的近似值来说，无论是计算量，还是误差控制，利用后面一个公式计算的优越性更多一些

步骤

  1. 记SINTNV(x,n)表示利用sin(x)/x的泰勒展式取前n+1项之和作为sin(x)/x的近似值

    得到：

    不难理解，在一定条件下，n值取得愈大，利用上面两式得到的sin(x)的近似值的精度愈高。
  2. 改写公式

    改写后得到以下递推公式

    而且有

例子

SINTNV

代码

double SINTNV(double x, int N)
{
  int K;
  double y = 1.0, xx = 0.0, temp;
  K = N * 2;
  xx = x * x;
  while (K > 0)
  {
    y = 1.0 - y * xx / K / (K + 1);
    K -= 2;
  }
  return y * x;
}

例题：

对数函数值计算方法

介绍

  1. 求对数函数值最关键的问题实际上是计算[1,2]内的数的自然对数。
  2. 对于计算大于2的数x的对数来说我们总可以找到整数K,使得2K<x< 2K+1，记u=x/2K,则有1<u<2,
    $ln(x) =ln(2^Ku)=Kln(2)+ln(u) $$
    此时还有1<u<2。

对于计算(0,1)内的数的对数来说，我们也可以利用

把它化为求大于1的数的对数问题。

求一般的对数来说，我们总可以利用换底公式转化为自然对数的计算问题。

计算方法

对数函数泰勒展式的处理方法

  1. 对数函数ln(1+x)在x=0处的泰勒展式

    只是在(-1,1)内处处收敛，但 收敛速度太慢，所以对于数值计算来说，实际意义不大。
    为此，我们在写出
    在x=0处的泰勒展式

    它的收敛区域为(-1,1)。
    将上面两式相减，再做点数学处理即得

    而且收敛区间为(-1,1)。
  2. 简单的数学处理方法
    1. 对适当给定的x，记LOGTNV(x,n)表示利用上面的(*)式右边取前n+1项和所得到的结果
    2. 假如要计算某个正数a(不妨假设a>1)的对数的近似值，我们有
    3. 而且对任意a>1，均有0<x<1，从而能保证(*)式右边收敛，所以对于适当的n值，我们有
  3. LOGTNV(x,n)的数值计算方法
    1. 模仿秦九韶算法的处理方法，可以先把LOGTNV(x,n)表示为

      不难得到下面的递推关系

代码

double LOGTNV(double x, int n)
{
  int K, NK;
  double xx, y;
  NK = n * 2 + 1;
  xx = x * x;
  y = 1.0 / NK;
  for (K = n; K > 0; K--)
  {
    NK -= 2;
    y = 1.0 / NK + xx * y;
  }
  return y;
}

int main()
{
  double a;
  double x, y0, y1, y2, y3, y4;
  a = 10;
  x = (a - 1) / (a + 1);
  y0 = 2.0 * x * LOGTNV(x, 84);
  y1 = log(a);
  printf("\ny0=%24.16e", y0);
  printf("\ny1=%24.16e", y1);
  getchar();
}

例题

指数函数与幂函数值计算方法

对于指数函数ex的数值计算问题来说，也可以利用ex在x=0处的泰勒展式

来帮助解决。由于

，所以不失一般性，可以假设0<x<2(
，
)

数值计算方法

  1. 改写成
  2. 可以立即得到如下的递推格式

代码

double EXPTNV(double x, int N)
{
  int K = 0;
  double y = 1.0;
  for (K = N; K > 0; K--)
    y = 1.0 + y * x / K;
  return y;
}

int main()
{
  double x, y0, y1;
  x = 1.50;
  y0 = EXPTNV(x, 24);
  y1 = exp(x);
  printf("\ny0=%24.16e", y0);
  printf("\ny1=%24.16e", y1);
  getchar();
}

对数误差

EXPTNR代码

int EXPTRN(double x, double eps)
{
  int N = 1;
  double y;
  N = (int)(x + 1);
  while (1)
  {
    if (EXPTNR(x, N) < eps)
      break;
    N++;
  }
  return N - 1;
}

double EXPTNR(double x, int N)
{
  int K;
  double y;
  y = 1.0 / (1.0 - x / N);
  for (K = 1; K <= N; K++)
    y *= x / K;
  return y;
}

double EXPTV(double x)
{
  int N = 23;
  int K = 0;
  double y = 1.0;
  if (x <= 0.5)
    N = 14; /*See Ch3 Table 7.1 */
  else if (x <= 1.0)
    N = 18;
  else
    N = 23;
  for (K = N; K > 0; K--)
    y = 1.0 + y * x / K;
  return y;
}

例题

求函数的零点与极值问题

函数的零点与极值问题概述

  • 零点
    设f(x)是定义在闭区间[a,b]上的连续函数，如果x∈[a,b]使得f(x)=0,则称x*是f(x)的一个零点。
  • 极值点
    1. 设y= f(x)是定义在区间[a,b]上的连续函数，对于x∈(a,b)，如果存在 δ > 0 使得当
      时，恒有
      ，则称 x*是f(x)在(a,b)内的一个极小值点。
    2. x 是极值点的必要条件是
      。稍晚些时我们将会看到，对于凸函数来说，极小值点也是最小值点，从而上面的必要条件变成充要条件。

介值定理的应用

由高等数学中的介值定理可知，如果f(x)是连续函数，而且f(a)·f(b)<0,则方程f(x)=0在[a,b]内一定有解。

如图所示，只要f(a)和(b)  异号，函数的图像就一定  在[a,b]内和x轴相交，也  就是 f(x)=0在[a,b]内有解,  从而找到了隔根区间。

求隔根区间的方法

  1. 分析法：利用高等数学(数学分析)中的函数作图方法画出y=f(x)的大致曲线即可确定隔根区间;(学生最适用)
  2. 试算法:试探性地给出两个数a,b,若f(a)·f(b)<0，则试算成功，否则修改a或b,重新试算；(适用于数学经验丰富者)
  3. 搜索法：从某一点开始，按一定的步长分别向左或向右逐点计算，直到出现相邻两点函数值异号。 (适用于工程师）

区间对分法

对分法适用的范围最大，数值稳定性性最好，算法特别简单，而且事先可以控制精度，

算法原理

例题

黄金分割法

设y= f(x)是定义在某区间[a,b]上的实函数，如果存在

使得f(x)在
上是单调减函数，而在
上是单调增函数，则称f(x)是区间[a,b]上的单峰函数。

设f(x)是定义在[a,b]上的严格的单峰函数，x*是f(x)是在[a,b]上的最小值点，则对任意满足a<x1<x2<b的

我们有：

  1. 若
    ，则
  2. 若
    ，则
  3. 若
    ，则

计算方法

  1. 带入A[k]与B[k]计算两端的黄金分割点（题目会给计算式子）
    1. x1=a+0.618(x2-a)=a+0.382(b-a)
    2. x2=a+0.618(b-a)
  2. 计算后x1 x2后代入到f(x)中计算，得出y1[k]、y2[k]，并同时计算AX[k]
  3. 判断y1[k]和y2[k]谁大，谁大就把对应的x1,x2换成A[k]或B[k]

例题

多项式计算

引言

  1. 一般的实系数n次多项式可以表为

    其中
    称为k次项的系数，k=0，1，…，n。

    1. 为了强调求A(x)是真正的n次多项式，也可以把A(x)记为
      ，在这种情况下，总是假定
      ≠0。
    2. 若n=0而
      ，则称A(x)为零次多项式，此时A(x)成为非零常数。
    3. 若A(x)≡0，则称之为-1次多项式。
  2. 由于多项式的运算主要是利用相应的系数进行计算，所以只需把多项式的系数存贮在计算机内即可。
  3. 我们约定用多项式的函数名作为存贮多项式的系数的数组名使用，数组的长度比多项式的次数多1。
    举例：多项式
    在计算机内可以说明为
    double A[4]={1.0,2.0,3.0,4.0};
    约定：对于多项式相关的各种数值计算来说，当多项式的次数 n≤0 时并无实际意义。所以，如果没有特别说明，我们总可以假定n>0。

多项式的基本运算

多项式的乘法

多项式带余除法

若A(x)和B(x)分别为n次和m次多项式，且n≥m，则商式q(x)为n-m次多项式，系数的个数为n-m+1；而余式R(x)最多为m-1次，系数的的个数最多为m，

线性方程求解

线性方程组概述

由n个变元,m个线性方程组成的线性方程组的一般形式为:

对于上面给出的一般形式有线性方程组,我们可以用矩阵和向量的记号把它简写为Ax=b,其中

它们分别称为系数矩阵、未知向量、以及常数向量.

为了简化我们的讨论,对于一般形式的线性方程组而言,我们总可以假定它的系数矩阵或者满行秩,或者满列秩.实际情况也大都是这样的.

在上面的假定下,我们有

  1. 如果m<n,那么方程组有无数组解,我们通常要找某种意义下的最优解,这是运筹学中的线性规划就是专门研究这个问题；
  2. 如果m>n;那么方程组无解,此时作为实际问题还是有意义的,在求近似解我们可以转向求它的最小二乘解,这是下一章中讨论的问题.
  3. 如果A是满秩的n阶方阵,那么方程有唯一解.

高斯消去法

中学里学过的加减消元法.它包括两个过程:首先把线性方程组化为上三角形线性方程组

选主元高斯约旦消去法

  1. 无论是高斯消去法还是约当消去法解线性方程组,我们都要进行除法运算.当某个aKK的绝对值非常小时,这两种方法的数值稳定性可能不好,为此.我们可用选主元消去法.
  2. 选主元的方法有选列主元和全主元之分，基本原理是互换任意两个方程的位置不影响方程组的解；互换任意两个变量的位置也不影响方程组的解。
  3. 选列主元的程序相对说来简单得多，所以我们只推荐选列主元。
  4. 选列主元的思想就是把aKK,aK+1K,…,aMK中绝对值最大的元素移到主对角线上来.
  5. 俩者计算方法
    • 高斯-约旦(G-J)消元法
      （1）从第一行开始，先顺序寻找到主元（不能为0），并将主元通过行变换移动到主对角线上
      （2）将主元所在的行内的所有元素除以主元，使得主元化为1（也可以放在第四部之后）
      （3）其他行减去 主元所在行 乘以一定的倍数，使得主元所在的列内除主元外的其他元素为0
      （4）在剩下的矩阵范围内寻找下一个主元，重复以上步骤
    • 列主元消去法
      （1）从第一列行开始，选出当前列中绝对值最大的元素作为这一列的主元，并将主元通过行变换移动到主对角线上
      （2）将主元所在的行内的所有元素除以主元，使得主元化为1（也可以放在第四部之后）
      （3）将剩下矩阵范围的其他行减去 主元所在行 乘以一定的倍数，使得剩下矩阵 主元所在的列内除主元外的其他元素为0
      （4）在剩下的矩阵范围内寻找下一个主元，重复以上步骤
      （5）回带求解

高斯消去法应用模块的编写

两个算法的对比分析

  1. 计算量:高斯消去法为O(0.33N3),约当消去法为O(0.5N3),从现代的观点看,两者的数量级相同；
  2. 算法简单:约当消去法占优；
  3. 通用性:  约当消去法占优；
  4. 数字稳定性:约当消去法更易于解决.结论:如果用手工或计算器求解线性方程组（比如应付考试）,
  5. 用高斯消去法较好,如果编程用计算机求解线性方程组,则用约当消去法更好.

最小二乘法与曲线拟合

线性函数拟合方法

  1. 为了让大家对最小二乘法适用的实际问题有一个感性认识，我们首先设想一个简单的工程问题：如何尽可能准确地测算出一个弹簧的长度与所受到的外界的拉力之间的函数关系。我们知道，在弹性限度以内，弹簧的长度的改变量与所受到的外界的拉力成正比，亦即弹簧的长度Y与所受到的外界的拉力X满足一个线性关系:     Y=a0+a1x 我们要作的就是综合物理和数学方法求出a0和a1。为了解决这个问题，我们可按下面的步骤来进行。

最小二乘法

下一篇