matlab 科学工程计算

matlab科学工程计算 同济大学数学系陈雄达 chenxiongda@gmail.com http://math.tongji.edu.cn/model/docs/matlabc.pps

teapotdemo

目录 • matlab入门 • matlab基础编程 • matlab作图 • matlab工具箱 • matlab上机操作 • matlab综合训练 • matlab命令速查

matlab入门 • matlab是什么 • matlab不是什么 • matlab集成环境 • matlab ABC

matlab入门 没有matlab就没有乐趣。 M.N. Nachtigal S.C. Reddy L.N. Trefethen 关于迭代法的Copper Mountain 论文集 1990

matlab是什么 • 一个可视化的计算程序, 广泛使用于从个人计算机到超级计算机范围内的各种计算机上 • 包括命令控制、可编程，上百个预先定义命令和函数 • 有许多强有力的命令, 能完成大量的高级矩阵处理 • 强有力的二维、三维图形工具 • 能与其他程序一起使用 • 25个(不断增加中)不同的工具箱应用于特殊的应用领域 • 工业研究与开发的有力工具 • 数学教学, 尤其线代, 数值分析, 科学计算方面的教研工具 • 电子学, 控制理论, 物理学等工程科学方面的教研工具 • 经济学, 化学和生物学等有计算问题的所有领域中的教学与研究 • 名字取自矩阵实验室(matrix laboratory)

matlab不是什么 • 不是万能的解决工具 • 不是最高性能的编程语言 • 受计算条件限制, 不能解决超大型实际问题 • 不能解决工具箱之外的问题种类 • -- 需要编写接口、算法甚至工具箱

集成环境 • 窗口系统 • View > Desktop Layout > Five Panel • 历史命令 / 变量和文件管理 / 命令窗口 • 菜单系统 • File / Edit / View / Web / Window / Help • 按钮

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];

数据结构 冒号(:) 矩阵操作入门标识符初等函数(elfun) 初等函数(exp) 函数fix,round, ceil,floor 逻辑判断关系运算运算的级别 matlab帮助结构化编程脚本文件函数文件 matlab基础编程

数据结构 • 最基本的数据结构-- 矩阵 • 数和向量看成为特殊的矩阵 • 矩阵以[ ]为定界符，与多维数组等同 • 字符串看成为每个元素都是单个字符的向量，也可以有字符矩阵 • 高维数组 • 细胞 (cell) • 结构体 (struct)

冒号(:) • 冒号 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 [ ] (空矩阵)

标识符 • 文件名、变量名、函数名 • 标识符的规定 • 最长不超过19个字符 • 为配合matlab的风格，采用见名知意的小写短名称 • 以字母开头，包含数字、大小写字母、下划线 • 变量命名规则 • 常用从简，专用从繁 • 矩阵大写 A, B, C • 向量小写 u, v, w, x, y, z • 函数小写 f, g, h, fun, f1, f2 • 常量小写 alpha, beta, a, b, c

内部变量 • pi 圆周率=3.1415926 • i,j 虚根 • inf 无限大 • eps 2.2204460e-16 • NaN 不定型(Not a Number) • 内部变量的运算规则 • inf 参与的运算 • NaN 参与的运算

数据及其显示格式 • format (short) 短格式，4位小数 • format long (e/g) 长格式，15位小数 • format + 只显示符号 • format bank 银行格式 • format rat 分数格式 • formatloose 大间距 • format compact 小间距

输出 >> A=magic(3); >> disp(A); >> s = ’hello world’; >> disp(s); >> v = 1:100; >> disp(v); >> fprintf(’%2d%2d%2d\n’,A); >> fprintf(’pi^2 = %12.6f\n’,pi^2); >> fid = fopen(’h.txt’,’w’); >> fprintf(fid,’s = %s\n’,’hello world’); >> fclose(fid);

点运算(向量运算) >> 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

矩阵操作入门 • 矩阵的抽取 • 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]

A(2)= A(n-1,:)= A(2+3, 2*3)= x= A(sqrt(4),x-3)= A(:)= v(end)= A(2) + 2 = x^2 + x = A(sqrt(4),x)-3 = A’ = inv(A) = (x>=0) = 赋值语句左值表达式: 指向某一存储空间的表达式，放在赋值语句左边。

初等函数(elfun) • 三角函数(sin,sinh,asin,asinh) • 指数对数(exp,log,log10,sqrt) • 复变函数(abs,conj,real,imag) • 取整取余(floor,ceil,round,mod) • 全部都具有向量功能: 输入向量, 返回对应函数值组成的向量

初等函数(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 向量功能

函数fix,round,ceil,floor

几个有用的函数 >> [a,ind] = sort([41 33 16 10 13 54] a = 10 13 16 33 41 54 ind = 4 5 3 2 1 6 >> a = find([41 33 16 10 13 54]>=35) a = 1 6 >> [a,ind] = max([sqrt(10) 22/7 pi]) a = 3.1623 ind = 1

简单矩阵操作 rand(n), eye(n), zeros(n) 生成 n 阶随机、单位、零矩阵 >> eye(2,3) ans = 1 0 0 0 1 0 >> zeros(3) ans = 0 0 0 0 0 0 0 0 0 >> 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

>> 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 diag

>> 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 triu/tril

>> 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 字符运算字符串

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 字符运算

逻辑判断 • 真值逻辑真(true,1)或假(false,0) • 逻辑运算 >, <, >=, <=, ==, ~= • 例如 • 又如 >> pi >= sqrt(10) ans = 0 >> [1:5] >= pi ans = 0 0 0 1 1 向量功能

关系运算 • 关系运算有： • 与(and) • and(a,b) • a & b • 或(or) • or(a,b) • a | b • 非(not) • not(a) • ~a • 异或(xor) • xor(a,b)

两个特殊关系运算(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

运算的级别 • 函数、括号 • 算术运算 • 逻辑运算 • 关系运算 • 每一个运算级别内仍有优先级的高低

matlab帮助 • help • help 命令名/工具箱名 • help help • doc • doc eig • what • 列出当前matlab文件 • who/whos • 列出当前变量 • which • 寻找matlab文件或命令所在目录 • lookfor • 比较 lookfor curve / help curve

结构化编程 • 分支结构(if,switch) • for循环 • while循环 • 嵌套结构 • 程序规范

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 分支结构(if)

分支结构(switch) switch month case {1,3,5,7,8,10,12} day = 31; case {4,6,9,11} day = 30; otherwise day = 28; % 假设不是闰年 end

for 循环变量 = 循环值 循环操作 end s = 0; for k = 1:100, s = s + 1/k^2; end for循环 for k = magic(3), k+1 end for k = ’hello world’, fprintf(’%c’,k); pause(0.1); end 向量功能

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 while循环 Collatz猜想

嵌套结构 for k = 1:9, for j = 1:9, fprintf(’%d*%d=%3d ’, k,j,k*j); end fprintf(’\n’); end 九九乘法表

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, 把它反过来读的数称为逆序数。如132的逆序数为231。写一个算法，输入n, 若它与自己的逆序数相等, 则停止；否则把它与自己的逆序数相加，重新判断。输入n 设置迭代指标为0 当迭代不成功时(=0), 保留当前n 计算n的逆序为rn 显示rn 重新设置迭代指标或 n <= rn + n end %迭代不成功源代码伪代码

程序规范 • 程序按照结构缩进(matlab会自行处理) • 按照逻辑关系分段 • 添加适当的注释(5%--30%) • 规范所有的标识符(文件名、变量名) • 适当增加空行、空格以保证可读性 • 办法: 参考Matlab内部程序风格

脚本文件 n = input(’years: ’); for k = n, if mod(k,400)==0 | ... (mod(k,4)==0&mod(k,100)~=0), disp([num2str(k) ’is a leap year’]); else disp([num2str(k) ’isn’’t a leap year’]); end end 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) ’isn’’t a leap year’]); elseif mod(k,4)==0, disp([num2str(k) ’is a leap year’]); else disp([num2str(k) ’isn’’t a leap year’]); end end 好像简单多了... 判别闰年不用这么烦吧？

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

>> [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)

函数的相互调用 寻找亲和数对, 即这样两个数, 它们的真因子之和等于对方。 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)); 子函数

递归调用 函数调用自己称为递归—因此必须有停止自调用的时刻。例：Frayer级数, 即[0,1]上分母不超过n的既约真分数从小到大排列。每一次把新出现的分数, 放在两个已有分数中, 这两个分数的分子和分母各自的和就是新分数的分子和分母。

递归调用 function [p,q] = frayer(n) if n>1, [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

递归调用 function v = fact(n) if n==0|n==1, v = 1; else v = n * fact(n-1); end 阶乘 >> v = prod(1:n);

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

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

matlab 科学工程计算