计算方法

笔记

误差

  • 计算机的浮点表示及算术运算
  • 误差来源
  • 误差的基本概念
  • 误差分析

计算机的浮点表示及算术运算

实数x在计算机中被表示为

给定的二进制浮点数只能表示一定范围的有限点集合

浮点数的四则运算

加减法先对阶

特点

  • 大数吃小数

  • 接近的数相减

  • 除以绝对值小的数

以上三种会造成比较大的误差,所以我们要合理安排顺序

误差来源

实际问题-> 数学模型-> 数值方法-> 程序设计-> 计算机计算-> 解答

观测误差 模型误差 方法误差 舍入误差

误差的基本概念

设$x$为准确值, 是$x$的一个近似值,称为近似值的绝对误差,或简称误差

,并满足

则称为近似值的绝对误差限,或简称误差限

设$x$为精确值,$x^*$为近似值,则称比值

为近似值的相对误差,记左 (实际常常用代替分母的)

是近似值的误差限,则称

为近似值$x^*$的相对误差限,此时,有

如果

则说$x^*$近视表示x准确到$10^{-n}$位,并从此位起j直到最左边的非零数值之间的一切数值都称为有效数字,并把有效数值的位数称为有效位数

将x的近似值$x^*$表示位十进制浮点数的标志形式

如果

则说近似值$x^*$具有$n$位有效数值.这里n是正整数,m是整数

若近似值$x^*$具有n位有效数值,则其相对误差满足

反之,如果的相对误差满足

则$x^*$至少具有n位有效数值。

误差分析

可以将$e^*$看成为微分

可以得到算术运算的误差

自变量有误差时计算函数值时也会参数误差,可以以用Taylor展开式进行估计

线性方程数值解

  • 直接法

    假设过程中不产生舍入误差,经过有限次运算可得方程组的精确解的方法

    低阶稠密矩阵方程组

    分为消元法和三角法

  • 迭代法

    从解的某个近似值出发,通过构造无穷序列去逼近精确解的方法

    大型稀疏矩阵方程组

高斯消元法

基本思想 用以下三种变换将其化为阶梯新式求解

  1. 用一个非零的数乘某一个方程
  2. 把一个方程的倍数加到另一个方程上
  3. 互换两个方程的位置

三个初等变换总是把方程组变为同解的方程组

我们可以化为上三角,这样可以简单的求解(其过程比较单一)

消元回代

实际上就是

增广矩阵

计算次数

回代 $\frac{n(n-1)}{2}$

矩阵的三角分解

矩阵的LU分解

可以把高斯消元的过程描述称矩阵

我们可以让

Doolittle分解

不能使用Gauss矩阵构造,相同与直接消元

待定系数法求LU

A=LU

待定系数法和矩阵相等

按顺序求解

可得

对称正定矩阵的平方根法和LDLT分解

当A时对称正定矩阵时,存在非奇异的下三角矩阵L,使得

当限定L的对角元素为正时,该分解是唯一的

为了避免开平方根

我们可以将对称正定矩阵A分解为

D是对角矩阵,而且都是正数

可以得到

然后和LU分解类似

解三对角方程组的追赶法

主对角和次对角非零,其他都为0

消元法在计算机上的实现

主元选择的必要性

按自然顺序

  • 可能无法选择(为0)
  • 误差
选主元的方法

选列主元

在变换到第k步时,选择绝对值最大的作为新的主元

然后交换到当前行

可以避免除以太小数的问题

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
function [x]=Gausscolumn(A,b)
% Gauss elimination method according to Column
[ni,nj]=size(A);
for j=1:ni %Find the max element in Column
[tv,ti]=max(A(j:ni,j));
A([ti+j-1,j],:)=A([j,ti+j-1],:);
b([ti+j-1,j])=b([j,ti+j-1]);
for i=j+1:ni %Elimination Process
d=-A(i,j)/A(j,j);
A(i,:)=A(i,:)+d*A(j,:);
b(i)=b(i)+d*b(j);
endfor
endfor
x=zeros(size(b)); %Back Substitute Process
x(ni)=b(ni)/A(ni,ni);
for i=ni-1:-1:1
x(i)=(b(i)-sum(x(i+1:ni).*A(i,i+1:ni)'))/A(i,i);
endfor
endfunction

定理2.2 设A非奇异,则存在转置矩阵P,以及单位下三角阵L和上三角阵U,使得

全选主元

整个区域选择主元

迭代改善

$r^{(1)}$需要采用高精度的字长

行列式和逆矩阵的计算

LU分解求解逆矩阵

通过LU分解求解向量方程组$AX=E$,得到的解即为矩阵A的逆矩阵

多个方程组求解

列主元

向量和矩阵范数

向量范数

定义 设$x\in R^n,N(x) \equiv |x|$为$x$的实直函数,若它满足下列条件

则称$N(x) \equiv |x| 为R^n$上的一个向量范数(或向量模),$|x|$的值成为向量$x$的范数

由定义可以推导出

三种范数

容易得到

等价性

设$|x|_s,|x|_t$为$R_n$上向量的任意两种范数,则存在常数$c_1,c_2>0$,使得对于一切$x\in R^n$ 有$c_1|x|\leq |x|_t \leq |x|_s$

向量序列收敛性

设$\{x^{(k)}=(x_{1}^{(k)} ,x_{2}^{(k)},x_{3}^{(k)} \cdots ,x_{n}^{(k)})^T\} $为$R^n$上的一个向量序列, 如果对于

则称向量序列$\{x^{(k)}\}$收敛于向量$x^*$

容易证明,$x^{(k)}$收敛于向量$x^*$的充要条件是

矩阵范数

对于实数,可以用绝对值度量其大小,对于复数可以用模度量它的大小,对于矩阵,可以用矩阵范数来度量它的大小

设$A\in R^{n\times n}\quad N(A)\equiv|A|$为A的实值函数,如果它满足下列条件

则称$N(A)\equiv|A|为R^{n\times n}$上的一个矩阵范数(或矩阵模),$|A|$的值成为矩阵$A$的范数

相容范数

算子范数

最常用的是利用向量范数来定义的矩阵范数

称之为矩阵A的算子范数,其中$p=1,2,\infty $

由上定义的矩阵范数为相容范数

算子范数满足

证明

2

谱半径

矩阵A的特征值按模最大值称为A的谱半径记作$\rho(A)$,即

由谱半径的定义,矩阵的2范数可以记为

当A是实对称矩阵时,有

矩阵的条件数与病态方程组

如果矩阵A或常数项b的微小变化,引起方程组$Ax=b$的矩阵变化,则称此方程组为‘病态’方程组,矩阵A称为‘病态’矩阵,否则称为良态

条件数

设$A$为非奇异矩阵,用$x$表示$Ax=b$的精确解,而$\widetilde{x}$时扰动方程组

的精确解.由$A=bx$两式相减,可知解的误差满足

利用矩阵的范数性质,有

又由$Ax=b$, 有

进一步考虑当A发生扰动时,方程组$Ax=b$的解的误差。设$\tilde{x}$是下列扰动方程组的

的精确解,由$Ax=b$和两式相减,可知解的误差满足方程

由此

于是有

定义 设A$\in R^{n\times n}$为可逆矩阵,称

为矩阵A在范数$|\cdot|_p$意义下的条件数

条件数越大,方程组的扰动误差对方程组的影响越严重

矩阵的平衡

对于条件数过大的矩阵,求解方程组

是困难的,为了改善系数矩阵的条件数可以用适当的非奇异对角矩阵$D_1$左乘上式的两端,得到

这叫做按行平衡

进一步还可以进行按列平衡,用非奇异对角矩阵$D_2$右乘矩阵$A$ 转而求解方程组

由$x=D_2y$求得x,选择对角阵$D_1,D_2$的标准就是使矩阵$D_1AD_2$的条件数尽可能的小

迭代法

比较适合打型稀疏系数矩阵

一般方法

考虑求解线性代数方程组

为了采用迭代法,首先要将方程组改写为等价形式

其中$H\in R^{n\times n}$ ,$H=[h_{ij}]$,$g\in R^n$ 为已知向量,

$x\in R$ 代表未知向量

给定初始近似值$x^{(0)}\in R^n$, 定义向量序列

则称$\{x^{(k)}\}$为迭代序列,并称$H$为迭代矩阵

如果当$k\rightarrow \infty$ 时,${x^{(k)}}$极限存在,设$x^*\in R^n$

则有

收敛性

利用迭代公式构造的序列$\{x^{(k)}\}$,以求得方程组的近似解的算法称为解式的简单迭代法,如果迭代序列$\{x^{()k}\}$收敛,久称此方法式收敛的

我们要关心的式是否能够收敛

即$x^{(k)} \rightarrow x^*$ 的充分必要条件是$H^{(k)}\rightarrow 0,k\rightarrow \infty$

引理1 ,对任意值收敛的充分必要条件是

引理2,

$H^{k}\rightarrow0(k\rightarrow\infty)$的充分必要条件是$H$的普半径

若$|H|< 1$ ,$(|\cdot|$允许为任何一种相容的矩阵范数$)$ 则迭代法收敛

一般步骤
  • 依据方程组分离$x$
  • 得到迭代格式
  • 判断迭代格式是否收敛
  • 迭代求解
  • 满足终止条件,迭代结束

定理

当$|H|\leq 1$ 时,由迭代式所定义的序列满足如下估计式

雅可比迭代法

考虑方程组

其中$A\in R^{n\times n}$ 是非奇异的,$b\in R$ 为已知向量,将矩阵A写成如下

其中$D=diag(a_{11},\cdots,a_{nn})$为对角阵,$-L,-U$分别为$A$的严格下,上三角部分构成的三角形

当D非奇异,即$a_{ii}\neq 0(i=1,2,\dots)$时,可将方程组$Ax=b$写成

可以得到迭代格式

此格式为求解方程组$Ax=b$的$jacobi$迭代法

注意到$L+U=D-A$ 也可以写为

$jacobi$方法的迭代矩阵为

$jacobi$迭代法的分量形式为

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
function [x]=jacobi(A,b)
n=size(A,1);
x=[0;0;0]
normVal=Inf;
JacobItr=0;
eps=1.0e-5;
while normVal>eps
xold=x;
for i=1:n
sigma=0;
for j=1:n
if j~=i
sigma=sigma+A(i,j)*xold(j);
%sigma=sigma+A(i,j)*x(j);
endif
endfor
x(i)=(1/A(i,i))*(b(i)-sigma);
endfor
JacobItr=JacobItr+1;
x=x
normVal=norm(xold-x);
end
endfunction
b=[25;10;14];
A=[20,2,3;1,8,1;2,-3,15];
x=jacobi(A,b);
高斯-塞德尔迭代法

简单迭代法的分量形式时

可以用这些新值来计算$x_i^{(k+1)}$于是可以得到迭代格式

此方法称为$Seidel$迭代法

是对于迭代的一种策略,而非具体的某种迭代方法

对Jacobi迭代运用$Seidel$技巧得到

称为$Gauss-Seidel$迭代法,其矩阵形式为

整理可得

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
 b=[25;10;14];
A=[20,2,3;1,8,1;2,-3,15];
function [x]=Gauss_Seidel(A,b)
n=size(A,1);
x=rand(n,1);
normVal=Inf;
eps=1.0e-5;
GaussItr=0;
for i=1:n
x(i)=0;
endfor
while normVal>eps
x_old=x;
for i=1:n
sigma=0;
for j=1:i-1
sigma=sigma+A(i,j)*x(j);
end
for j=i+1:n
sigma=sigma+A(i,j)*x_old(j);
end
x(i)=(1/A(i,i))*(b(i)-sigma);
end
GaussItr=GaussItr+1;
x
normVal=norm(x_old-x);
end
endfunction
Gauss_Seidel(A,b)

两种方法收敛性无关

解非线性方程(组)的数值法

引入

非线性方程组问题

的求根问题,其中f(x)为非线性函数

给定如下非线性方程组

引入向量,向量函数记号

则方程组可改为

定义 设有使得则称为方程的根或零点,若存在正整数m,使得

,则称为m重根。当m=1时,为单根,这时满足条件

二分法

介值定理 如果函数$f(x)$在$[a,b]$上连续,且$f(a)f(b)\leq0$,则至少有一个数$\xi\in[a,b]$使得$f(\xi)=0$

若同时$f(x)$的一阶导数$f’(x)$在$[a,b]$内存在且保持定号(单调),则这样的$\xi$在$[a,b]$内唯一

基本思想

首先确定有根区间,将区间二等分,通过判断f(x)的符号,逐步将有根区间缩小,直到有限根区间足够地小,便可以求出满足精度要求的近似根

1
2
3
4
5
6
7
8
9
10
11
12
13
14
function [mid] = Bisection(f,a,b,eps)
% Bisection method: find zero point of nonlinear function f
% within [a, b], with precision of eps
if a<b
while (b-a)>eps
mid = (a+b)/2;
if (f(a)*f(mid)<0)
b = mid;
else
a = mid;
endif
endwhile
endif
endfunction

不动点迭代法及其收敛性

不动点迭代法是一种逐次逼近的方法,它的基本思想是通过构造一个递推关系式(迭代格式),计算出根的近似值序列,并要求该序列收敛于方程的根

将方程$f(x)$改写成等价形式

称所有满足方程$x=\varphi(x)$的点$x$为次方程的不动点

取初始迭代值$x_0$,记$x_1=\varphi(x_0),x_2=\varphi(x_1),\dots$, 一般有

记数列$x_0,x_1,\dots,x_k,\dots$为$\{x_k\}$,当k充分大时,如果$\{x_k\}$有极限存在,即

则对两边同时取极限,得到

这表明为等价方程$x=\varphi(x)$的不动点,即$x^*$为方程$f(x)=0$的根

  • 化等价方程,构造迭代函数
  • 研究迭代数列$\{x_k\}$的收敛性,收敛速度及误差估计

收敛性定理

假定方程$x=\varphi(x)$满足:

1) $\varphi(x)$在$[a,b]$上连续;

2) 当$x\in [a,b]$时,$\varphi(x)\in[a,b]$;

3) $\varphi’(x)$存在,且对任取$x\in [a,b]$有

则方程$x=\varphi(x)$在$[a,b]$上存在唯一不动点$x^*$,且对任取$x_0\in [a,b]$,由迭代格式

所定义$\{x_k\}$收敛于唯一不动点,且对任取$x_0\in [a,b]$,即同时下列误差估计式成立

证明:

条件三可以放宽为

若$\varphi(x)$在$[a,b]$满足$Lipschiz$条件,即对$[a,b]$上的任意两点$x_1,x_2$,有

且$Lipschiz$常数$L<1$,则定理依然成立

1
2
3
4
5
6
7
8
9
10
%fixed-point
function [x2] = FixedPoint(phi,a,b,eps)
% Newtow method to find the zero point of function f within [a, b]
% with the precision of eps. fd is the derivation funciton of f
while (abs(b-a)>eps)
a = b;
b = phi(a);
endwhile
x2 = b;
endfunction

收敛阶

设迭代格式

收敛于方程$x=\varphi(x)$的不动点$x^*$,如果迭代误差

满足渐进关系

则称次迭代式是p阶收敛的

Newton迭代法

将非线性函数$f(x)$逐步线性化,从而将非线性方程$f(x)=0$近似地转换为求线性方程求解

设已直方程$f(x)=0$有近似根$x_k$(假定$f’(x_k)\neq 0$) 将函数在点$x_k$处泰勒展开

方程$f(x)=0$可以近似地表示成

其根为

记根为$x_{k+1}$,则有

称上式为求方程$f(x)=0$根的Newton迭代法,也称为切线法

设函数$f(x)$在根附近有二阶连续导数,且$f’(x^)\neq 0$,则存在$x^$的临域

使对于任意$x_0\in D$,由迭代公式所定义的所定义的序列$\{x_k\}$收敛到方程$f(x)=0$的根,且由

平方收敛

证明

假定方程$f(x)=0$在区间$[a,b]$上有二阶连续函数,且满足条件
1.$f(a)f(b)<0$

2.$f’(x)\neq 0,x\in [a,b]$

3.$f’’(x)$在$[a,b]$上不变好。那么对于任意$x_0\in [a,b]$,只要满足

则以$x_0$为初值所产生的Newton迭代序列$\{x_k\}$必定收敛到$f(x)=0$的唯一实根

1
2
3
4
5
6
7
8
9
function [x2] = Newton(f, fd, a,b,eps)
% Newtow method to find the zero point of function f within [a, b]
% with the precision of eps. fd is the derivation funciton of f
while (abs(b-a)>eps)
a = b;
b = a - f(a)/fd(a);
endwhile
x2 = b;
endfunction

割线法

为了避开求导数,可以考虑用差商代替微商,从而避免了复杂的导数计算,利用相邻两次迭代的函数值做差商得

带入Newton迭代公式后,得到

1
2
3
4
5
6
7
8
9
function [x2] = Secant(f, a,b,eps)
% Secant method to find the zero point of function f, with precision of eps.
% a and b is the initial two points
while (abs(b-a)>eps)
x2 = b - f(b)*(b-a)/(f(b)-f(a));
a = b;
b = x2;
endwhile
endfunction

解非线性方程组的Newton迭代法

考虑非线性方程组

是它的一个精确解

作业

image-20200517144549130

image-20200517144621460

image-20200528201803308

image-20200603145901857

可以构造出

1
2
3
4
5
a =  2
a = 1.5000
a = 1.4167
a = 1.4142
ans = 1.4142
1
2
3
4
5
6
7
8
9
10
11
12
phi=@(x)(0.5*(x+2/x))
%fixed-point
function [x2] = FixedPoint(phi,a,b,eps)
% Newtow method to find the zero point of function f within [a, b]
% with the precision of eps. fd is the derivation funciton of f
while (abs(b-a)>eps)
a = b
b = phi(a);
endwhile
x2 = b;
endfunction
FixedPoint(phi,1,2,0.00005)

Newton迭代法可得

1
2
3
4
5
a =  3
a = 2.2917
a = 1.9655
a = 1.8844
a = 1.8794
1
2
3
4
5
6
7
8
9
10
11
12
f=@(x)(x*x*x-3*x-1);
f_1=@(x)(3*x*x-3);
function [x2] = Newton(f, fd, a,b,eps)
% Newtow method to find the zero point of function f within [a, b]
% with the precision of eps. fd is the derivation funciton of f
while (abs(b-a)>eps)
a = b
b = a - f(a)/fd(a);
endwhile
x2 = b;
endfunction
x=Newton(f,f_1,1,3,0.00005);

用割线法可得

1
2
3
a =  1.9000
a = 1.8811
a = 1.8794
1
2
3
4
5
6
7
8
9
10
11
12
f=@(x)(x*x*x-3*x-1);
f_1=@(x)(3*x*x-3);
function [x2] = Secant(f, a,b,eps)
% Secant method to find the zero point of function f, with precision of eps.
% a and b is the initial two points
while (abs(b-a)>eps)
x2 = b - f(b)*(b-a)/(f(b)-f(a));
a = b
b = x2;
endwhile
endfunction
x=Secant(f,2,1.9,0.00005);

image-20200603150927965

所以1 发散,其他应该收敛

实验可得

代码

1
2
3
4
5
%deri=@(x)(4*x-3+cos(x)+exp(x));
phi_1=@(x)((x*x+2)./3);
phi_2=@(x)(sqrt(3*x-2));
phi_3=@(x)(3-2./x);
phi_4=@(x)((x*x-2)./(2*x-3));
1
2
3
4
5
6
7
8
9
10
%fixed-point
function [x2] = FixedPoint(phi,a,b,eps)
% Newtow method to find the zero point of function f within [a, b]
% with the precision of eps. fd is the derivation funciton of f
while (abs(b-a)>eps)
a = b;
b = phi(a);
endwhile
x2 = b;
endfunction

选择

1
2
3
a=1.1;
b=3.0;
eps=0.0001;
1
2
3
4
x_1=FixedPoint(phi_1,a,b,eps);
x_2=FixedPoint(phi_2,a,b,eps);
x_3=FixedPoint(phi_3,a,b,eps);
x_4=FixedPoint(phi_4,a,b,eps);
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
%deri=@(x)(4*x-3+cos(x)+exp(x));
phi_1=@(x)((x*x+2)./3);
phi_2=@(x)(sqrt(3*x-2));
phi_3=@(x)(3-2./x);
phi_4=@(x)((x*x-2)./(2*x-3));
%fixed-point
function [x2] = FixedPoint(phi,a,b,eps)
% Newtow method to find the zero point of function f within [a, b]
% with the precision of eps. fd is the derivation funciton of f
while (abs(b-a)>eps)
a = b;
b = phi(a);
endwhile
x2 = b;
endfunction
a=1.1;
b=3.0;
eps=0.0001;
x_1=FixedPoint(phi_1,a,b,eps);
x_2=FixedPoint(phi_2,a,b,eps);
x_3=FixedPoint(phi_3,a,b,eps);
x_4=FixedPoint(phi_4,a,b,eps);

得到同样的结果

1
2
3
4
x_1 =  Inf
x_2 = 2.0002
x_3 = 2.0001
x_4 = 2.0000