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

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

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

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

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

帮助窗口

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

MATLAB的程序构成：

常变量及其命名规则 变量名可以有数字、字母、下划线构成； 变量的首字符必须是字母； 区分变量名的大小写 每个变量名最长只能包含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 4 5 6 7 8 9 10 11 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 …

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指令和必要的注解。需要在工作空间中创建并获取变量，也就是说处理的数据为命令窗口中的数据，没有输入参数，也不会返回参数。 程序运行时只需在工作空间中键入其名称即可。

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

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

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）：分割图形显示窗口 9、semilogx：绘制以x轴为对数坐标（以10为底），y轴为线性坐标的半对数坐标图形。 semilogy：绘制以y轴为对数坐标（以10为底），x轴为线性坐标的半对数坐标图形。 10、了解应用型绘图指令：可用于数值统计分析或离散 数据处理 bar（x,y）；hist（y,x） stairs（x,y）；stem（x,y）

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

曲线拟合 %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 %一次多项是拟合 %已知离散点 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

微分方程 Van der Pol Equation

标准形式改写

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

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

且系数矩阵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称为迭代矩阵。用分量形式可写成

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

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

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

赛德尔(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,…

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,… 例 用赛德尔迭代法解方程组 解 将原方程组写成等价形式并构造赛德尔迭代公式

由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而得到。若在修正项的前面加上一个参数ω,便得到松驰法的迭代公式 =(1-ω)x(k)+ω(Bx(k+1)+Cx(k)+f) 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收敛,且有误差估计式

CG： 关于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. 用Jacobi或G-S求解; 用CG求解方程组.