《数值分析第3讲(MATLAB简介2).ppt》由会员分享，可在线阅读，更多相关《数值分析第3讲(MATLAB简介2).ppt（36页珍藏版）》请在三一办公上搜索。
1、第三讲 MATLAB基础与程序设计,函数定义说明:function是函数定义关键词，y，n是输出变量，mean是函数名，(x)是输入变量。它们可以是0到若干个,多个时用逗号或空格间隔。输出一个变量时，方括号可以省略，输入变量必须用圆括号。%开头的行为注释行，不执行。第2行为H1行实际上为帮助行。函数帮助文本：在命令行用help mean时显示从H1行到第一个非%行结束。函数体实现函数功能，有返回参数必须在返回前赋值，否则出错或给出不期望的输出值。return语句可以省略，没有它时系统将自动返回调用函数。函数调用可以嵌套，被调用函数称为子函数，调用函数名必须与M文件名相同。函数参数均是局部的，作
2、用范围仅在函数内，除非用函数内外已经用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
3、,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=on
4、es(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
5、)牛顿科斯特法,同上.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)求矩阵的核(零空间
6、)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
7、(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,
8、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,
9、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),t
10、aylor(f,a),taylor(f,x)求泰勒展式dsolve(equ1,equ2,x)求解微分方程,equ可能是初边值条件,符号变量名x缺省时自动用内部符号变量tv=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
11、(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)，上下：
12、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:-
13、实线，-虚线，:点线.，句点线，-.点划线，*星线，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：绘出心脏线（红实线）和正玄四叶玫瑰线（蓝点线）在同一和图上。,解：它们的极坐标方程
14、分别为：,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*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*s
15、in(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
16、 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
17、,.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
18、: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 ttfigure(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,
19、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);endfprintf(fid,);fprintf(fid,%-20.8en,u);fclose(fid);,练习题已知矩阵 11 12 13 14 21 22 23 24 A=31 32 33 34 41 42 43 44A(:,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
20、,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命令完成下列矩阵运算
21、(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,45 已知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
22、),(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
23、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+bsi
24、nx)(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/x13 求下列函数的定积分（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,
25、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)=115 解线性方程组 a1x+b1y+c1z=d1 a2x+b2y+c2z=d2 其中ai,bi,ci,di为常数,x,y,z为变量 a3x+b3y+d3z=d316 求方程 tan(x)+sin(x)=2在区间-2,2上的全部 实数解并用图示。,17 求解以下微分方程.(1)dy/dx=y+sinx y|x=/2=1(2)d2
26、y/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)=218 求解以下微分方程(用符号解法）.(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/d
27、x)+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-oscillat
28、ing 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的区域,首先求y1的最小值与最大值,然后再 求次大值与次小值,21 用MATLAB命令构造一个10
29、1个元素的矢量,它的元素交替为1和-1,使用plot,comb 和stairs 命令分别绘制该矢量的图形22(1)利用prod命令生成与列写N的阶乘(N!),并对N=110进行验证(2)使用for命令代替prod 命令,完成(1)项中的功能(3)使用sum命令求N个数的连续和.并对N=110进行验证(4)使用for命令代替sum命令,完成(3)项中的功能23 绘制分段连续函数图 x2|x|5-5x+50 5x10 5x+50-10 x-5 0|x|10,24 求以下定积分：25 请绘出以下函数的三维面图：f(x,y)=3(1-x)2exp(-x2-(y+1)2)-10(x/5-x3-y5)ex
30、p(-x2-y2)-(1/3)exp(-(x+1)2-y2)26 请练习以下命令 plot semilogx semilogy loglog polar plotyy ezplot ezpolar fplot feather pie bar contour plot3 mesh surf 27 已知微分方程 dy/dx=(x+y)2,初始值为 y(0)=1,在x的取值范围 0,1内求方程的数值解并绘出函数图形。,0,1,0,1+x2,1+x4,dx,28 已知sin0.32=0.314567,sin0.34=0.333487,sin0.36=0.352274,用线性插值、抛物线插值和三次样条插
31、值计算sin0.3367的值并作出误差。给定函数f(x)=,x-5,5，节点xk=-5+k,用三次样条插值求s(xk).30 给出 f(x)=lnx数值表用线性插值及二次插值计算ln0.54的近似值。X 0.4 0.5 0.6 0.7 0.8 lnx-0.916291-0.693147-0.510826-0.357765-0.22314431 给出数据表如下 X 0.25 0.30 0.39 0.45 0.53 Y 0.5000 0.5477 0.6245 0.6708 0.7280 试求三次样条插值s(x)，并满足条件：S(0.25)=1.000,s(0.53)=0.6868;s(0.25)=s(0.53)=0,
