MATLAB小结、 经典迭代法、CG.

Presentation on theme: "MATLAB小结、 经典迭代法、CG."— Presentation transcript:

1 MATLAB小结、 经典迭代法、CG

2 1 .MATLAB 代表MATrix LABoratory
它的首创者是美国新墨西哥大学计算机系的系主任Cleve Moler博士，他在教授线性代数课程发现其他语言很不方便，篇构思开发了MATLAB。最初采用FORTRAN语言编写，20世纪80年代后出现了MATLAB的第二版，全部采用C语言编写. 1984年Moler博士和一批数学家及软件专家创建了MathWorks公司，专门开发MATLAB。 1993年出现了微机版，到2003年是6.5版

3 2 .一种演草纸式的科学计算语言 3 .MATLAB 是一高性能的技术计算语言. 强大的数值计算和工程运算功能 符号计算功能
强大的科学数据可视化能力 多种工具箱  

4 MATLAB 能干什么？ MATLAB可以进行： MATLAB广泛应用于： 数学计算、算法开发、数据采集 建模、仿真、原型
数据分析、开发和可视化 科学和工程图形应用程序的开发，包括图形用户界面的创建。 MATLAB广泛应用于： 数值计算、图形处理、符号运算、数学建模、系统辨识、小波分析、实时控制、动态仿真等领域。

5 掌握 MATLAB …… MATLAB的构成： MATLAB开发环境：进行应用研究开发的交互式平台
图形化开发：二维、三维图形开发的工具 应用程序接口 (API)：用于与其他预言混编 面向专门领域的工具箱：小波工具箱、神经网络工具箱、信号处理工具箱、图像处理工具箱、模糊逻辑工具箱、优化工具箱、鲁棒控制工具箱等几十个不同应用的工具箱。

6 开发环境 包括：命令窗口、图形窗口、编辑窗口、帮助窗口。

7 命令窗口 可在提示符后输入交互式命令 结果会自动的产生 例如： command (typed at prompt)
MATLAB prompt (>>) and cursor (|) command (typed at prompt) MATLAB output

8 图形窗口 在窗口中输入： Plot([1,2,4,9,16],[1,2,3,4,5]) MATLAB 划出如下图形:

9 编辑窗口 用来创建和修改M-files (MATLAB 脚本)

10 帮助窗口

11 The MATLAB Language MATLAB 语言的特点 Matlab的基本数据单元是不需指定维数的矩阵。
double 是一个双精度浮点数，每个存储的双精度数用64位。 char用于存储字符，每个存储的字符用16位。

12 MATLAB的程序构成：

13 常变量及其命名规则 变量名可以有数字、字母、下划线构成； 变量的首字符必须是字母； 区分变量名的大小写 每个变量名最长只能包含19个字符。

14 Matlab中预定义变量 ans 分配最新计算表达式的值，这个表达式并没有给定一个名字 eps 返回机器精度
realmax 返回计算机能处理的最大浮点数 realmin 返回计算机能处理的最小的非零浮点数 pi , inf 定义为1/0 。当出现被零除时，Matlab就返回inf，并不中断执行而继续计算 NaN 定义为“Not a Number”，这个非数值要么是％类型，要么是inf/inf

15 向量的创建 在matlab的命令窗口键入以下字符 >> a = [1 2 3 4 5 6 9 8 7] a =
希望得到元素从0到20，步距为2的一个向量，只需键入以下命令即可 >> t = [0:2:20] t =

16 矩阵的创建 输入矩阵时每一行元素有分号或者回车键分隔。例如： B = [1 2 3 4;5 6 7 8;9 10 11 12] B =

17 各种运算符

18 语句 变量名列表=表达式 Matlab语言最基本的赋值语句结构为： 注1：整个赋值语句以；结束，则不在屏幕上返回结果，否则立即返回结果。
注2：多个语句可在同一行，用逗号分开。 注3：表达是太长可以用续行符号…

19 函数 Matlab由包括许多标准函数，每个函数都完成某一特定功能的代码组成。
Matlab也允许用户编写自己所需的函数，其扩展名为.m，其中必须以关键字function开头.

20 流程控制 循环语句 for, while 条件转移 if end, if elseif else end 开关语句 switch case
注释语句 % 续行 … 中断语句 break 暂停语句 pause 回显语句 echo on/off

21 1、for循环语句 基本格式 for 循环变量＝起始值：步长：终止值 循环体 end 步长缺省值为1，可以在正实数或负实数范围内任意指定。对于正数，循环变量的值大于终止值时，循环结束；对于负数，循环变量的值小于终止值时，循环结束。循环结构可以嵌套使用。

22 若表达式为真，则执行循环体的内容，执行后再判断表达式是否为真，若不为真，则跳出循环体，向下继续执行。
2、while循环语句 基本格式 while 表达式 循环体 end 若表达式为真，则执行循环体的内容，执行后再判断表达式是否为真，若不为真，则跳出循环体，向下继续执行。 While循环和for循环的区别在于，while循环结构的循环体被执行的次数不是确定的，而for结构中循环体的执行次数是确定的。

23 3、if，else，elseif语句 （1）if 逻辑表达式 执行语句 end （2）if 逻辑表达式 （3） if 逻辑表达式1 执行语句 执行语句1 else elseif 逻辑表达式2 执行语句 执行语句2 end …

24 4、switch语句 switch 表达式（可以是标量或字符串） case 值1 语句1 case 值2 语句2 …. otherwise 语句3 end

25 MATLAB程序的基本组成结构 ％说明 清除命令：清除workspace中的变量和图形（clear,close）
定义变量：包括全局变量的声明及参数值的设定 逐行执行命令：指MATLAB提供的运算指令或工具箱 … … … 提供的专用命令 控制循环 ： 包含for,if then,switch,while等语句 逐行执行命令 … … … end 绘图命令：将运算结果绘制出来 当然更复杂程序还需要调用子程序，或与simulink以及其他应用程序结合起来。

26 MATLAB的程序类型 MATLAB的程序类型有三种，一种是在命令窗口下执行的脚本M文件；另外一种是可以存取的M文件，也即程序文件；最后一种是函数（function）文件。 1、脚本M文件 在命令窗口中输入并执行，它所用的变量都要在工作空间中获取，不需要输入输出参数的调用，退出MATLAB后就释放了。 2、程序M文件 以.m格式进行存取，包含一连串的MATLAB指令和必要的注解。需要在工作空间中创建并获取变量，也就是说处理的数据为命令窗口中的数据，没有输入参数，也不会返回参数。 程序运行时只需在工作空间中键入其名称即可。

27 与在命令窗口中输入命令一样，函数接受输入参数，然后执行并输出结果。用help命令可以显示它的注释说明。 具有标准的基本结构。
3、函数文件 与在命令窗口中输入命令一样，函数接受输入参数，然后执行并输出结果。用help命令可以显示它的注释说明。 具有标准的基本结构。 （1）函数定义行（关键字function） function[out1,out2,..] = filename(in1,in2,..) 输入和输出（返回）的参数个数分别由nargin和nargout两个MATLAB保留的变量来给出。 （2）第一行帮助行，即H1行 以（%）开头，作为lookfor指令搜索的行 （3）函数体说明及有关注解 以（%）开头，用以说明函数的作用及有关内容 （4）函数体语句 函数体内使用的除返回和输入变量这些在function语句中直接引用的变量以外的所有变量都是局部变量，即在该函数返回之后，这些变量会自动在MATLAB的工作空间中清除掉。如果希望这些中间变量成为在整个程序中都起作用的变量，则可以将它们设置为全局变量。

28 Graphics MATLAB提供了丰富的绘图功能 help graph2d可得到所有画二维图形的命令

29 1、基本的绘图命令 plot（x1,y1,option1,x2,y2,option2,…） x1,y1给出的数据分别为x,y轴坐标值，option1为选项参数，以逐点连折线的方式绘制1个二维图形；同时类似地绘制第二个二维图形。 这是plot命令的完全格式，在实际应用中可以根据需要进行简化。比如： plot(x,y)；plot(x,y,option) 选项参数option定义了图形曲线的颜色、线型及标示符号，它由一对单引号括起来。

30 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坐标轴的单位刻度大小调整为一样。

31 6、文字标示 text(x,y,’字符串’) 在图形的指定坐标位置(x,y)处，标示单引号括起来的字符串。 title(‘字符串’) 在所画图形的最上端显示说明该图形标题的字符串。 xlabel(‘字符串’)，ylabel(‘字符串’) 设置x，y坐标轴的名称。 输入特殊的文字需要用反斜杠（\）开头。 7、legend(‘字符串1’,‘字符串2’,…,‘字符串n’) 在屏幕上开启一个小视窗，然后依据绘图命令的先后次序，用对应的字符串区分图形上的线。

32 8、subplot（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）

33 三维的绘图命令 例: 二次型对应的曲面 [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

34 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) ;

35 Examples 绘图实例 函数分析 矩阵运算 线性方程组 曲线拟合 微分方程

36 绘图实例

37 函数分析 fplot('func',[-1 1.5]) %作图 result = func(0) %求函数值
xsolve = fzero('func',3) %求解 Xmin = fminbnd('func',0.5,1) %求最小值

38 矩阵运算 A = [1 2 3 ; ; 7 8 9]; B = [1 2 3 ; 4 5 6]; C = [1 0 1 ; ; 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)

39 线性方程组与特征值 A = [ ; ; ]; b = [ 3.6 ; 2.1 ; -1.4 ]; X = A\b [EigenVectors,EigenValues] = eig(A)

40 曲线拟合 %7次多项是拟合 %已知离散点 x = [1 1.5 3 4 5 6 6.5 7 8];
y = [ ]; %最小二乘拟合 p7 = polyfit(x,y,7); xi = 1:0.25:8; yi = polyval(p7,xi); plot(x,y,'*r',xi,yi); grid on %一次多项是拟合 %已知离散点 x = [ ]; y = [ ]; %最小二乘拟合 p1 = polyfit(x,y,1); y1 = polyval(p1,x); plot(x,y1); hold on plot(x,y,'ro') grid on

41 微分方程 Van der Pol Equation

42 标准形式改写

43 程序实现 function dydt = DifferentialCoe(t,y)
dydt = [y(2);(1-y(1)^2)*y(2)-y(1)]; 再调用ode23,ode23s,ode23t,ode23tb,ode45等

44 简单迭代法 设所给的线性方程组为 Mx=g, 其中

45 且系数矩阵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

46 当然选取的B应该便于求逆,如B为单位矩阵或对角矩阵等。
对任意给定的初始向量x(0),构造迭代公式 x(k+1)=Ax(k)+f, k=0,1,2,… 这里A称为迭代矩阵。用分量形式可写成

47 算出迭代公式的解的各次近似值x(1),x(2), …,x(k),…。这种迭代求解的方法称为简单迭代法。
特别当M的对角元素均不为零且按绝对值来说较大时,常取 则 A=D-1 C, f=D-1 g

48 例 用简单迭代法解下列方程组 解将方程组写成等价形式

49 取初始值x(0)= 0,按迭代公式

50 赛德尔(Seidel)迭代法 从简单迭代法看到,已知x(k)计算x(k+1)时,需要保留x(k)和x(k+1)两个分量。逐次用前面算出的新分量来计算下一个分量,这就是赛德尔迭代法的基本思想。 设方程组的等价形式为 x=Ax+f 将矩阵A分解为 A=B+C 其中

51 于是 x=Bx+Cx+f x(k+1)=Bx(k+1)+Cx(k)+f, k=0,1,2,…

52 x(k+1)=Bx(k+1)+Cx(k)+f, k=0,1,2,…
于是, x=Bx+Cx+f x(k+1)=Bx(k+1)+Cx(k)+f, k=0,1,2,… 例 用赛德尔迭代法解方程组 解 将原方程组写成等价形式并构造赛德尔迭代公式

53

54 由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 由此看出, 赛德尔迭代法实际上是某种简单迭代法。

55 松驰法 赛德尔迭代公式为 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

56 x(k+1)可以看作在向量x(k)上加修正项Δx而得到。若在修正项的前面加上一个参数ω,便得到松驰法的迭代公式
=(1-ω)x(k)+ω(Bx(k+1)+Cx(k)+f) k=0,1,2,… 或者用分量形式写成

57 其中ω叫做松驰因子,当ω＞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不同。

58 定理 对于任意初始向量x(0)和常数项f,由迭代公式
x(k+1)=Ax(k)+f, k=0,1,2,…  收敛的充要条件是 ρ(A)＜1 定理 若迭代矩阵A的范数‖A‖＜1,则迭代公式x(k+1)=Ax(k)+f收敛,且有误差估计式

59 CG： 关于CG方法请参考Shewchuk的 An introduction to CG without the agonizing pain.

60 练习: 课本112页. a_{ij}=max{i,j} (i~=j), a_{ii}=10*n*i. 右端项b=A*ones(n,1),取n=10. 用Jacobi或G-S求解; 用CG求解方程组.