matlab与科学计算(简单易学老师的ppt很不错).ppt

《matlab与科学计算(简单易学老师的ppt很不错).ppt》由会员分享，可在线阅读，更多相关《matlab与科学计算(简单易学老师的ppt很不错).ppt（90页珍藏版）》请在三一办公上搜索。

1、matlab科学工程计算,同济大学数学系 陈雄达,2,teapotdemo,3,目录,matlab入门matlab基础编程matlab作图matlab工具箱matlab上机操作matlab综合训练matlab命令速查,4,matlab入门,matlab是什么matlab不是什么matlab集成环境matlab ABC,5,matlab入门,没有matlab就没有乐趣。M.N.Nachtigal S.C.Reddy L.N.Trefethen 关于迭代法的Copper Mountain 论文集 1990,6,matlab是什么,一个可视化的计算程序,广泛使用于从个人计算机到超级计算机范围内的各种

2、计算机上包括命令控制、可编程，上百个预先定义命令和函数有许多强有力的命令,能完成大量的高级矩阵处理强有力的二维、三维图形工具能与其他程序一起使用25个(不断增加中)不同的工具箱应用于特殊的应用领域工业研究与开发的有力工具数学教学,尤其线代,数值分析,科学计算方面的教研工具电子学,控制理论,物理学等工程科学方面的教研工具经济学,化学和生物学等有计算问题的所有领域中的教学与研究名字取自矩阵实验室(matrix laboratory),7,matlab不是什么,不是万能的解决工具不是最高性能的编程语言受计算条件限制,不能解决超大型实际问题不能解决工具箱之外的问题种类-需要编写接口、算法甚至工具箱,8

3、,集成环境,窗口系统View Desktop Layout Five Panel历史命令/变量和文件管理/命令窗口菜单系统File/Edit/View/Web/Window/Help按钮,9,matlab ABC,提示符 注释符%续行.123+45-67+8 9x=3*sin(pi/4)2 向量(数组)v=1 3 5 2 4 6;矩阵(二维数组)A=1 3;5,2 4 6;,10,matlab基础编程,数据结构冒号(:)矩阵操作入门标识符初等函数(elfun)初等函数(exp)函数fix,round,ceil,floor逻辑判断关系运算,运算的级别matlab帮助结构化编程脚本文件函数文件,1

4、1,数据结构,最基本的数据结构-矩阵数和向量看成为特殊的矩阵矩阵以 为定界符，与多维数组等同字符串看成为每个元素都是单个字符的向量，也可以有字符矩阵高维数组细胞(cell)结构体(struct),12,冒号(:),冒号 a:s:b 从a开始,步长为s,界为b例如 1:2:10 1 3 5 7 9 7:-2:0 7 5 3 1 1:6 1 2 3 4 5 6 6:2(空矩阵),13,标识符,文件名、变量名、函数名标识符的规定最长不超过19个字符为配合matlab的风格，采用见名知意的小写短名称以字母开头，包含数字、大小写字母、下划线变量命名规则常用从简，专用从繁矩阵大写 A,B,C向量小写 u,

5、v,w,x,y,z函数小写 f,g,h,fun,f1,f2常量小写 alpha,beta,a,b,c,14,内部变量,pi 圆周率=3.1415926i,j 虚根inf 无限大eps 2.2204460e-16NaN 不定型(Not a Number)内部变量的运算规则inf 参与的运算NaN 参与的运算,15,点运算(向量运算),a=1 2 3 4;b=0 1 2 3;a+bans=1 3 5 7 b./aans=0 0.5000 0.6667 0.7500 b.*aans=0 2 6 12 a.bans=1 2 9 64,16,矩阵操作入门,矩阵的抽取 A=8 1 6;3 5 7;4 9

6、2;A(2,3)得到 7 A(2,3 2 1)得到 7 5 3 A(3 2,2:-1:1)得到 9 4;5 3 A(2,:)得到 3 5 7 A(:,end)得到 6;7;2,17,赋值语句,A(2)=A(n-1,:)=A(2+3,2*3)=x=A(sqrt(4),x-3)=A(:)=v(end)=,A(2)+2=x2+x=A(sqrt(4),x)-3=A=inv(A)=(x=0)=,左值表达式:指向某一存储空间的表达式，放在赋值语句左边。,18,初等函数(elfun),三角函数(sin,sinh,asin,asinh)指数对数(exp,log,log10,sqrt)复变函数(abs,conj

7、,real,imag)取整取余(floor,ceil,round,mod)全部都具有向量功能:输入向量,返回对应函数值组成的向量,19,初等函数(exp),x=0:0.2:1x=0 0.2 0.4 0.6 0.8 1.0 y=exp(x)y=1.00 1.22 1.49 1.82 2.22 2.71,向量功能,20,函数fix,round,ceil,floor,21,简单矩阵操作,rand(n),eye(n),zeros(n)生成 n 阶随机、单位、零矩阵,rand(2)ans=0.9501 0.6068 0.2311 0.4860 rand(3,2)ans=0.8913 0.0185 0.7

8、621 0.8214 0.4564 0.4447,eye(2,3)ans=1 0 0 0 1 0 zeros(3)ans=0 0 0 0 0 0 0 0 0,22,diag,diag(1 2 3;4 5 6 7 8 9)ans=1 5 9 diag(1 2 3;4 5 6 7 8 9,1)ans=2 6,diag(1 5 9)ans=1 0 0 0 5 0 0 0 9 diag(4 8,-1)ans=0 0 0 4 0 0 0 8 0,23,triu/tril,A=magic(3)A=8 1 6 3 5 7 4 9 2 U=triu(A,-1)U=8 1 6 3 5 7 0 9 2,L1=tr

9、il(A,-1)L1=0 0 0 3 0 0 4 9 0 L2=tril(A)L2=8 0 0 3 5 0 4 9 2,24,字符运算,char(4,1,8,8,11+100)ans=hello abs(matlab)ans=109 97 116 108 97 98 eval(exp(-pi*i)ans=-1+0i,feval(sin,pi/2)ans=1 fun=cos;feval(fun,pi/2)ans=0,字符串,25,字符运算,for k=1:9,st=fun.num2str(k);x(k)=.feval(st,x(k);endfunction y=fun1(x)function y

10、=fun2(x),str2num(2*pi)ans=6.2832str2num(cos(2*pi)ans=1eval(cos(2*pi)ans=1,26,逻辑判断,真值 逻辑真(true,1)或假(false,0)逻辑运算,=,=,=,=例如又如,pi=sqrt(10)ans=0,1:5=pians=0 0 0 1 1,向量功能,27,关系运算,关系运算有：与(and)and(a,b)a&b或(or)or(a,b)a|b非(not)not(a)a异或(xor)xor(a,b),28,两个特殊关系运算(any,all),仅当输入向量 v 每个分量为真(非零值),all(v)返回1；否则返回0。仅

11、当输入向量 v 每个分量为假(零),any(v)返回0；否则返回1。any(0 1 2 3)ans=1 all(0 1 2 3)ans=0,29,运算的级别,函数、括号算术运算逻辑运算关系运算每一个运算级别内仍有优先级的高低,30,matlab帮助,helphelp 命令名/工具箱名help helpdocdoc eigwhat列出当前matlab文件who/whos列出当前变量which寻找matlab文件或命令所在目录lookfor比较 lookfor curve/help curve,31,结构化编程,分支结构(if)for循环while循环嵌套结构程序规范,32,分支结构(if),if

12、 判断条件1 执行语句1elseif 判断条件2 执行语句2elseif 判断条件3 执行语句3else 执行语句n+1end,if temp 10,disp(cold!);elseif temp 20,disp(cool);elseif temp 27,disp(warm);elseif temp 31,disp(hot);else disp(too hot);end,33,for循环,for 循环变量=循环值 循环操作end,s=0;for k=1:100,s=s+1/k2;end,for k=hello world,fprintf(%c,k);pause(0.1);end,for k=m

13、agic(3),k+1end,向量功能,34,while循环,while 循环条件 循环操作end while 1,if 循环条件,break;end 循环操作end,n=input(n=);while n=1,if mod(n,2)=1,n=n*3+1;else n=n/2;endend,Collatz猜想,35,嵌套结构,for k=1:9,for j=1:9,fprintf(%d*%d=%3d,k,j,k*j);end fprintf(n);end,九九乘法表,36,嵌套结构,n=input(n=0:);flag=0;while flag,nn=n;rn=0;while nn=0,d=m

14、od(nn,10);rn=rn*10+d;nn=floor(nn/10);end disp(rn)if rn=n,flag=1;else n=rn+n;endend,输入n 设置迭代指标为0 当迭代不成功时(=0),保留当前n 计算n的逆序为rn 显示rn 重新设置迭代指标 或 n=rn+n end%迭代不成功,伪代码,源代码,对一个整数n,把它反过来读的数称为逆序数。如132的逆序数为231。写一个算法，输入n,若它与自己的逆序数相等,则停止；否则把它与自己的逆序数相加，重新判断。,37,程序规范,程序按照结构缩进(matlab会自行处理)按照逻辑关系分段添加适当的注释(5%-30%)规范所

15、有的标识符(文件名、变量名)适当增加空行、空格以保证可读性办法:参考Matlab内部程序风格,38,脚本文件,n=input(years:);for k=n,if mod(k,400)=0,disp(num2str(k)is a leap year);elseif mod(k,100)=0,disp(num2str(k)isnt a leap year);elseif mod(k,4)=0,disp(num2str(k)is a leap year);else disp(num2str(k)isnt a leap year);endend,判别闰年不用这么烦吧？,n=input(years:)

16、;for k=n,if mod(k,400)=0|.(mod(k,4)=0 endend,好像简单多了.,39,函数文件,function p,t=fun(n)%n:input integer%p:sum of factors%t:number of factors p=0;t=0;for k=1:n,if mod(n,k)=0,p=p+k;t=t+1;end end,function p,t=fun(n)%n:input integer%p:sum of factors%t:number of factors k=1:n;z=mod(n,k);t=sum(z=0);p=sum(k(z=0);

17、,保存为fun.m,40,函数的调用,z,p=fun(15)z=24p=4 x=fun(6)x=12 fun(28)ans=56,p,z=funInput argument n is undefined.p,z=fun(20,18)Too many input arguments.p,t,b=fun(15)Too many output arguments.,function p,t=fun(n),41,函数的相互调用,寻找亲和数对,即这样两个数,它们的真因子之和等于对方。,function X=amic(n)s=1;for k=2:n,if k=fun(fun(k),X(s,1 2)=k,f

18、un(k);s=s+1;end end,function p=fun(n)k=1:n-1;z=mod(n,k);p=sum(k(z=0);,调用,子函数,42,递归调用,函数调用自己称为递归因此必须有停止自调用的时刻。,例：Frayer级数,即0,1上分母不超过n的既约真分数从小到大排列。,每一次把新出现的分数,放在两个已有分数中,这两个分数的分子和分母各自的和就是新分数的分子和分母。,43,递归调用,function p,q=frayer(n)if n1,p,q=frayer(n-1);ind=find(q(1:end-1)+q(2:end)=n);for k=ind(end:-1:1),q

19、=q(1:k)n q(k+1:end);p=p(1:k)p(k)+p(k+1)p(k+1:end);end elseif n=1,p=0,1;q=1,1;end,44,递归调用,function v=fact(n)if n=0|n=1,v=1;else v=n*fact(n-1);end,v=prod(1:n);,阶乘,45,算法的复杂度,s=0;for k=1:n,p=1;for j=1:k,p=p*j;end s=s+1/(p+1);end,op=0;for k=1:n,op=op+1;for j=1:k,op=op+1;end op=op+3;end,op=4n+n(n+1)/2,46,

20、算法的复杂度,s=0;p=1;for k=1:n,p=p*k;s=s+1/(p+1);end,op=2;for k=1:n,op=op+1;op=op+3;end,op=4n+2,47,matlab作图,二维作图三维作图其它作图,48,二维作图(plot),plot(x,y,s)例如 x=linspace(0,2*pi);y=sin(x);plot(x,y,r-);,49,二维作图(plot),plot(x,y,s)-s的规定,50,二维作图(plot),xlabel(x);ylabel(y);title(y=sin(x);,51,二维作图(polar),theta=linspace(0,6*

21、pi);rho=theta;polar(theta,rho,r-);,52,三维作图(plot3),plot3(x,y,z,s)t=0:pi/50:10*pi;plot3(sin(t),cos(t),t);,53,三维作图(comet3),t=0:pi/500:20*pi;comet3(sin(t),cos(t),t);,54,三维作图(meshgrid),X,Y=meshgrid(x,y)例如 X，Y=meshgrid(1:3,4:7)则 X=1 2 3 且 Y=4 4 4 1 2 3 5 5 5 1 2 3 6 6 6 1 2 3 7 7 7,55,三维作图(mesh),x=linspac

22、e(-5,5);X,Y=meshgrid(x);Z=X.2 Y.2;mesh(X,Y,Z);view(-60,28);,如果X,Y由meshgrid(x,y)生成，可写 mesh(x,y,Z)surf有同样的格式,向量功能,56,其它作图(bar3),A=82 82 63 81 82 81 84 61 84 85 80 79 65 84 86;bar3(A,grouped);,57,其它作图(pie3),pie3(2 4 3 5,0 1 1 0,.North,South,East,West),58,surf,contour,例:山区地貌-在某山区测得一些地点的高程如下表。平面区域为 1200=

23、x=4000,1200=y=3600)。试作出该山区的地貌图和等高线图。,59,surf,contour,A=1130 1250 1280 1230 1040 900 500 700 1320 1450 1420 1400 1300 700 900 850 1390 1500 1500 1400 900 1100 1060 950 1500 1200 1100 1350 1450 1200 1150 1010 1500 1200 1100 1550 1600 1550 1380 1070 1500 1550 1600 1550 1600 1600 1600 1550 1480 1500 155

24、0 1510 1430 1300 1200 980;x=1200:400:4000;y=1200:400:3600;surf(x,y,A);contour(x,y,A,100);,60,matlab工具箱(toolbox),矩阵计算(help matfun)最优化工具箱(help optim)微分方程工具箱(help ode45,help pde)多项式计算(help polyfun)统计工具箱(help stats),61,矩阵计算(inv),inv(A)求矩阵A的逆例 A=magic(3)A=8 1 6 3 5 7 4 9 2 inv(A)ans=0.1472-0.1444 0.0639-

25、0.0611 0.0222 0.1056-0.0194 0.1889-0.1028,62,矩阵计算(),A=hilb(3)A=1.0000 0.5000 0.3333 0.5000 0.3333 0.2500 0.3333 0.2500 0.2000 b=ones(3,1);%3行1列全是1的向量 x=inv(A)*bx=3.0000-24.0000 30.0000,A=hilb(3)A=1.0000 0.5000 0.3333 0.5000 0.3333 0.2500 0.3333 0.2500 0.2000 b=ones(3,1);%3行1列全是1的向量 x=A bx=3.0000-24.

26、0000 30.0000,63,矩阵的秩(rank),A=magic(4)A=16 2 3 13 5 11 10 8 9 7 6 12 4 14 15 1 rank(A)ans=3,64,矩阵的行列式(det),A=pascal(4)A=1 1 1 1 1 2 3 4 1 3 6 10 1 4 10 20 det(A)ans=1,65,矩阵的迹(trace),A=1:4*1:4A=1 2 3 4 2 4 6 8 3 6 9 12 4 8 12 16 trace(A)ans=30,66,矩阵的Cholesky分解(chol),A=pascal(4)A=1 1 1 1 1 2 3 4 1 3 6

27、10 1 4 10 20 B=chol(A)B=1 1 1 1 0 1 2 3 0 0 1 3 0 0 0 1,67,矩阵的特征值(eig),A=gallery(3)A=-149-50-154 537 180 546-27-9-25 eig(A)ans=1.00000000001122 1.99999999999162 2.99999999999700,V,D=eig(A)V=0.3162-0.4041-0.1391-0.9487 0.9091 0.9740-0.0000 0.1010-0.1789D=1.0000 0 0 0 2.0000 0 0 0 3.0000,68,最优化工具箱(lin

28、prog),什么是线性规划？,f=1 1;A=-2-1;-5-8;b=-10-40;Aeq=;beq=;lb=0 4;ub=7 inf;linprog(f,A,b,Aeq,beq,lb,ub)Optimization terminated successfully.ans=3.0000 4.0000,69,最优化工具箱(fminu),x,opts=fminu(f,x0,options,g,p1,p2,),最优解向量,求解情况,工具箱函数,目标函数字符串,初值,算法参数,目标函数梯度字符串,目标函数其他参数,helpfoptions,70,最优化工具箱(fminu),function v=f(x

29、)v=100*(x(2)-x(1)2)2+.(x(1)-1)2;x=fminu(f,0 0)x=1.0000 1.0000,71,最优化工具箱(fzero),x,fv=fzero(f,x0,options,p1,p2,),零点向量,函数值,工具箱函数,目标函数字符串,初值,算法参数,目标函数其他参数,helpfoptions,72,最优化工具箱(fzero),f=inline(x+exp(x)-3);x,fv=fzero(f,2)x=0.79205996843068fv=0只能求单变量函数的根！,73,最优化工具箱(fmincon),x=fmincon(fun,x0,.1 1 1,8,.1 2

30、 0,6,.0 0 0,.cons.);function v=fun(x)v=-3*x(1)-x(2)2+x(3);function c,ceq=cons(x)c(1)=x*x-16;c(2)=x(2)*(x(1)+x(3)-12;ceq=;,最优解6 0 0,74,微分方程工具箱(ode45),t,y=ode45(odefun,tspan,y0,options,p1,p2.),求解函数,工具箱函数,算法参数,helpodeset,求解区间,初值,解的向量可用plot(x,y)画出解,目标函数其他参数,75,微分方程工具箱(ode45),function dydx=f(x,y)dydx=y x

31、+1;x,y=.ode45(f,0,2,1);plot(x,y,r-);,76,多项式计算,多项式表示：向量,p=1 0 2 表示,多项式求值 polyval(p,1:5)多项式求根 r=roots(p)用根生成多项式 p=poly(r),多项式相乘 p=1 0 2;q=2 1;s=conv(p,q);多项式相除 q,r=deconv(s,p);多项式求导 d=polyder(p),77,多项式计算,多项式拟合(polyfit)x=linspace(0,pi,8)x=0 0.4488 0.8976 1.3464 1.7952 2.2440 2.6928 3.1416 y=sin(x)y=0 0

32、.4339 0.7818 0.9749 0.9749 0.7818 0.4339 0.0000 p=polyfit(x,y,2)p=-0.4011 1.2601-0.0179 xx=linspace(0,pi,200);yy=polyval(p,xx);plot(xx,yy,r-,x,y,b-,linewidth,2),78,插值(interp1),yi=interp1(x，y，xi，method),nearest:最邻近插值linear:线性插值spline:三次样条插值 cubic:立方插值缺省时:分段线性插值,79,插值(interp2),z=interp2(x0,y0,z0,x,y,m

33、ethod),nearest 最邻近插值linear 双线性插值cubic 双三次插值缺省时,双线性插值,80,统计工具箱,参数估计betafit,binofit,mle,概率密度函数betapdf,binopdf,hygepdf,各种分布的均值、方差betastat,binostat,chi2stat,线性、非线性的各种回归模型anova1,anova2,polyfit,regress,lsqnonneg,nlinfit,nlintool,81,matlab上机操作,遵守机房各项规定保持机房卫生整洁matlab系统若有问题，请重新开机，在选择进入系统前按Ctrl-R,82,matlab练习题

34、,求下面表达式的值 其中,83,matlab练习题,已知矩阵A,B如下,计算并比较 A*B,B*A,A.*B,B.*A,A/B,AB,A./B,A.B,84,matlab练习题,已知矩阵A,B如下,给出这两个矩阵的(1)逆(2)行列式(3)迹(4)特征值和特征向量(5)最简阶梯形矩阵,85,matlab练习题,已知矩阵A及多项式如下,给出多项式在矩阵上的值；比较结果矩阵和原矩阵的特征值和特征向量。,86,matlab综合训练,矩阵:写一个程序，对某一输入的对称矩阵A进行分解,A=LDL,其中L为单位下三角矩阵，D为对角矩阵。参考数值计算中的Cholesky分解。画图:给出某一曲面方程及一垂直于

35、XOY平面的平面方程，利用平面把曲面割成两片，画出这两片，中间的缝隙要逐渐变大。算法:输入平面上的至少四个点，挑选当中的部分点，并把它们连成一个多边形，要求面积尽量大，给出其面积之值。,87,matlab命令速查,:.*./.abs all and any asin asinh bar3 ceil cell char chol comet3 conj contour conv deconv det diag disp doc else elseif end eig exp eval eye false feval fix floor fmincon fminu for fprintf func

36、tion fzero help hilb i if imag inf input interp1 interp2 inv j linprog linspace log log10 lookfor magic mesh meshgrid NaN not num2str ode45 ones or pascal pi pie3 plot plot3 polar poly polyder polyfit polyval prod rand rank real roots round sin sinh sqrt str2num struct surf title trace tril triu true view what which while who whos xlabel xor ylabel zeros,88,作业:三维画图,在某海域测得一些点(x,y)处的水深z由下表给出，船的吃水深度为5英尺，在矩形区域（75，200）*（-50，150）里的哪些地方船要避免进入。,89,作业:三维画图,90,THANK YOU,