Transcript E:WHU教学2010Lecture2008Lec6MATLAB_CG

MATLAB小结、
经典迭代法、CG
1 .MATLAB 代表MATrix LABoratory
• 它的首创者是美国新墨西哥大学计算机系的系主任
Cleve Moler博士，他在教授线性代数课程发现其
他语言很不方便，篇构思开发了MATLAB。最初采用
FORTRAN语言编写，20世纪80年代后出现了MATLAB
的第二版，全部采用C语言编写.
•
1984年Moler博士和一批数学家及软件专家创建了
MathWorks公司，专门开发MATLAB。
•
1993年出现了微机版，到2003年是6.5版
2 .一种演草纸式的科学计算语言
3 .MATLAB 是一高性能的技术计算语言.
–
–
–
–
强大的数值计算和工程运算功能
符号计算功能
强大的科学数据可视化能力
多种工具箱
MATLAB 能干什么？
MATLAB可以进行：
• 数学计算、算法开发、数据采集
• 建模、仿真、原型
• 数据分析、开发和可视化
• 科学和工程图形应用程序的开发，包括图形用户
界面的创建。
MATLAB广泛应用于：
• 数值计算、图形处理、符号运算、数学建模、系
统辨识、小波分析、实时控制、动态仿真等领域。
掌握 MATLAB ……
MATLAB的构成：
• MATLAB开发环境：进行应用研究开发的交互式平台
• MATLAB 数学与运算函数库：用于科学计算的函数
• MATLAB 语言：进行应用开发的编程工具
• 图形化开发：二维、三维图形开发的工具
• 应用程序接口 (API)：用于与其他预言混编
• 面向专门领域的工具箱：小波工具箱、神经网络工具
箱、信号处理工具箱、图像处理工具箱、模糊逻辑工
具箱、优化工具箱、鲁棒控制工具箱等几十个不同应
用的工具箱。
开发环境
包括：命令窗口、图形窗口、编辑窗口、
帮助窗口。
命令窗口
– 可在提示符后输入交互式命令
– 结果会自动的产生
– 例如：
command (typed at prompt)
MATLAB output
MATLAB prompt (>>) and cursor (|)
图形窗口
在窗口中输入：
• Plot([1,2,4,9,16],[1,2,3,4,5])
• MATLAB 划出如下图形:
编辑窗口
– 用来创建和修改M-files (MATLAB 脚本)
帮助窗口
The MATLAB Language
MATLAB 语言的特点
– Matlab的基本数据单元是不需指定维数的矩阵。
– Matlab的所有计算都是通过双精度进行的，在
内存中的数都是双精度的。
– double 是一个双精度浮点数，每个存储的双
精度数用64位。
– char用于存储字符，每个存储的字符用16位。
MATLAB的程序构成：
程序
•M文件与m函数
•流程控制
•函数
•语句
变量
•各种运算符
•图形显示
其它输出
常变量及其命名规则
• 变量名可以有数字、字母、下划线构成；
• 变量的首字符必须是字母；
• 区分变量名的大小写
• 每个变量名最长只能包含19个字符。
Matlab中预定义变量
•
•
•
•
•
•
ans 分配最新计算表达式的值，这个表达式并没有
给定一个名字
eps
返回机器精度
realmax 返回计算机能处理的最大浮点数
realmin 返回计算机能处理的最小的非零浮点数
pi ,3.14159265
inf 定义为1/0 。当出现被零除时，Matlab就返回inf，
并不中断执行而继续计算
NaN 定义为“Not a Number”，这个非数值要么是
％类型，要么是inf/inf
向量的创建
• 在matlab的命令窗口键入以下字符
• >> a = [1 2 3 4 5 6 9 8 7]
•
a=
•
1 2 3 4 5 6 9 8 7
• 希望得到元素从0到20，步距为2的一个向量，只需
•
•
•
键入以下命令即可
>> t = [0:2:20]
t=
2 4 6 8 10 12 14 16 18 20
矩阵的创建
输入矩阵时每一行元素有分号或者回车键
分隔。例如：
•
•
•
•
•
•
B = [1 2 3 4;5 6 7 8;9 10 11 12]
B=
1 2 3
5 6 7
9 10 11
4
8
12
各
种
运
算
符
语句
Matlab语言最基本的赋值语句结构为：
变量名列表=表达式
注1：整个赋值语句以；结束，则不在屏幕上返回
结果，否则立即返回结果。
注2：多个语句可在同一行，用逗号分开。
注3：表达是太长可以用续行符号…
函数
• Matlab由包括许多标准函数，每个函数都
完成某一特定功能的代码组成。
• Matlab也允许用户编写自己所需的函数，
其扩展名为.m，其中必须以关键字
function开头.
流程控制
•
•
•
•
•
•
•
•
循环语句
条件转移
开关语句
注释语句
续行 …
中断语句
暂停语句
回显语句
for, while
if end, if elseif else end
switch case
%
break
pause
echo on/off
1、for循环语句
基本格式
for 循环变量＝起始值：步长：终止值
循环体
end
步长缺省值为1，可以在正实数或负实数范围
内任意指定。对于正数，循环变量的值大于终
止值时，循环结束；对于负数，循环变量的值
小于终止值时，循环结束。循环结构可以嵌套
使用。
2、while循环语句
基本格式
while 表达式
循环体
end
• 若表达式为真，则执行循环体的内容，执行后
再判断表达式是否为真，若不为真，则跳出循
环体，向下继续执行。
While循环和for循环的区别在于，while循环结构的循环体被执
行的次数不是确定的，而for结构中循环体的执行次数是确定的。
3、if，else，elseif语句
（1）if
逻辑表达式
执行语句
end
（2）if 逻辑表达式 （3） if 逻辑表达式1
执行语句1
执行语句1
else
elseif 逻辑表达式2
执行语句2
执行语句2
end
…
end
4、switch语句
switch 表达式（可以是标量或字符串）
case
值1
语句1
case
值2
语句2
….
otherwise
语句3
end
MATLAB程序的基本组成结构
％说明
清除命令：清除workspace中的变量和图形（clear,close）
定义变量：包括全局变量的声明及参数值的设定
逐行执行命令：指MATLAB提供的运算指令或工具箱
… … …
提供的专用命令
控制循环 ：
包含for,if then,switch,while等语句
逐行执行命令
… … …
end
绘图命令：将运算结果绘制出来
• 当然更复杂程序还需要调用子程序，或与simulink以及其
他应用程序结合起来。
MATLAB的程序类型
MATLAB的程序类型有三种，一种是在命令窗口下执行的
脚本M文件；另外一种是可以存取的M文件，也即程序文
件；最后一种是函数（function）文件。
1、脚本M文件
在命令窗口中输入并执行，它所用的变量都要在工作空间
中获取，不需要输入输出参数的调用，退出MATLAB后就
释放了。
2、程序M文件
• 以.m格式进行存取，包含一连串的MATLAB指令和必要的
注解。需要在工作空间中创建并获取变量，也就是说处理
的数据为命令窗口中的数据，没有输入参数，也不会返回
参数。
• 程序运行时只需在工作空间中键入其名称即可。
3、函数文件
与在命令窗口中输入命令一样，函数接受输入参数，然后执行
并输出结果。用help命令可以显示它的注释说明。
具有标准的基本结构。
（1）函数定义行（关键字function）
• function[out1,out2,..] = filename(in1,in2,..)
• 输入和输出（返回）的参数个数分别由nargin和nargout两个
MATLAB保留的变量来给出。
（2）第一行帮助行，即H1行
• 以（%）开头，作为lookfor指令搜索的行
（3）函数体说明及有关注解
• 以（%）开头，用以说明函数的作用及有关内容
（4）函数体语句
• 函数体内使用的除返回和输入变量这些在function语句中直接
引用的变量以外的所有变量都是局部变量，即在该函数返回之
后，这些变量会自动在MATLAB的工作空间中清除掉。如果希
望这些中间变量成为在整个程序中都起作用的变量，则可以将
它们设置为全局变量。
Graphics
MATLAB提供了丰富的绘图功能
help graph2d可得到所有画二维图形的命令
help graph3d可得到所有画三维图形的命令
1、基本的绘图命令
plot（x1,y1,option1,x2,y2,option2,…）
x1,y1给出的数据分别为x,y轴坐标值，option1为选
项参数，以逐点连折线的方式绘制1个二维图形；同
时类似地绘制第二个二维图形。
这是plot命令的完全格式，在实际应用中可以根据需
要进行简化。比如：
plot(x,y)；plot(x,y,option)
选项参数option定义了图形曲线的颜色、线型及标示
符号，它由一对单引号括起来。
2、选择图像
figure（1）；figure（2）；…；figure(n)
打开不同的图形窗口，以便绘制不同的图形。
3、grid on：在所画出的图形坐标中加入栅格
grid off：除去图形坐标中的栅格
4、hold on：把当前图形保持在屏幕上不变，同时
允许在这个坐标内绘制另外一个图形。
hold off：使新图覆盖旧的图形
5、设定轴的范围
axis（[xmin xmax ymin ymax]）
axis(‘equal’)：将x坐标轴和y坐标轴的单位刻
度大小调整为一样。
6、文字标示
text(x,y,’字符串’)
在图形的指定坐标位置(x,y)处，标示单引号括起来的
字符串。
• title(‘字符串’)
在所画图形的最上端显示说明该图形标题的字符串。
• xlabel(‘字符串’)，ylabel(‘字符串’)
设置x，y坐标轴的名称。
• 输入特殊的文字需要用反斜杠（\）开头。
7、legend(‘字符串1’,‘字符串2’,…,‘字符串n’)
• 在屏幕上开启一个小视窗，然后依据绘图命令的先
后次序，用对应的字符串区分图形上的线。
8、subplot（m，n，k）：分割图形显示窗口
m:上下分割个数，n:左右分割个数，k:子图编号
9、semilogx：绘制以x轴为对数坐标（以10为底），
y轴为线性坐标的半对数坐标图形。
semilogy：绘制以y轴为对数坐标（以10为底），
x轴为线性坐标的半对数坐标图形。
10、了解应用型绘图指令：可用于数值统计分析或离散
数据处理
bar（x,y）；hist（y,x）
stairs（x,y）；stem（x,y）
三维的绘图命令
 3 2
 2
1 T
T
例: 二次型对应的曲面 f ( x)  x Ax  b x, A  
, b   .

2
2 6
  8
[X Y]=meshgrid(-10:0.5:10);
Z=0.5*(3*X.^2 + 4*X.*Y+6*Y.^2)-2*X+8*Y;
surfc(X,Y,Z); colormap hsv
x = (a*(1-v/(2*pi)).*(1+cos(u)) + c) .* cos(n*v) ;
y = (a*(1-v/(2*pi)).*(1+cos(u)) + c) .* sin(n*v) ;
z = b*v/(2*pi) + a*(1-v/(2*pi)) .* sin(u) ;
Examples
• 绘图实例
• 函数分析
• 矩阵运算
• 线性方程组
• 曲线拟合
• 微分方程
绘图实例
函数分析
fplot('func',[-1 1.5])
%作图
result = func(0)
%求函数值
xsolve = fzero('func',3)
%求解
Xmin = fminbnd('func',0.5,1)
%求最小值
矩阵运算
• A = [1 2 3 ; 4 5 6 ; 7 8 9];
• B = [1 2 3 ; 4 5 6];
• C = [1 0 1 ; 0 2 3 ; 4 5 0];
•
•
•
•
•
•
•
•
•
expC = exp(C)
expM = expm(C)
logM = logm(expM)
detA = det(A)
traceA = trace(A)
BT = B'
invA = inv(A)
rankA = rank(A)
[EigenVectors,EigenValues] = eig(A)
线性方程组与特征值
A = [ 3 1 -1 ; 1 2 4 ; -1 4 5 ];
b = [ 3.6 ; 2.1 ; -1.4 ];
X = A\b
[EigenVectors,EigenValues] = eig(A)
曲线拟合
•
•
•
•
•
•
•
•
•
•
•
%一次多项是拟合
%已知离散点
x = [1 1.5 3 4 5 6 6.5 7 8];
y = [1.2 1 1.7 2.5 2 2.3 2.5 3
3.1];
%最小二乘拟合
p1 = polyfit(x,y,1);
y1 = polyval(p1,x);
plot(x,y1);
hold on
plot(x,y,'ro')
grid on
• %7次多项是拟合
• %已知离散点
• x = [1 1.5 3 4 5 6 6.5 7
•
•
•
•
•
•
•
8];
y = [1.2 1 1.7 2.5 2 2.3
2.5 3 3.1];
%最小二乘拟合
p7 = polyfit(x,y,7);
xi = 1:0.25:8;
yi = polyval(p7,xi);
plot(x,y,'*r',xi,yi);
grid on
微分方程
• Van der Pol Equation
2
d y
dy
2

(
y

1
)

y

0
2
dt
dt
标准形式改写
dy1
 y2
dt
dy2
2
 y2 (1  y1 )  y1
dt
程序实现
• function dydt = DifferentialCoe(t,y)
• dydt = [y(2);(1-y(1)^2)*y(2)-y(1)];
• 再调用ode23,ode23s,ode23t,ode23tb,ode45等
简单迭代法
设所给的线性方程组为
Mx=g, 其中
 m11 m12
M   m21 m22
m
 n1 mn 2
m1n 
m2 n 
mnn 3 
 x1 
 g1 
x 
g 
2
2


x
,g 
 
 
 
 
 xn 
 gn 
且系数矩阵M为非奇异,g≠0。将方程组改
写成等价形式
x=Ax+f
这种改写的方法很多,例如将M分解为两个
矩阵之差
M=B-C
其中矩阵B可逆,于是方程组成为
Bx=Cx+g
x=B-1 Cx+B-1g
令
A=B-1 C,f=B-1 g
当然选取的B应该便于求逆,如B为单位矩阵
或对角矩阵等。
对任意给定的初始向量x(0),构造迭代公式
x(k+1)=Ax(k)+f,
k=0,1,2,…
这里A称为迭代矩阵。用分量形式可写成
( k 1)
i
x
n
  aij x
j 1
(k )
j
i  1,2, , n
 fi , 
k  0,1,2,
算出迭代公式的解的各次近似值
x(1),x(2), …,x(k),…。这种迭代求解的方法称
为简单迭代法。
特别当M的对角元素均不为零且按绝对值
来说较大时,常取
 m11
B  D  


则
A=D-1 C, f=D-1 g



mnn 
例 用简单迭代法解下列方程组
10 x1  x2  2 x3  7.2

  x1  10 x2  2 x3  8.3
  x  x  5 x  4.2
 1
2
3
解将方程组写成等价形式
 x1  0.1x2  0.2 x3  0.72

 x2  0.1x1  0.2 x3  0.83
 x  0.2 x  0.2 x  0.84
 3
1
2
取初始值x(0)= 0,按迭代公式
 x1( k 1)  0.1x2( k )  0.2 x2( k )  0.72
 ( k 1)
(k )
(k )
x

0.1
x

0.2
x
 2
1
3  0.83
 x ( k 1)  0.2 x ( k )  0.2 x ( k )  0.84
1
2
 3
赛德尔(Seidel)迭代法
从简单迭代法看到,已知x(k)计算x(k+1)时,需要保留
x(k)和x(k+1)两个分量。逐次用前面算出的新分量
来计算下一个分量,这就是赛德尔迭代法的基本思
想。
设方程组的等价形式为
x=Ax+f
将矩阵A分解为
A=B+C
其中
于是
x=Bx+Cx+f
x(k+1)=Bx(k+1)+Cx(k)+f,
k=0,1,2,…
( k 1)
i
x
i 1
  aij x
j 1
( k 1)
j
n
  aij x
j i
(k )
j
i  1,2, , n
 f j,
k  0,1,2,
于是, x=Bx+Cx+f
x(k+1)=Bx(k+1)+Cx(k)+f, k=0,1,2,…
( k 1)
i
x
i 1
  aij x
j 1
( k 1)
j
n
  aij x
j i
(k )
j
i  1,2, , n
 f j,
k  0,1,2,
例 用赛德尔迭代法解方程组
10 x1  x2  2 x3  7.2

  x1  10 x2  2 x3  8.3
  x  x  5 x  4.2
 1
2
3
解 将原方程组写成等价形式并构造赛德尔迭代公式
 x1( k 1)  0.1x2( k )  0.2 x3( k )  0.72
 ( k 1)
( k 1)
(k )
x

0.1
x

0.2
x
 2
1
3  0.83
 x ( k 1)  0.2 x ( k 1)  0.2 x ( k 1)  0.84
1
2
 3
•
•
•
由G-S可得
(I-B)x(k+1)=Cx(k)+f
因为矩阵I-B是非奇异矩阵,则I-B
可逆,所以迭代公式可写成
•
x(k+1)=(I-B)-1
• Cx(k)+(I-B)-1 f,k=0,1,2,…
•
这里迭代矩阵
•
A=(I-B)-1 C
• 由此看出, 赛德尔迭代法实际上是某种简
单迭代法。
• 松驰法
• 赛德尔迭代公式为
•
x(k+1)=Bx(k+1)+Cx(k)+f
• 现令
•
Δx=x(k+1)-x(k)=Bx(k+1)+Cx(k)+f-x(k)
• 于是
•
x(k+1)=x(k)+Δx
x(k+1) 可以看作在向量x(k) 上加修正项Δx而
得到。若在修正项的前面加上一个参数ω,
便得到松驰法的迭代公式
x(k+1)=x(k)+ωΔx
=(1-ω)x(k)+ω(Bx(k+1)+Cx(k)+f)
k=0,1,2,…
或者用分量形式写成
i 1
n
j 1
j i
xi( k 1)  (1   ) xi( k )   (  aij x (jk 1)   aij x (jk )  f i )
i  1, 2, n

 k  0,1, 2,
• 其中ω叫做松驰因子,当ω＞1时叫超松
驰,ω＜1时叫低松驰,ω=1时就是赛德尔
迭代法。松驰因子ω的选取直接影响到松
弛法的收敛性。
• 因为(I-ωB)-1存在,所以还可以改写成
• x(k+1)=(I-ωB)-1((1-ω)I+ωC)x(k)+ω(IωB)-1f 这是简单迭代公式,迭代矩阵为
A=(I-ωB)-1 ((1-ω)I+ωC)
• 实际上,赛德尔迭代法或松弛法都是某种
简单迭代法,它们仅是迭代矩阵A不同。
定理 对于任意初始向量x(0)和常数项f,由
迭代公式
x(k+1)=Ax(k)+f,
k=0,1,2,…
收敛的充要条件是
ρ(A)＜1
定理 若迭代矩阵A的范数‖A‖＜1,则迭代公
式x(k+1)=Ax(k)+f收敛,且有误差估计式
x
(k )
A
x 
x ( k 1)  x ( k )
1 A
*
k
x
(k )
A
x 
x (0)  x (1)
1 A
*
CG：
x0，
r0  b  Ax0 ,
p0  r0 ,
rkT rk
k  T
,
pk Apk
xk 1  xk   k pk ,
rk 1  rk   k Apk ,
rkT1rk 1
k  T ,
rk rk
pk 1  rk 1   k pk ,
关于CG方法请参考Shewchuk的 An introduction to CG
without the agonizing pain.
练习:
课本112页.
a_{ij}=max{i,j} (i~=j), a_{ii}=10*n*i.
右端项b=A*ones(n,1),取n=10.
1.用Jacobi或G-S求解;
2.用CG求解方程组.