matlabmatlab科学工程计算科学工程计算科学工程计算科学工程计算同济大学数学系同济大学数学系同济大学数学系同济大学数学系 陈雄达陈雄达陈雄达陈雄达chenxiongdagmail.comchenxiongdagmail.comteapotdemo2目录目录matlab入门matlab基础编程matlab作图matlab工具箱matlab上机操作matlab综合训练matlab命令速查3matlab入门matlab是什么matlab不是什么matlab集成环境matlab ABC4matlab入门 没有matlab就没有乐趣。 M.N. Nachtigal S.C. Reddy L.N. Trefethen 关于迭代法的Copper Mountain 论文集 19905matlab是什么一个可视化的计算程序, 广泛使用于从个人计算机到超级计算机范围内的各种计算机上包括命令控制、可编程，上百个预先定义命令和函数有许多强有力的命令, 能完成大量的高级矩阵处理强有力的二维、三维图形工具能与其他程序一起使用25个(不断增加中)不同的工具箱应用于特殊的应用领域工业研究与开发的有力工具数学教学, 尤其线代, 数值分析, 科学计算方面的教研工具电子学, 控制理论, 物理学等工程科学方面的教研工具经济学, 化学和生物学等有计算问题的所有领域中的教学与研究名字取自矩阵实验室(matrix laboratory)6matlab不是什么不是万能的解决工具不是最高性能的编程语言受计算条件限制, 不能解决超大型实际问题不能解决工具箱之外的问题种类- 需要编写接口、算法甚至工具箱7集成环境窗口系统View Desktop Layout Five Panel历史命令 / 变量和文件管理 / 命令窗口菜单系统File / Edit / View / Web / Window / Help按钮8matlab 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;9matlab基础编程数据结构冒号(:)矩阵操作入门标识符初等函数(elfun)初等函数(exp)函数fix,round,ceil,floor逻辑判断关系运算运算的级别matlab帮助结构化编程脚本文件函数文件10数据结构最基本的数据结构- 矩阵数和向量看成为特殊的矩阵矩阵以 为定界符，与多维数组等同字符串看成为每个元素都是单个字符的向量，也可以有字符矩阵高维数组细胞 (cell)结构体 (struct)11冒号(:)冒冒号号 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 (空矩空矩阵阵)12标识符文件名、变量名、函数名文件名、变量名、函数名标识符的规定标识符的规定最长不超过最长不超过1919个字符个字符为配合为配合matlabmatlab的风格，采用见名知意的小写短名称的风格，采用见名知意的小写短名称以字母开头，包含数字、大小写字母、下划线以字母开头，包含数字、大小写字母、下划线变量命名规则变量命名规则常用从简，专用从繁常用从简，专用从繁矩阵大写矩阵大写 A, B, CA, B, C向量小写向量小写 u, v, w, x, y, zu, v, w, x, y, z函数小写函数小写 f, g, h, fun, f1, f2f, g, h, fun, f1, f2常量小写常量小写 alpha, beta, a, b, calpha, beta, a, b, c13内部变量pi 圆周率圆周率=3.1415926i,j 虚根虚根inf 无限大无限大eps 2.2204460e-16NaN 不定型不定型(Not a Number)内部变量的运算规则内部变量的运算规则inf 参与的运算参与的运算NaN 参与的运算参与的运算14点运算(向量运算) 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 6415矩阵操作入门矩阵的抽取矩阵的抽取 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; 28 81 16 63 35 57 74 49 92 216赋值语句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) = 左值表达式左值表达式: 指向某一存储空间的表达式，放在赋值语句左边。指向某一存储空间的表达式，放在赋值语句左边。17初等函数(elfun)三角函数三角函数(sin,sinh,asin,asinh)指数对数指数对数(exp,log,log10,sqrt)复变函数复变函数(abs,conj,real,imag)取整取余取整取余(floor,ceil,round,mod)全部都具有全部都具有向量功能向量功能: : 输入向量输入向量, , 返回对应函数值组成的向量返回对应函数值组成的向量18初等函数(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向量功向量功能能19函数fix,round,ceil,floorfix,round,ceil,floorfixroundceilfloor功能功能朝零取整朝零取整四舍五入四舍五入向上取整向上取整向下取整向下取整3.141533430.61800110-2.7181-2-3-2-3-1.4142-1-1-1-220简单矩阵操作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 021diag 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 022triu/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 223字符运算 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字符串24字符运算for k = 1:9, st = fun . num2str(k); x(k) = . feval(st, x(k);endfunction y = fun1(x)function y = fun2(x)str2num(2*pi)ans = 6.2832str2num(cos(2*pi)ans = 1eval(cos(2*pi)ans = 125逻辑判断真真值值 逻辑逻辑真真(true,1)(true,1)或假或假(false,0)(false,0)逻辑逻辑运运算算 , =, pi = sqrt(10) ans = 0 1:5 = pians = 0 0 0 1 1向量功向量功能能26关系运算关关系系运运算有：算有：与与(and)(and)and(a,b)a & b或或(or)(or)or(a,b)a | b非非(not)(not)not(a)a异异或或( (xorxor) )xor(a,bxor(a,b) )a0011b0101a&b0001a|b0111a1100xor(a,b)011027两个特殊关系运算(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 = 028运算的级别函数、括号函数、括号算术运算算术运算逻辑运算逻辑运算关系运算关系运算每一个运算级别内仍有优先级的高低每一个运算级别内仍有优先级的高低29matlab帮助helphelp 命令名命令名/工具箱名工具箱名help helpdocdoc eigwhat列出列出当当前前matlab文件文件who/whos列出列出当当前前变变量量which寻寻找找matlab文件或命令所在目文件或命令所在目录录lookfor比比较较 lookfor curve / help curve30结构化编程分支结构(if)for循环while循环嵌套结构程序规范31分支结构(if)if 判判断条断条件件1 执执行行语语句句1elseif 判判断条断条件件2 执执行行语语句句2elseif 判判断条断条件件3 执执行行语语句句3else 执执行行语语句句n+1endif temp 10, disp(cold!);elseif temp 20, disp(cool);elseif temp 27, disp(warm);elseif temp =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; endend 输入输入n 设置迭代指标为设置迭代指标为0 当迭代不成功时当迭代不成功时(=0), 保留当前保留当前n 计算计算n的逆序为的逆序为rn 显示显示rn 重新设置迭代指标重新设置迭代指标 或或 n 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)40函数的相互调用寻找亲和数对寻找亲和数对, , 即这样两个数即这样两个数, , 它们的真因子之和等于对方。它们的真因子之和等于对方。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 endfunction p = fun(n) k = 1:n-1; z = mod(n,k); p = sum(k(z=0);调用调用子函数子函数41递归调用函数调用自己称为递归函数调用自己称为递归因此必须有停止自调用的时刻。因此必须有停止自调用的时刻。0 01 1例：例：Frayer级数级数, 即即0,1上分母不超过上分母不超过n的既约真分数的既约真分数从小到大排列。从小到大排列。1 14 41 13 31 12 22 23 34 45 53 34 41 11 11 15 52 25 53 35 5每一次把新出现的分数每一次把新出现的分数, 放在两个已有分数中放在两个已有分数中, 这两个这两个分数的分子和分母各自的和就是新分数的分子和分母。分数的分子和分母各自的和就是新分数的分子和分母。42递归调用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; end43递归调用function v = fact(n) if n=0|n=1, v = 1; else v = n * fact(n-1); end v = prod(1:n);阶乘阶乘44算法的复杂度s = 0;for k = 1:n, p = 1; for j = 1:k, p = p * j; end s = s + 1/(p+1);endop = 0;for k = 1:n, op = op + 1; for j = 1:k, op = op + 1; end op = op + 3;endop=4n+n(n+1)/245算法的复杂度s = 0; p = 1;for k = 1:n, p = p*k; s = s + 1/(p+1);endop = 2;for k = 1:n, op = op + 1; op = op + 3;endop=4n+246matlab作图二维作图三维作图其它作图47二维作图(plot)plot(x,y,s)例如例如 x = linspace(0,2*pi); y = sin(x); plot(x,y,r-);48二维作图(plot)plot(x,y,s) - s的规定的规定colorpoint styleline styleRed. point- solidGreenx x-mark: dottedBlueo circle-. dashdotCyan+ plus- dashedMagenta* starYellows,p,h 正正456边边形形blacKd diamondv, 三角形三角形49二维作图(plot)xlabel(x);ylabel(y);title(y=sin(x);50二维作图(polar) theta = linspace(0,6*pi); rho = theta; polar(theta,rho,r-);51三维作图(plot3) plot3(x,y,z,s) t = 0:pi/50:10*pi; plot3(sin(t),cos(t),t);52三维作图(comet3) t = 0:pi/500:20*pi; comet3(sin(t),cos(t),t);53三维作图(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 754三维作图(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有同样的格式有同样的格式向量功向量功能能55其它作图(bar3) A = 82 82 63 81 82 81 84 61 84 85 80 79 65 84 86; bar3(A,grouped);56其它作图(pie3)pie3(2 4 3 5,0 1 1 0,. North,South,East,West)57surf, contour例例: 山区地貌山区地貌- 在某山区测得一些地点的高程如下表。平在某山区测得一些地点的高程如下表。平面区域为面区域为 1200=x=4000,1200=y 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);59matlab工具箱(toolbox)矩阵计算矩阵计算(help matfun)最优化工具箱最优化工具箱(help optim)微分方程工具箱微分方程工具箱(help ode45,help pde)多项式计算多项式计算(help polyfun)统计工具箱统计工具箱(help stats)60矩阵计算(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.102861矩阵计算() 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.0000 30.000062矩阵的秩(rank) A = magic(4)A = 16 2 3 13 5 11 10 8 9 7 6 12 4 14 15 1 rank(A)ans = 363矩阵的行列式(det) A = pascal(4)A = 1 1 1 1 1 2 3 4 1 3 6 10 1 4 10 20 det(A)ans = 164矩阵的迹(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 = 3065矩阵的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 166矩阵的特征值(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.000067最优化工具箱(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.000068最优化工具箱(fminu)x,opts=fminu(f,x0,options,g,p1,p2,)最优解最优解向量向量求解求解情况情况工具箱工具箱函数函数目标函数目标函数字符串字符串 初值初值算法算法参数参数目标函数目标函数梯度梯度字符串字符串 目标函数目标函数其他参数其他参数 helpfoptions69最优化工具箱(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.000070最优化工具箱(fzero)x,fv=fzero(f,x0,options,p1,p2,)零点向量零点向量函数值函数值工具箱工具箱函数函数目标函数目标函数字符串字符串 初值初值算法算法参数参数目标函数目标函数其他参数其他参数 helpfoptions71最优化工具箱(fzero) f = inline(x+exp(x)-3); x,fv=fzero(f,2)x = 0.79205996843068fv = 0只能求单变量函数的根！只能求单变量函数的根！72最优化工具箱(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 073微分方程工具箱(ode45)t,y = ode45(odefun,tspan,y0,options,p1,p2.)求解函数求解函数工具箱工具箱函数函数算法算法参数参数helpodeset求解区间求解区间初值初值解的向量解的向量可用可用plot(x,y)画出解画出解目标函数目标函数其他参数其他参数 74微分方程工具箱(ode45)function dydx =f(x,y) dydx = y x + 1; x,y = . ode45(f,0,2,1); plot(x,y,r-);75多项式计算多项式表示：向量多项式表示：向量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)76多项式计算多项式拟合多项式拟合(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)77插值(interp1)yi=interp1(x，y，xi，method)插值方法插值方法被插值点被插值点插值节点插值节点插值结果插值结果nearest: 最邻近插值最邻近插值linear: 线性插值线性插值spline: 三次样条插值三次样条插值 cubic: 立方插值立方插值缺省时缺省时: 分段线性插值分段线性插值78插值(interp2)z=interp2(x0,y0,z0,x,y,method)被插值点插值方法插值节点nearestnearest 最邻近插值最邻近插值linearlinear 双线性插值双线性插值cubiccubic 双三次插值双三次插值缺省时缺省时, , 双线性插值双线性插值被插值点的函数值79统计工具箱参数估计参数估计betafitbetafit, , binofitbinofit, , mlemle, , 概率密度函数概率密度函数betapdfbetapdf, , binopdfbinopdf, , hygepdfhygepdf, , 各种分布的均值、方差各种分布的均值、方差betastatbetastat, , binostatbinostat, chi2stat, , chi2stat, 线性、非线性的各种回归模型线性、非线性的各种回归模型anova1, anova2, anova1, anova2, polyfitpolyfit, regress, , regress, lsqnonneglsqnonneg, , nlinfitnlinfit, , nlintoolnlintool, , 80matlab上机操作遵守机房各项规定遵守机房各项规定保持机房卫生整洁保持机房卫生整洁matlabmatlab系统若有问题，请重新开机，在选择系统若有问题，请重新开机，在选择进入系统进入系统前前按按Ctrl-RCtrl-R81matlab练习题求下面表达式的值求下面表达式的值 其中其中82matlab练习题已知矩阵已知矩阵A,BA,B如下如下, , 计算并比较计算并比较 A*B, B*A, A.*B, B.*A, A/B, AB, A./B, A.B83matlab练习题已知矩阵已知矩阵A,BA,B如下如下, , 给出这两个矩阵的给出这两个矩阵的 (1)(1)逆逆 (2)(2)行列式行列式 (3)(3)迹迹 (4)(4)特征值和特征向量特征值和特征向量 (5) (5)最简阶梯形矩阵最简阶梯形矩阵84matlab练习题已知矩阵已知矩阵A A及多项式如下及多项式如下, , 给出多项式在矩阵上的值；给出多项式在矩阵上的值； 比较结果矩阵和原矩阵的特征值和特征向量。比较结果矩阵和原矩阵的特征值和特征向量。85matlab综合训练矩阵矩阵: : 写一个程序，对某一输入的对称矩阵写一个程序，对某一输入的对称矩阵A A进行进行分解分解, A=LDL, , A=LDL, 其中其中L L为单位下三角矩阵，为单位下三角矩阵，D D为对角为对角矩阵。矩阵。 参考参考数值计算数值计算中的中的CholeskyCholesky分解。分解。画图画图: : 给出某一曲面方程及一垂直于给出某一曲面方程及一垂直于XOYXOY平面的平平面的平面方程，利用平面把曲面割成两片，画出这两片，面方程，利用平面把曲面割成两片，画出这两片，中间的缝隙要逐渐变大。中间的缝隙要逐渐变大。算法算法: : 输入平面上的至少四个点，挑选当中的部分输入平面上的至少四个点，挑选当中的部分点，并把它们连成一个多边形，要求面积尽量大，点，并把它们连成一个多边形，要求面积尽量大，给出其面积之值。给出其面积之值。86matlab命令速查 : .* ./ . 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 function 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 87作业: 三维画图 在某海域测得一些点在某海域测得一些点(x,y)(x,y)处的水深处的水深z z由下由下表给出，船的吃水深度为表给出，船的吃水深度为5 5英尺，在矩形区域英尺，在矩形区域（7575，200200）* *（-50-50，150150）里的哪些地方船）里的哪些地方船要避免进入。要避免进入。x129.0 140.0 103.588.0185.5 195.0 105.5y7.5141.523.0147.022.5137.585.5z4868688x157.5 107.577.081.0162.0 162.0 117.5y-6.5-81.03.056.5-66.584.0-33.5z998894988作业: 三维画图89THANK YOU90