数值分析第二讲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.1416shortE短指（底5指3）:-3.1416e+000shortG（5个有效数字的常规数值）:-3.1416LonglongEe+000longG（16Hex十六进制（16个有效16进制数）:c00921fb54442d18Bank银行（银行格式）:-3.14,+符号（数的符号）:-Rational有理数（分数格式）:-355/113变量命名规则:与C语言类似.必须以字母开始中间混用字母与数字区分大小写变量名,函数名,文件名均同以上法则大约不超过32个字符MATLAB库函数均以小写字母开始常数(名称及值)anspi圆周率eps 机器最小容许差realmin 最小可用实数realmax 最大可用实数flops 浮点运算,用于统计计算量i,j 默认的虚数单位inf 无穷大NaN 不定量获得帮助?关键词help关键词 lookfor关键词doc关键词demo关键词tour关键词进入HELP菜单(比较详细,例子多)常用符号意义逗号，,空格：分隔数据，变量，语句；分号；：语句换行，抑制屏幕输出，数据换行点号.:小数点，结构域标志，点乘运算双点.父目录，续行标志单引号：标志字符串百分号%：注释标志A圈：函数柄（函数指针，函数地址）美圆号$：系统命令调用前缀波符号：逻辑非运算等号=：赋值操作冒号:：标定矩阵下标范围运算符数值运算+加-减*乘/右除 左除 幂.+点加，.-点减，.*点乘，./点右除，.点左除，.点幂逻辑运算（还有相应的函数）=等于，=不等于，大于，=大于或等于，&逻辑与，|逻辑或，逻辑非，位运算（还有相应的函数）&位与，|位或，位非括号运算（）改变运算优先，函数定义或调用的变量标志部分构造矩阵（元素间隔用逗号或空格，行间隔号用分号）构造字符或字符串常数构造枚举集合，a，b，c，dMATLAB常用函数名数学函数（中文意义略）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+CA1=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，cubicyi=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 onTitle(指数函数回归)数值积分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,infc=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缺省时自动用内部符号变量tv=suns(F,s,a)符号变量s替换为常量a代入计算F的值vsolve(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*piplot3(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%matlab programm to compute partitial%derivative equation%derivative equation%with the difference method%with the difference method%Autor:XXX%Autor:XXX%date:2003-3-25,time:8.00pm%date:2003-3-25,time:8.00pm%equation:du/dt-integal(t-s)2*d2u/dx2,0,t=1%equation:du/dt-integal(t-s)2*d2u/dx2,0,t=1%bound:u(0,t)=u(1,t)=0%bound:u(0,t)=u(1,t)=0%initial:u(x,0)=x(x-1)%initial:u(x,0)=x(x-1)%0=x=1,0=t=1%0=x=1,0=t=1clear all;clcclear all;clc%global J ux0 uxN;%global J ux0 uxN;t0=0;t1=1;x0=0;x1=1;t0=0;t1=1;x0=0;x1=1;J=input(Please input the x grid numer:n)J=input(Please input the x grid numer:n)N=input(Please input the t grid numer:n)N=input(Please input the t grid numer:n)dt=(t1-t0)/N;dx=(x1-x0)/J;lambda=dt/dx;dt=(t1-t0)/N;dx=(x1-x0)/J;lambda=dt/dx;t=t0:dt:t1;x=x0:dx:x1;t=t0:dt:t1;x=x0:dx:x1;u=zeros(N+1,J+1);u=zeros(N+1,J+1);u(1,:)=x.*(x-ones(1,J+1);u(1,:)=x.*(x-ones(1,J+1);u(:,1)=zeros(N+1,1);u(:,1)=zeros(N+1,1);u(:,J+1)=zeros(N+1,1);u(:,J+1)=zeros(N+1,1);A=eye(J-1)*(1+2*lambda2/sqrt(dt)+sparse(2:J-1,.A=eye(J-1)*(1+2*lambda2/sqrt(dt)+sparse(2:J-1,.1:J-2,ones(1,J-2),J-1,J-1)*(-lambda2/sqrt(dt)1:J-2,ones(1,J-2),J-1,J-1)*(-lambda2/sqrt(dt)+(sparse(2:J-1,1:J-2,ones(1,J-2),J-1,J-1)+(sparse(2:J-1,1:J-2,ones(1,J-2),J-1,J-1)*(-lambda2/sqrt(dt);*(-lambda2/sqrt(dt);B=eye(J-1)*(-2)+sparse(2:J-1,1:J-2,ones(1,J-2),B=eye(J-1)*(-2)+sparse(2:J-1,1:J-2,ones(1,J-2),J-1,J-1)+(sparse(2:J-1,1:J-2,ones(1,J-2),J-1,J-1);J-1,J-1)+(sparse(2:J-1,1:J-2,ones(1,J-2),J-1,J-1);IA=inv(A);tic;IA=inv(A);tic;for i=2:N+1for i=2:N+1 tt=ones(J-1,1)*dt;tt=ones(J-1,1)*dt;u0=u(i-1,2:J);u0=u(i-1,2:J);ut=zeros(J-1,1);ut=zeros(J-1,1);for j=2:i-1 for j=2:i-1 u1=u(j,2:J);u1=u(j,2:J);uf=(B*u1)*lambda2/sqrt(i-j+1)*dt);uf=(B*u1)*lambda2/sqrt(i-j+1)*dt);ut=ut+uf;ut=ut+uf;end end u0=u0+tt;u0=u0+tt;ut=ut+u0;ut=ut+u0;u1=IA*ut;u1=IA*ut;u(i,2:J)=u1;u(i,2:J)=u1;endendt0=toc;t1=fix(clock);t0=toc;t1=fix(clock);clear ut0 ux0 uxN A B u0 u1 uf ut ttclear ut0 ux0 uxN A B u0 u1 uf ut ttfigure(1);figure(1);title(The solution surface);title(The solution surface);mesh(u);grid on;mesh(u);grid on;fid=fopen(difds1.txt,w+t);fid=fopen(difds1.txt,w+t);fprintf(fid,%8fn,t0);fprintf(fid,%8fn,t0);fprintf(fid,%20dn,t1);fprintf(fid,%20dn,t1);fprintf(fid,NxJ=%-5dx%5dn,N,J);fprintf(fid,NxJ=%-5dx%5dn,N,J);fprintf(fid,);fprintf(fid,);for i=1:N+1for i=1:N+1 fprintf(fid,);fprintf(fid,);for j=1:J+1 for j=1:J+1 fprintf(fid,%20.8e,u(i,j);fprintf(fid,%20.8e,u(i,j);end end fprintf(fid,n);fprintf(fid,n);endendfprintf(fid,);fprintf(fid,);fprintf(fid,%-20.8en,u);fprintf(fid,%-20.8en,u);fclose(fid);fclose(fid);