数值分析第二讲MATLAB简介.ppt

数值分析第2讲 MATLAB简介,主讲：王礼广 数理学院 2009年9月,MATLAB简介,MATLAB是“Matrix Laboratery”的缩写，意为“矩阵实验室”，是当今最流行的科学计算软件。在众多的领域，如控制论、时间序列分析、系统仿真、图象信号处理、各种数值模拟等有广泛的用户。由美国Mathwork公司推出至今有10多版了。,MATLAB以矩阵运算作为基本对象，具有丰富的函数库和图形绘制功能，具备交互式和批处理程序的多种用户接口，具备专门个领域的工具箱，这些工具箱均由某领域的专家编写，是其领域最先进的算法。,MATLAB的基本用法,一、数与变量,MATLAB的内部本质上只有两种类型的数据：实数和字符。一个实数占8个字节（64比特），一个字符占用1字节（8比特）。复数由两个实数构成，矩阵元素由字符、实数、复数构成。由矩阵再构成“struct”（结构）数据类型，字符串是以字符为变量的行向量。所以表现形式有：字符、字符串、整数、实数、复数、结构。 characters: char 表有符号字符 8bits, uchar 无符号字符 8bits. integers :short 短整数16 bits, long 长整数 32 bits, floating-point: single- 表单精度浮点数32bits, double表双精度浮点数:64bits,进一步，实数显示的格式又分(以-为例): Shot短（5位有效数字）: -3.1416 shortE短指（底5指3） :-3.1416e+000 shortG（5个有效数字的常规数值）:-3.1416 Long长（16个有效数字的）:-3.14159265358979 longE长指（底16指3）:-3.141592653589793e+000 longG（16个有效数字常规数值）:-3.14159265358979 Hex十六进制（16个有效16进制数）:c00921fb54442d18 Bank银行（银行格式）:-3.14, +符号（数的符号）:- Rational有理数（分数格式）:-355/113,变量命名规则:与C语言类似. 必须以字母开始 中间混用字母与数字 区分大小写 变量名,函数名,文件名均同以上法则 大约不超过32个字符 MATLAB库函数均以小写字母开始,常数(名称及值),ans pi圆周率 eps 机器最小容许差 realmin 最小可用实数 realmax 最大可用实数 flops 浮点运算,用于统计计算量 i, j 默认的虚数单位 inf 无穷大 NaN 不定量,获得帮助, ? 关键词 help关键词 lookfor关键词 doc关键词 demo关键词 tour关键词 进入HELP菜单(比较详细,例子多),常用符号意义,逗号，,空格 ：分隔数据，变量，语句； 分号；：语句换行，抑制屏幕输出，数据换行 点号.:小数点，结构域标志，点乘运算 双点 父目录，续行标志 单引号：标志字符串 百分号%：注释标志 A圈：函数柄（函数指针，函数地址） 美圆号$：系统命令调用前缀 波符号：逻辑非运算 等号=：赋值操作 冒号:：标定矩阵下标范围,运算符,数值运算 + 加 - 减 * 乘 / 右除 左除 幂 .+点加，.-点减， .*点乘， ./点右除， .点左除，.点幂,逻辑运算（还有相应的函数） =等于，=不等于， 大于， =大于或等于， & 逻辑与，|逻辑或，逻辑非， 位运算（还有相应的函数） &位与，|位或，位非 括号运算 （）改变运算优先，函数定义或调用的变量标志部分 构造矩阵（元素间隔用逗号或空格，行间隔号用分号） 构造字符或字符串常数 构造枚举集合，a，b，c，d,MATLAB常用函数名,数学函数（中文意义略） abs，acos，acosh，angle，asin，asinh，atan，atan2，atanh，ceil，conj，cos，cosh，exp，fix，floor，gcd，imag，lcm，log，log10，real，rem，round，sign，sin，sinh，sqrt，tan，tanh， 构造矩阵函数 eye(n)，ones(n)，zeros(n)：产生n阶单位、全1，全0矩阵 eye(m,n)，ones(m.n)，zeros(m.n)：产生mn的单位、全1，全0矩阵 Sparse：构造稀疏矩阵 Size(M)取矩阵M的各维大小(结果为向量) size (M,n) 取矩阵M的第n维大小,矩阵构造函数,随机矩阵 rand(n),rand(m,n),rand(n1,n2,)：构造nn，mn，n1n2的随机矩阵 randn(n),randn(m,n),rand(n1,n2,)：同上,特用正态分布 rand(size(A),rand(seed,a),rand(seed)：同上,带种子a或取机器种子 特殊矩阵 compan(A)取伴随阵， diag(d)取对角阵，gallery取试验矩阵，handmard取Handmard阵 ，hankel取Hankel阵，hilb取希尔伯特阵，invhilb取逆希尔伯特阵，kron取克罗内克张量积，magic取魔方矩阵，pascal取Pascal 阵，toeplotz取Toeplotz矩阵，vander取范得蒙矩阵 ，wilkinson取Wilkinson特性实验矩阵,例1:构造两个可加减乘的矩阵A，B，C,求A+B，A*C，sin(A)+cos(B),以及取A,B,C的子矩阵作同样的运算。,解:A= 1,2,3,4;5 6 7 8; 9 10 11 12; B=ones(3,4) C= pascal(3) D=A+B E=sin(A)+cos(B) (3) F=A*C 错误操作：A*B，A+C,A1=A(1:2,2:4),B1=(2:3,1:3),C1=(1:3,2:3) D1=A1+B1 E1=sin(A1)+cos(B1) F1=A1*C1,例2: 用两种方法求线性方程组Ax=b的解x,并比较它们的差,其中,解:A=2+0.3,sin(pi/4),0,0,0;cos(pi/3),exp(2)+0.22,pi*10,0,0; 0,sinh(1),abs(-4),1e-4,0;0,0,84,2+exp(2),1e-5;0,0,0,0.4,1+log(100) T=0:pi/5:4*pi/5;b=sin(T) x=Ab,x1=inv(A)*b; de=any(x-1);,MATLAB编程入门,MATLAB的编程风格与C语言相同。 常用流程控制语句 if if语句条件 elseif if语句条件 end 终止作用域 for 指定次数的循环 while 不指定次数的循环 break 终止循环 switch 开关语句 case 列出语句 otherwise否则语句 return 返回调用函数 function 函数定义语句,函数定义语句的格式,MATLAB的M文件有两类：脚本文件和函数文件，他们可以用任何的TXT文件编辑器编写。我们将MATLAB环境下（“”提示符下）直接输入的语句放在一个以.m为后缀的文本(TXT)文件中,这样的文件就称为脚本(script)文件。另一类文件是函数文件, 是用于给脚本文件或其他函数文件调用的, 系统调用的规则是:首先在当前目录搜索与被调用函数同名的M文件,然后按照MATLAB的FILE/PREFERENCE菜单项中设置的MATLAB search path 去搜索这个文件名,如果没有找到,则报告出错。函数文件由5部分构成： -函数定义行 - H1行 -函数帮助文件 -函数体 -注释,例如，函数文件mean.m是用于求平均值的,该函数输入数据变量x，输出两个变量，平均值y和数据个数n，它可如下定义： function y，n=mean(x) 函数定义行 % y,n=Mean(x) average or mean value H1行 % for vector x, mean return y, the mean value of x ,函数帮 % and n, the data length of x. 助文件 n,m=size(x); if n=1 函 n=m; 数 end 体 y=sum(x)/n; return ； ,函数定义说明: function是函数定义关键词，y，n是输出变量，mean是函数名，(x)是输入变量。它们可以是0到若干个,多个时用逗号或空格间隔。输出一个变量时，方括号可以省略，输入变量必须用圆括号。 %开头的行为注释行，不执行。第2行为H1行实际上为帮助行。 函数帮助文本：在命令行用help mean时显示从H1行到第一个非%行结束。 函数体实现函数功能，有返回参数必须在返回前赋值，否则出错或给出不期望的输出值。 return语句可以省略，没有它时系统将自动返回调用函数。 函数调用可以嵌套，被调用函数称为子函数，调用函数名必须与M文件名相同。 函数参数均是局部的，作用范围仅在函数内，除非用函数内外已经用global语句定义为全局变量。,常用的一些数值计算函数,插值函数 yi=interp1(x,y,xi,method),一维插值 x,y=n维向量,给定的数据对. xi,yi=m维向量,xi要插值的点向量,yi插值得到的向量 Method=方法字串,nearest，linear，spline，cubic分别指定为最邻近插值，线性内插，三次样条插值，三次插值。 二维，三维插值，类同，仅x为二维和三维数组，y为与x同维的函数值。 yi=interp2(x,y,xi,method)，二维插值。 nearest，linear，spline，cubic yi=interp3(x,y,xi,method)，三维插值。 nearest，linear，spline，cubic,回归及拟合 MATLAB的多项式拟合函数拟合 例:x=0:0.1:1; y=-0.447,1.978,3.28,6.16,7.08,7.34,7.66,9.56,9.48,9.30,11.2; P=polyfit(x,y,2);,拟合过程 t =0 0.3 0.8 1.1 1.6 2.3; y=0.5 0.82 1.14 1.25 1.35 1.40; X1=ones(size(t) t t.2; A=X1y;,例:自建拟合函数拟合。根据已知的数据,分别用下面的两个函数拟合已知数据。,用二次多项式,用指数函数拟合,X2=ones(size(t) exp(-t) t.*exp(-t); B=X2y;,T=0:0.1:2.5; Y1=ones(size(T) T T.2*A; Y2=ones(T) exp(-T) T.*exp(-T)*B; fige(1); Subplot(2,2,1); Plot(T,Y1,-,t,y,o); Title(多项式回归);,Subplot(2,2,2); Plot(T,Y2,-,t,y,o),grid on Title(指数函数回归),数值积分 S=quad(,a,b,tol,trace)辛普生法trace=1展现积分过程图形, trace= 0无图,缺省为无图. S=quad8(fname,a,b,tol,trace)牛顿科斯特法,同上. Z=trapz(X,Y)梯形法,X,Y必须同行数或同维. Sc=cumsum(Y)欧拉法,Sc,Y同维,Sc各列给出Y的积分(求和), 但Sc还需另外乘步长dx(该函数一般用于等步长积分). dlquad(f(x,y),x1,x2,y1,y2)二重积分.,矩阵分析函数,c = cond(X),c = cond(X,p)求矩阵的条件数,p=1,2,inf c = condeig(A),V,D,s = condeig(A)求条件数对应的特征值 d = det(X),求矩阵对应的行列式值 n = norm(A),n = norm(A,p)求矩阵的2(P)范数 Z = null(A),Z = null(A,r)求矩阵的核(零空间) k = rank(A),k = rank(A,tol)求矩阵的秩 b = trace(A)求矩阵的迹 Y = inv(X)求矩阵的逆 /, x=AB等价于求方程Ax=B的解, A/B等价于A*inv(B),d = eig(A),d = eig(A,B),V,D = eig(A),V,D = eig(A,nobalance) V,D = eig(A,B),V,D = eig(A,B,flag)求特征值和特征向量 s = svd(X),U,S,V = svd(X),U,S,V = svd(X,0)求奇异值 L,U = lu(X)求LU分解 Q,R = qr(A) 求QR分解 R = chol(X),R,p = chol(X)求对称正定矩阵的Cholesky分解 expm(X), logm(X), sqrtm(A)求矩阵X的指数,对数,平方函数 T=schur(A),T=schur(A,flag),U,T=schur(A)求矩阵的schur分解 U,V,X,C,S = gsvd(A,B)求广义奇异值分解 P,H = hess(A),H = hess(A)求矩阵的Hessenberg矩阵 B = pinv(A),B = pinv(A,tol)求矩阵的Moore-Penrose广义逆.,线性方程组AX=B求解函数,X=AB 左除法 X=inv(A)*B 高斯法 L,U= lu(X) LU法 L,U = luinc(X,0) 不完全LU分解法 x,flag,relres,iter=cgs(A,b,tol,maxit,M1,M2,x0)共轭梯度法 x,flag,relres,iter=pcg(A,b,tol,maxit,M1,M2,x0)预共轭梯度法 x,flag,relres,iter=bicg(A,b,tol,maxit,M1,M2,x0)双共轭梯度法 x,flag,relres,iter=bicgstab(A,b,tol,maxit,M1,M2,x0)双共稳轭梯度法 x,flag,relres,iter=lsqr(A,b,tol,maxit,M1,M2,x0)最小二乘法 x,flag,relres,iter=gmres(A,b,restart,tol,maxit,M1,M2,x0)广义残余法, x,flag,relres,iter= minres(A,b,tol,maxit,M1,M2,x0)最小残余法 x,flag,relres,iter=symmlq(A,b,tol,maxit,M1,M2,x0)对称最小二乘 x,flag,relres,iter= qmr(A,b,tol,maxit,M1,M2,x0) 拟最小残余法,分析运算 符号变量定义:sym v or syms v1, v2, df=diff(f,n)微分(差分),s必须为符号变量 b=limit(F,x,a) ,F函数,x符号变量,a常数,b极限 int(F,s), int(F,s,a,b)求不定积分或定积分,s为符号变量. taylor(f), taylor(f,n), taylor(f,a), taylor(f,x)求泰勒展式 dsolve(equ1,equ2,x)求解微分方程,equ可能是初边值条件, 符号变量名x缺省时自动用内部符号变量t v=suns(F,s,a)符号变量s替换为常量a代入计算F的值v,solve(eqn1,eqn2,.,eqnN), solve(eqn1,eqn2,.,eqnN,var1,var2,.,varN), solve(eqn1,eqn2,.,eqnN,var1,var2,.varN) 代数方程的符号解,eqn为方程组 SYMSUM(S), SYMSUM(S,v), SYMSUM(S,a,b) and SYMSUM(S,v,a,b)符号求和 F=sym(f(x);创建抽象函数 例: f=sym(f(x);df=(subs(f,x,x+h)-f)/h syms x h df=(subs(f,x,x+h)/h ; %df=(f(x+h)-f(x)/h,代数方程求解,矩阵操作,取矩阵元素：x=A(i,j),矩阵元素赋值：A(i,j)=x 取矩阵子块B=A(i:j,k:l);子块赋值：A(i:j,k:l)=B; 将矩阵元素重新排列成mn阵：reshape(A,m,n) 复制按mn平铺矩阵A：repmat(A,m,n) 按反时针方向旋转k个90度：rot90(A,k) 取A的第k条对角线元素： diag(A),diag(A,k) 矩阵翻转，左右： fliplr(A)，上下：flipud(A) 取矩阵第k条对角线的下三角阵： tril(A) tril(A,k), 取矩阵第k条对角线的上三角阵： triu(A), triu(A,k) 按第dim维合并矩阵： C=cat(dim,A,B,),绘图函数,plot(x1,y1,s1,x2,y2,s2,)绘二维图形，x，y是同规模的横纵坐标数据，s为绘图选择项，可多组重叠绘制。 polar(theta,p,s)以极坐标方式绘二维图形，theta，p是同规模极角与极径数据，s为绘图选项，可多组重叠绘制。 plot3(x,y,z,s)绘制三维图形，x,y,z为同规模三维坐标数据，s为绘图选项 绘图选项s为字串XX，由点形和颜色各取1字符组成： 点X: -实线，-虚线，:点线.，句点线，-.点划线，*星线，o圈线，x叉线，+加号线，s正方线，D菱形线，V下三角线，右三角线，H六角形线，P五角形线 色X: r红，g绿， b蓝， y黄， m洋红， c青色， w白色， k黑色 其他函数:hold on,hold off, axis image,axis(x1,x2,y1,y2),例1：在-上以点距0.1用红绿蓝颜色实线、虚线、点线绘出正弦、余弦、直线y=x/2在同一图上。,解：t=-pi:0.1:pi; plot(t,sin(t),r-,t,cos(t),g-,t,t/2,b:);,例2：绘出心脏线（红实线）和正玄四叶玫瑰线（蓝点线）在同一和图上。,解：它们的极坐标方程分别为：,a=3, t=0:0.01:2*pi; polar(t,a*(1-cos(t),r-); hold on; Polar(t,2*a*sin(2*t),b:);,例3:用点距0.01绘出半径为0.5，中心在原点的下半球面。,解:因为该球面方程为,R=0.5, s=pi/2:0.01:pi,t=0:0.01:2*pi plot3(R*sin(s)*cos(t),R*sin(s)*sin(t),R*cos(s)*ones(size(t),axis image;,例4:用点距0.01绘出螺旋线图形三周。,解：因为螺线方程为,R=1,b=2,t=0:0.01:6*pi; plot3(R*cos(t),R*sin(t),b*t,b-),例5 用欧拉差分法在MATLAB环境编程求解下面的方程，并绘出解图形。,解：程序如下,% matlab programm to compute partitial % derivative equation % with the difference method % Autor:XXX % date: 2003-3-25,time:8.00pm % equation:du/dt-integal(t-s)2*d2u/dx2,0,t=1 % bound: u(0,t)=u(1,t)=0 % initial: u(x,0)=x(x-1) % 0=x=1,0=t=1,clear all;clc %global J ux0 uxN; t0=0;t1=1;x0=0;x1=1; J=input(Please input the x grid numer:n) N=input(Please input the t grid numer:n) dt=(t1-t0)/N; dx=(x1-x0)/J;lambda=dt/dx; t=t0:dt:t1;x=x0:dx: