数值分析第3讲(MATLAB简介2).ppt

第三讲 MATLAB基础与程序设计,函数定义说明: 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),rH); 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:x1; u=zeros(N+1,J+1); u(1,:)=x.*(x-ones(1,J+1); u(:,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,. 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) *(-lambda2/sqrt(dt); 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); IA=inv(A);tic; for i=2:N+1 tt=ones(J-1,1)*dt; u0=u(i-1,2:J); ut=zeros(J-1,1);,for j=2:i-1 u1=u(j,2:J); uf=(B*u1)*lambda2/sqrt(i-j+1)*dt); ut=ut+uf; end u0=u0+tt; ut=ut+u0; u1=IA*ut; u(i,2:J)=u1; end,t0=toc;t1=fix(clock); clear ut0 ux0 uxN A B u0 u1 uf ut tt figure(1); title(The solution surface); mesh(u);grid on; fid=fopen(difds1.txt,w+t); fprintf(fid,%8fn,t0); fprintf(fid,%20dn,t1); fprintf(fid,NxJ=%-5dx%5dn,N,J); fprintf(fid,);,for i=1:N+1 fprintf(fid,); for j=1:J+1 fprintf(fid,%20.8e,u(i,j); end fprintf(fid,n); end fprintf(fid,); fprintf(fid,%-20.8en,u); fclose(fid);,练习题 已知矩阵 11 12 13 14 21 22 23 24 A= 31 32 33 34 41 42 43 44 A(:,1) (2) A(2,:) (3) A(:,2:3) (4) A(2:3,2:3) (5) A(:,1:2:3) (6) A(2:3) (7) A(:) (8) A(:,:) (9) ones(2,2) (10) eye(2) (11) A,ones(2,2);eye(2) (12) diag(A) (13) diag(A,1) (14) diag(A,-1) (15) diag(A,2),2 输入如下矩阵 /3 /6 /2 (2) 求矩阵B1, B1中每一元素为对应矩阵中每一元素的正弦函数 (3) 求矩阵B2, B2中每一元素为对应矩阵中每一元素的余弦函数 (4) 求 B12+B22 (5) 求矩阵的特征值与特征矢量：称特征矢量为，而特征值矩阵为 (6) 求Msin(L)M-1 (7)使用funm命令求矩阵A的正弦函数（结果应该与（）同） (8)求cosA (9) 证明 sin2A+cos2A=I,3 按题目要求用MATLAB命令完成下列矩阵运算 (1) 使用rand命令产生5个2x2随机矩阵A,B,C,D,E (2) 求矩阵F(使用和不使用inv命令两种情况) F=A-1B+C-1(D-1E) 4 手算和上机分别求 A.*B A.B A.B 其中A=1;1;1 B=2,3,4 5 已知A=2 7 6;9 0 10;3 0.5 6; B=8 0.2 0;3 2 5;4 0 7;求 (1)A|B, A (3)求每一个研究生的平均成绩.,7 已知矩阵A=1.2 3 5 0.9;5 1.7 5 6;3 9 0 1;1 2 3 4,求 (1)A的特征多项式 (2)特征多项式中未知数为20 时的值 (3)特征多项式的根 (4)特征多项式的导数 8 已知五个数据点: (1,5.5), (2,43.1), (3,128), (4,290.7), (5,498.4) (1)用三次曲线拟合上述数据点 (2)在同一图形中绘出数据点和拟合曲线 (3)用适当图形表示拟合精度 9 在实验中测得如下6组数据: (0,1.5), (1.5,3.4), (2.8,13), (3.8,36), (4.5,63), (4.9,78) 请用三次曲线拟合以上数据并给出以下结果： (1)三次多项式 的各项系数； (2)将数据点和拟合曲线以最佳效果在同一图中绘出； (3)指出具有最大拟合误差的数据点和最大误差绝对值。,10 求下列函数的极限 （1）lim(x2/sin2(x/3) x 0 （2）lim(tanx-sinx)/sin3x) x 0 （3）lim(sin(a+x)-sin(a-x)/x x 0 （4） limxcos(1/x) x 0 （5） lim(1+mx)n-(1+nx)m)/x2 x 0 （6）lim(1+1/n)(n+5) n （7）lim(1-2/x)x x （8）lim(1+cosx) 3secx x /2 11 求下列函数的积分 （1）x2/sin2(x/3) （2） (tanx-sinx)/sin3x （3）(sin(a+x)-sin(a-x)/x （4）xcos(1/x) （5） (1+mx)n-(1+nx)m)/x2 (6) cos2x (7) sinaxcosbx (8)cosaxsinbx (9) arcsin(x/a) (10)1/(a+bsinx) (11) xarcsin(x/a) (12) x2arcsin(x/a),12 求下列函数的反函数 （1）y=2sin3x （2）y=1+ln(x+2) （3）y=2x/(2x+1) （4） y=(ex+e(-x)/2 （5） y=1/2(arccos(x/2) （6）y=x+1/x 13 求下列函数的定积分 （1）(x+sinx)/(1+cosx) 0,pi/2 （2） ln(1+tanx ) 0,pi/4 （3）1/(1+cos2x ) 0,pi/2 （4）cos5xsinx 0,pi/2 （5） (3x4+3x2+1)/(x2+1) -1,0 (6) x2+1/x4 1,2 (7) tan2 x 0,pi/4 (8) 4cos4x -pi/2,pi/2 (9) 1-sin3x 0,pi (10)1/(11+5x) 3 -2,1 (11) cosxcos2x -pi/2 , pi/2 (12) (x 3 sin 2 x )/(x 4 + 2x 2 + 1) -5,5,14 求解代数方程 (1) ax2+bx+c=0 (2) cos(2x)+sin(x)=1 15 解线性方程组 a1x+b1y+c1z=d1 a2x+b2y+c2z=d2 其中ai, bi, ci, di为常数, x, y, z为变量 a3x+b3y+d3z=d3 16 求方程 tan(x)+sin(x)=2在区间-2, 2上的全部 实数解并用图示。,17 求解以下微分方程. (1) dy/dx=y+sinx y |x= /2=1 (2) d2y/dx2+bdy/dx+a2y=0 dy/dx|x= /a=0 y|x=0=1 (3) dx/dt=2x+y dy/dt=-x+3y x(0)=1 y(0)=2 (4) dy/dx=y+cosx y |x= /2=1 (5) d2y/dx2+bdy/dx+a2y=0 dy/dx|x= /a=1 y|x=0=2 (6) dx/dt=3x+y dy/dt=-2x+5y x(0)=1 y(0)=2 18 求解以下微分方程(用符号解法）. (1) dy/dx=y/x+tany/x (2) (3x2+6xy2) dx +(6x2y+4y3)dy=0 (3) dx/dt=-xcost+(1/2 )sin2t (4)y= (dy/dx) 2 -x (dy/dx)+x 2 /2 (5) dy/dx=(x-y+1)/(x+y-3),19 绘制下列各种函数图形 (1)绘制下列极坐标图形 r=3(1-cos) r=2(1+cos ) r=2(1+sin ) r=cos3 r=exp(4 ) (2)求函数z的三维图形.定义区间与z的函数表达式如下: -5 x 5 , -5 y 5 1 1 z= (x+1)2+(y+1)2+1 (x-1)2+(y-1)2+1,20 用MATLAB命令绘制下列数学函数的图形 (1) 绘制如下函数图形 y(t)=1-2exp(-t)sin(t) (0 t 8) 且在x轴写上“Time”,y轴写上“Amplitude”标号,图形标题为 “Decaying-oscillating Exponential” (2) 绘制如下图形 y(t)=5exp(-0.2t)cos(0.9t-30O)+0.8exp(-2t) (0 t 30) (3)在0 t10区间内绘制下面图形 y(t)=1.23cos(2.83t+240O)+0.625 x(t)=0.625t (4)在0 t20区间内,且在同一图中绘制下述函数图形 y1(t)=2.62exp(-0.25t)cos(2.22t+174O)+0.6 y2(t)=2.62exp(-0.25t)+0.6 y3(t)=0.6 在y值对应的-2到3的区域,首先