第1章 数值计算中的误差
1.1 实验目的
了解MATLAB基本操作、数值计算中误差的传递,避免误差的方法。
1.2 MATLAB基本操作
矩阵运算
%%行向量
>> x=[-3 1 7] %也可用逗号分隔
x =
-3 1 7
%%列向量
>> y=[0 2 9]' %行向量转置
y =
0
2
9
%%矩阵
>> A=[1 2 3;4 5 6;7 8 9]%分号换行
A =
1 2 3
4 5 6
7 8 9
%%矩阵转置
>> B=A'
B =
1 4 7
2 5 8
3 6 9
%%矩阵加减法运算
>> C=A+B
C =
2 6 10
6 10 14
10 14 18
>> C=A-1
C =
0 1 2
3 4 5
6 7 8
%%矩阵乘法运算
>> C=C*2
C =
0 2 4
6 8 10
12 14 16
>> C=A*C
C =
48 60 72
102 132 162
156 204 252
%%逆
>> A=[1 2;3 5];
>> inv(A)
ans =
-5.0000 2.0000
3.0000 -1.0000
%%矩阵除法运算
>> A=[1 2 3;4 5 6;7 8 0];
>> b=[1 1 1]';
>> x=A\b %左除
x =
-1.0000
1.0000
-0.0000
此处矩阵的除法运算相当于下面方程组的解:
A
x
=
b
Ax=b
Ax=b
即:
x
=
A
−
1
b
x=A^{-1}b
x=A−1b
数组运算
数组的运算符号是点运算,即 .* ./ 或 .\ 还有 .^ 。(加减运算无点运算)
左除和右除的区别在于交换了被除数和除数。
%%数组输入
>> a=0:2:10
a =
0 2 4 6 8 10
>> a=linspace(0,10,6) %构造一维等差数组
a =
0 2 4 6 8 10
%%数组加减法与矩阵加减法运算相同
%%数组乘法
>> a=0:2:10
a =
0 2 4 6 8 10
>> b=0:2:10
b =
0 2 4 6 8 10
>> c=a.*b
c =
0 4 16 36 64 100
%%数组除法
>> a./b
ans =
NaN 1 1 1 1 1
%%数组乘幂
>> c=a.^2
c =
0 4 16 36 64 100
逻辑运算
与:&,或:|,非:~。
>> a=[-1 2 4;5 4 -8];
>> b=[-2 1 -1;3 -1 2];
>> c=(a>0)&(b>0)
c =
2×3 logical array
0 1 0
1 0 0
>> d=(a>0)|(b>0)
d =
2×3 logical array
0 1 1
1 1 1
>> e=~(a>0)
e =
2×3 logical array
1 0 0
0 0 1
基本函数
基本初等函数
sin(x) %正弦
cos(x) %余弦
tan(x) %正切
asin(x) %反正弦
exp(x) %e的指数
log(x) %自然对数
sqrt(x) %平方根
abs(x) %绝对值
与矩阵有关的常用函数
二范数:二范数指矩阵A的2范数,就是A的转置共轭矩阵与矩阵A的积的最大特征根的平方根值,是指空间上两个向量矩阵的直线距离。类似于求棋盘上两点间的直线距离。
条件数:数值分析中,一个问题的条件数是该数量在数值计算中的容易程度的衡量,也就是该问题的适定性。一个低条件数的问题称为良态的,而高条件数的问题称为病态(或者说非良态)的。条件数是线性方程组Ax=b的解对b中的误差或不确定度的敏感性的度量。数学定义为矩阵A的条件数等于A的范数与A的逆的范数的乘积,即cond(A)=‖A‖·‖A的逆‖,对应矩阵的3种范数,相应地可以定义3种条件数。
>> A=[1 -2;-3 4];
>> norm(A) %求二范数
ans =
5.4650
>> cond(A) %求条件数
ans =
14.9330
>> rank(A) %求矩阵的秩
ans =
2
>> zeros(2,3) %生成一个2行3列的零矩阵
ans =
0 0 0
0 0 0
>> ones(2,3) %生成一个2行3列元素为1的矩阵
ans =
1 1 1
1 1 1
>> eye(3) %生成3阶单位阵
ans =
1 0 0
0 1 0
0 0 1
>> length(A) %求向量长度维数(列数)
ans =
2
控制流语句
for循环语句
for 循环变量 = 初值:增量:终值 %增量缺省值为1
语句
end
while循环语句
while 条件 %条件成立时执行下面语句,否则跳过
语句
end
if条件语句(三种形式)
if 条件
语句
end
if 条件
语句1
else
语句2
end
if 条件1
语句1
else if 条件2
语句2
else
语句3
end
M文件
命令文件
假设有一个tset.m文件如下:
n=6;
pp=1;
for i=1:n
pp=pp*i;
end
在命令输入窗口执行以下操作会出现相应结果:
>>test
>>pp
pp =
720
函数文件
function 因变量 = 函数(自变量)
%说明语句
语句
%计算n!的M文件
function tt = prod(n)
tt = 1;
for i = 1:n
tt = tt*i;
end
简单绘图
画函数曲线以及散点图
命令格式:plot(x,y,‘s’)
x代表横坐标,y代表纵坐标,s是可选择参数(默认为蓝色实线),具体意义如下:
例1:画出正余弦曲线
>> title('sin(x)&cos(x)') %图形标题
>> t=0:0.01:2*pi; %变量范围
>> plot(t,sin(t),'K',t,cos(t),'R')
>> xlabel('x');ylabel('y'); %坐标轴题注
结果如下:
例2:画出y=x^2 (0 <= x <= 1)的散点图
>> x=0:0.1:1;
>> y=x.^2;
>> plot(x,y,'.');
>> xlabel('x');ylabel('y');
结果如下:
例3:在区域[-pi,pi]×[-pi,pi]上画以下函数的网格曲面。
z
=
s
i
n
x
c
o
s
y
z=sinxcosy
z=sinxcosy
>> x=-pi:0.3:pi;y=-pi:0.3:pi;
>> [X,Y]=meshgrid(x,y);Z=sin(X).*cos(Y);
>> mesh(X,Y,Z);xlabel('x');ylabel('y');zlabel('z');
结果如下:
1.3 实验1例题:数值计算中的误差分析
例1:防止大数吃小数。设a=1×10308,b=1.1×10308,c=-1.001×10^308,用MATLAB计算a+(b+c),(a+b)+c。
>> a=1e+308;b=1.1e+308;c=-1.001e+308;
>> r1=(a+b)+c
r1 =
Inf
>> r2=a+(b+c)
r2 =
1.0990e+308
在MATLAB中能表示的最大的正实数是realmax,即1.7977e+308。凡是比这个数大的数,MATLAB就用Inf表示(正无穷)。这个题中就是因为这个原因使得结合律不成立。
例2:分析下面程序,循环语句是否会结束?
a=2;
b=2;
while a+b~=a
b=b/2;
end
按理来说这个程序不会结束,因为b永远不可能为0。但是实际上这个程序很快就会结束,返回b的值为b=2.2204e-016,说明a远大于b的值时,b的值就会被忽略,这时a+b=a,程序结束。
例3:举例说明太小的数不能做除数。
计算以下表达式:
(
2
+
x
)
−
2
x
\frac{(2+x)-2}{x}
x(2+x)−2
>> x=3;y=((2+x)-2)/x
y =
1
>> x=1e-12;y=((2+x)-2)/x
y =
1.0001
>> x=1e-16;y=((2+x)-2)/x
y =
0
该表达式只要x不为零结果就一定等于1。当x越小的时候可以发现结果和1的偏差越大,所以太小的数不宜做除数。
分析下面程序,循环语句是否会结束?
a=2;
b=2;
while a+b~=a
b=b/2;
end
按理来说这个程序不会结束,因为b永远不可能为0。但是实际上这个程序很快就会结束,返回b的值为b=2.2204e-016,说明a远大于b的值时,b的值就会被忽略,这时a+b=a,程序结束。
例3:举例说明太小的数不能做除数。
计算以下表达式:
(
2
+
x
)
−
2
x
\frac{(2+x)-2}{x}
x(2+x)−2
>> x=3;y=((2+x)-2)/x
y =
1
>> x=1e-12;y=((2+x)-2)/x
y =
1.0001
>> x=1e-16;y=((2+x)-2)/x
y =
0
该表达式只要x不为零结果就一定等于1。当x越小的时候可以发现结果和1的偏差越大,所以太小的数不宜做除数。