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

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是什么,一个可视化的计算程序, 广泛使用于从个人计算机到超级计算机范围内的各种计算机上 包括命令控制、可编程，上百个预先定义命令和函数 有许多强有力的命令, 能完成大量的高级矩阵处理 强有力的二维、三维图形工具 能与其他程序一起使用 25个(不断增加中)不同的工具箱应用于特殊的应用领域 工业研究与开发的有力工具 数学教学, 尤其线代, 数值分析, 科学计算方面的教研工具 电子学, 控制理论, 物理学等工程科学方面的教研工具 经济学, 化学和生物学等有计算问题的所有领域中的教学与研究 名字取自矩阵实验室(matrix laboratory),7,matlab不是什么,不是万能的解决工具 不是最高性能的编程语言 受计算条件限制, 不能解决超大型实际问题 不能解决工具箱之外的问题种类 - 需要编写接口、算法甚至工具箱,8,集成环境,窗口系统 View Desktop Layout Five Panel 历史命令 / 变量和文件管理 / 命令窗口 菜单系统 File / Edit / View / Web / Window / Help 按钮,9,matlab ABC,提示符 注释符 % 续行 . 123 + 45 - 67 + 8 9 x = 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帮助 结构化编程 脚本文件 函数文件,11,数据结构,最基本的数据结构- 矩阵 数和向量看成为特殊的矩阵 矩阵以 为定界符，与多维数组等同 字符串看成为每个元素都是单个字符的向量，也可以有字符矩阵 高维数组 细胞 (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, v, w, x, y, z 函数小写 f, g, h, fun, f1, f2 常量小写 alpha, beta, a, b, c,14,内部变量,pi 圆周率=3.1415926 i,j 虚根 inf 无限大 eps 2.2204460e-16 NaN 不定型(Not a Number) 内部变量的运算规则 inf 参与的运算 NaN 参与的运算,15,点运算(向量运算), a = 1 2 3 4; b = 0 1 2 3; a + b ans = 1 3 5 7 b./a ans = 0 0.5000 0.6667 0.7500 b.*a ans = 0 2 6 12 a.b ans = 1 2 9 64,16,矩阵操作入门,矩阵的抽取 A = 8 1 6; 3 5 7; 4 9 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,real,imag) 取整取余(floor,ceil,round,mod) 全部都具有向量功能: 输入向量, 返回对应函数值组成的向量,19,初等函数(exp), x = 0:0.2:1 x = 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.7621 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 = tril(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); end function y = fun1(x) function y = fun2(x),str2num(2*pi) ans = 6.2832 str2num(cos(2*pi) ans = 1 eval(cos(2*pi) ans = 1,26,逻辑判断,真值 逻辑真(true,1)或假(false,0) 逻辑运算 , =, =, =, = 例如 又如, pi = sqrt(10) ans = 0, 1:5 = pi ans = 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。 仅当输入向量 v 每个分量为假(零), any(v)返回0；否则返回1。 any(0 1 2 3) ans = 1 all(0 1 2 3) ans = 0,29,运算的级别,函数、括号 算术运算 逻辑运算 关系运算 每一个运算级别内仍有优先级的高低,30,matlab帮助,help help 命令名/工具箱名 help help doc doc eig what 列出当前matlab文件 who/whos 列出当前变量 which 寻找matlab文件或命令所在目录 lookfor 比较 lookfor curve / help curve,31,结构化编程,分支结构(if) for循环 while循环 嵌套结构 程序规范,32,分支结构(if),if 判断条件1 执行语句1 elseif 判断条件2 执行语句2 elseif 判断条件3 执行语句3 else 执行语句n+1 end,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 = magic(3), k+1 end,向量功能,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; end end,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 = mod(nn,10); rn = rn * 10 + d; nn = floor(nn/10); end disp(rn) if rn=n, flag = 1; else n = rn + n; end end,输入n 设置迭代指标为0 当迭代不成功时(=0), 保留当前n 计算n的逆序为rn 显示rn 重新设置迭代指标 或 n = rn + n end %迭代不成功,伪代码,源代码,对一个整数n, 把它反过 来读的数称为逆序数。 如132的逆序数为231。 写一个算法，输入n, 若 它与自己的逆序数相等, 则停止；否则把它与自 己的逆序数相加，重新 判断。,37,程序规范,程序按照结构缩进(matlab会自行处理) 按照逻辑关系分段 添加适当的注释(5%-30%) 规范所有的标识符(文件名、变量名) 适当增加空行、空格以保证可读性 办法: 参考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); end end,判别闰年不用 这么烦吧？,n = input(years: ); for k = n, if mod(k,400)=0 | . (mod(k,4)=0 end end,好像简单 多了.,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);,保存为fun.m,40,函数的调用, z,p = fun(15) z = 24 p = 4 x = fun(6) x = 12 fun(28) ans = 56, p,z = fun Input 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, fun(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 = 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,算法的复杂度,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*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 = linspace(-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=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 1550 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 -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) * b x = 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 b x = 3.0000 -24.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:4 A = 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 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.1789 D = 1.0000 0 0 0 2.0000 0 0 0 3.0000,68,最优化工具箱(linprog),什么是线性规划？, 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,),最优解 向量,求解 情况,工具箱 函数,目标函数 字符串,初值,算法参数,目标函数 梯度 字符串,目标函数 其他参数,help foptions,70,最优化工具箱(fminu),function v=f(x) 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,),零点向量,函数值,工具箱 函数,目标函数 字符串,初值,算法参数,目标函数 其他参数,help foptions,72,最优化工具箱(fzero), f = inline(x+exp(x)-3); x,fv=fzero(f,2) x = 0.79205996843068 fv = 0 只能求单变量函数的根！,73,最优化工具箱(fmincon),x = fmincon(fun,x0,. 1 1 1,8,. 1 2 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.),求解函数,工具箱 函数,算法参数,help odeset,求解区间,初值,解的向量 可用plot(x,y)画出解,目标函数 其他参数,75,微分方程工具箱(ode45),function dydx =f(x,y) dydx = y x + 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.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,method),nearest 最邻近插值 linear 双线性插值 cubic 双三次插值 缺省时, 双线性插值,80,统计工具箱,参数估计 betafit, binofit, mle, 概率密度函数 betapdf, binopdf, hygepdf, 各种分布的均值、方差 betastat, binostat, chi2stat, 线性、非线性的各种回归模型 anova1, anova2, polyf