Transcript matlab科学工程计算 同济大学数学系 陈雄达 [email protected] 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 论文集5 matlab是什么  一个可视化的计算程序, 广泛使用于从个人计算机到超级计算机 范围内的各种计算机上  包括命令控制、可编程，上百个预先定义命令和函数  有许多强有力的命令, 能完成大量的高级矩阵处理  强有力的二维、三维图形工具  能与其他程序一起使用  25个(不断增加中)不同的工具箱应用于特殊的应用领域 

matlab科学工程计算
同济大学数学系 陈雄达
[email protected]
http://math.tongji.edu.cn/model/docs/matlabc.pps
teapotdemo
2
目录
matlab入门
matlab基础编程
matlab作图
matlab工具箱
matlab上机操作
matlab综合训练
matlab命令速查
3
matlab入门
 matlab是什么
 matlab不是什么
 matlab集成环境
 matlab ABC
4
matlab入门
没有matlab就没有乐趣。
M.N. Nachtigal
S.C. Reddy
L.N. Trefethen
关于迭代法的Copper Mountain 论文集
1990
5
matlab是什么
 一个可视化的计算程序, 广泛使用于从个人计算机到超级计算机
范围内的各种计算机上
 包括命令控制、可编程，上百个预先定义命令和函数
 有许多强有力的命令, 能完成大量的高级矩阵处理
 强有力的二维、三维图形工具
 能与其他程序一起使用
 25个(不断增加中)不同的工具箱应用于特殊的应用领域
 工业研究与开发的有力工具
 数学教学, 尤其线代, 数值分析, 科学计算方面的教研工具
 电子学, 控制理论, 物理学等工程科学方面的教研工具
 经济学, 化学和生物学等有计算问题的所有领域中的教学与研究
 名字取自矩阵实验室(matrix laboratory)
6
matlab不是什么
不是万能的解决工具
不是最高性能的编程语言
受计算条件限制, 不能解决超大型实际问题
不能解决工具箱之外的问题种类
-- 需要编写接口、算法甚至工具箱
7
集成环境
 窗口系统
– View > Desktop Layout > Five Panel
– 历史命令 / 变量和文件管理 / 命令窗口
 菜单系统
 File / Edit / View / Web / Window / Help
 按钮
8
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];
9
matlab基础编程
 数据结构
 冒号(:)
 矩阵操作入门
 标识符
 初等函数(elfun)
 初等函数(exp)
 函数fix,round,
ceil,floor
 逻辑判断
 关系运算
 运算的级别
 matlab帮助
 结构化编程
 脚本文件
 函数文件
10
数据结构
 最基本的数据结构-- 矩阵
 数和向量看成为特殊的矩阵
 矩阵以[ ]为定界符，与多维数组等同
 字符串看成为每个元素都是单个字符的向量，
也可以有字符矩阵
 高维数组
 细胞 (cell)
 结构体 (struct)
11
冒号(:)
 冒号
a:s:b
 例如
1:2:10
7:-2:0
1:6
6:2
从a开始, 步长为s, 界为b
[1 3 5 7 9]
[7 5 3 1]
[1 2 3 4 5 6]
[ ]
(空矩阵)
12
标识符
 文件名、变量名、函数名
 标识符的规定
 最长不超过19个字符
 为配合matlab的风格，采用见名知意的小写短名称
 以字母开头，包含数字、大小写字母、下划线
 变量命名规则
 常用从简，专用从繁
 矩阵大写 A, B, C
 向量小写 u, v, w, x, y, z
 函数小写 f, g, h, fun, f1, f2
 常量小写 alpha, beta, a, b, c
13
内部变量
 pi
 i,j
 inf
 eps
 NaN
圆周率=3.1415926
虚根
无限大
2.2204460e-16
不定型(Not a Number)
 内部变量的运算规则
inf 参与的运算
NaN 参与的运算
14
数据及其显示格式
 format
 format
 format
 format
 format
(short)
long (e/g)
+
bank
rat
 format loose
 format compact
短格式，4位小数
长格式，15位小数
只显示符号
银行格式
分数格式
大间距
小间距
15
输出
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
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);
16
点运算(向量运算)
>> 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
>> b.*a
ans =
0
2
6
12
>> a.^b
ans =
1
2
9
64
0.7500
17
矩阵操作入门
 矩阵的抽取
 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]
8
1
6
3
5
7
4
9
2
18
赋值语句
左值表达式: 指向某一存储空间的表达式，放在赋值语句左边。
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) =
19
初等函数(elfun)
 三角函数(sin,sinh,asin,asinh)
 指数对数(exp,log,log10,sqrt)
 复变函数(abs,conj,real,imag)
 取整取余(floor,ceil,round,mod)
 全部都具有向量功能:
输入向量, 返回对应函数值组成的向量
20
初等函数(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
向量功
能
21
函数fix,round,ceil,floor
fix
round
ceil
floor
功能
朝零取整
四舍五入
向上取整
向下取整
3.1415
3
3
4
3
0.6180
0
1
1
0
-2.7181
-2
-3
-2
-3
-1.4142
-1
-1
-1
-2
22
几个有用的函数
>> [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
23
简单矩阵操作
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
24
diag
>> diag([1
7
ans = 1
5
9
>> diag([1
7
ans = 2
6
2 3; 4 5 6
8 9])
2 3; 4 5 6
8 9],1)
>> 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
25
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
26
字符运算
>> 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
字符串
27
字符运算
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
28
逻辑判断
 真值 逻辑真(true,1)或假(false,0)
 逻辑运算 >, <, >=, <=, ==, ~=
 例如
>> pi >= sqrt(10)
ans = 0
 又如
>> [1:5] >= pi
ans = 0 0 0 1 1
向量功能
29
关系运算
 关系运算有：
 与(and)
 and(a,b)
a & b
 或(or)
 or(a,b)
a | b
 非(not)
 not(a)
 ~a
 异或(xor)
 xor(a,b)
a
0
0
1
1
b
0
1
0
1
a&b
0
0
0
1
a|b
0
1
1
1
~a
1
1
0
0
xor(a,b)
0
1
1
0
30
两个特殊关系运算(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
31
运算的级别
 函数、括号
 算术运算
 逻辑运算
 关系运算
 每一个运算级别内仍有优先级的高低
32
matlab帮助
 help
 help 命令名/工具箱名
 help help
 doc
 doc eig
 what
 列出当前matlab文件
 who/whos
 列出当前变量
 which
 寻找matlab文件或命令所在目录
 lookfor
 比较 lookfor curve
/
help curve
33
结构化编程
 分支结构(if,switch)
 for循环
 while循环
 嵌套结构
 程序规范
34
分支结构(if)
判断条件1
执行语句1
elseif 判断条件2
执行语句2
elseif 判断条件3
执行语句3
else
执行语句n+1
end
if
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
35
分支结构(switch)
switch month
case {1,3,5,7,8,10,12}
day = 31;
case {4,6,9,11}
day = 30;
otherwise
day = 28; % 假设不是闰年
end
36
for循环
for 循环变量 = 循环值
循环操作
end
for k = magic(3),
k+1
end
向量功
能
s = 0;
for k = 1:100,
s = s + 1/k^2;
end
for k = ’hello world’,
fprintf(’%c’,k);
pause(0.1);
end
37
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
猜想
38
嵌套结构
for k = 1:9,
for j = 1:9,
fprintf(’%d*%d=%3d ’, k,j,k*j);
end
fprintf(’\n’);
九九乘法表
end
39
嵌套结构
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 %迭代不成功
源代码
伪代码
40
程序规范
 程序按照结构缩进(matlab会自行处理)
 按照逻辑关系分段
 添加适当的注释(5%--30%)
 规范所有的标识符(文件名、变量名)
 适当增加空行、空格以保证可读性
 办法: 参考Matlab内部程序风格
41
脚本文件
n = input(’years: ’);
for k = n,
mod(k,400)==0,| ...
if mod(k,400)==0
disp([num2str(k)
’is a leap year’]);
(mod(k,4)==0&mod(k,100)~=0),
elseif
mod(k,100)==0,
disp([num2str(k)
’is a leap year’]);
disp([num2str(k) ’isn’’t a leap year’]);
else
elseif
mod(k,4)==0,’isn’’t a leap year’]);
disp([num2str(k)
disp([num2str(k) ’is a leap year’]);
end
else
end
disp([num2str(k) ’isn’’t a leap year’]);
end
end
好像简单
判别闰年不用
多了...
这么烦吧？
42
函数文件
function [p,t] = fun(n)
% n: input integer
% p: sum of factors
% t: number of factors
p = 0;
保存为fun.m
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));
43
函数的调用
>> [z,p] = fun(15)
z = 24
p = 4
>> x = fun(6)
x = 12
>> fun(28)
ans = 56
>> [p,z] = fun
Input argument ‘n’ is undefin
ed.
>> [p,z] = fun(20,18)
Too many input arguments.
>> [p,t,b] = fun(15)
Too many output arguments.
function [p,t] = fun(n)
44
函数的相互调用
寻找亲和数对, 即这样两个数, 它们的真因子之和等于对方。
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));
45
递归调用
函数调用自己称为递归—因此必须有停止自调用的时刻。
例：Frayer级数, 即[0,1]上分母不超过n的既约真分数
从小到大排列。
0
1
1
1
2
1
3
2
3
4
1
1
5
4
3
5
2
5
3
4
5
1
每一次把新出现的分数, 放在两个已有分数中, 这两个
分数的分子和分母各自的和就是新分数的分子和分母。
46
递归调用
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
47
递归调用
function v = fact(n)
if n==0|n==1,
v = 1;
else
v = n * fact(n-1);
end
>> v = prod(1:n);
阶乘
48
算法的复杂度
s = 0;
for k = 1:n,
p = 1;
n
for j = 1:k,
p = p * j;
k 1
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
1
求
k!1
op=4n+n(n+1)/2
49
算法的复杂度
s = 0; p = 1;
for k = 1:n,
p = p*k;
s = s + 1/(p+1);
end
n
1
求
k 1 k!1
op = 2;
for k = 1:n,
op = op + 1;
op = op + 3;
end
op=4n+2
50
matlab作图
 二维作图
 三维作图
 其它作图
51
二维作图(plot)
 plot(x,y,s)
 例如
x = linspace(0,2*pi);
y = sin(x);
plot(x,y,’r-’);
1
0.8
0.6
0.4
0.2
0
-0.2
-0.4
-0.6
-0.8
-1
0
1
2
3
4
5
6
7
52
二维作图(plot)
 plot(x,y,s) -- s的规定
color
point style
line style
Red
.
point
-
solid
Green
x
x-mark
:
dotted
Blue
o
circle
-.
dashdot
Cyan
+
plus
--
dashed
Magenta
*
star
Yellow
s,p,h
正456边形
d
diamond
blacK
v,^,<,>
三角形
53
二维作图(plot)
 xlabel(’x’);
 ylabel(’y’);
 title(’y=sin(x)’);
y=sin(x)
1
0.8
0.6
0.4
y
0.2
0
-0.2
-0.4
-0.6
-0.8
-1
0
1
2
3
4
5
6
7
x
54
二维作图(polar)
>> theta = linspace(0,6*pi);
>> rho = theta;
>> polar(theta,rho,‘r-’);
90
20
120
150
60
15
30
10
5
180
0
210
330
240
300
270
55
三维作图(plot3)
 plot3(x,y,z,s)
>> t = 0:pi/50:10*pi;
>> plot3(sin(t),cos(t),t);
35
30
25
20
15
10
5
0
1
1
0.5
0.5
0
0
-0.5
-0.5
-1
-1
56
三维作图(comet3)
>> t = 0:pi/500:20*pi;
>> comet3(sin(t),cos(t),t);
70
60
50
40
30
20
10
0
1
1
0.5
0.5
0
0
-0.5
-0.5
-1
-1
57
三维作图(meshgrid)
 [X,Y] = meshgrid(x,y)
 例如 [X，Y] = meshgrid(1:3, 4:7)
则 X =[1
1
1
1
2
2
2
2
且 Y = [4
5
6
7
3
3
3
3]
8
7.5
7
6.5
6
5.5
5
4
5
6
7
4
5
6
7]
4.5
4
3.5
3
0
0.5
1
1.5
2
2.5
3
3.5
4
58
三维作图(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有同样的格式
59
参数曲面
theta = linspace(0,2*pi);
psi
= linspace(-pi/2+0.5,pi/2-0.5);
[Theta,Psi] = meshgrid(theta,psi);
X = 4 * sec(Psi) .* cos(Theta);
Y = 4 * sec(Psi) .* sin(Theta);
Z = 3 * tan(Psi);
surf(X,Y,Z)
view(-41,22);
shading interp;
colormap spring
60
其它作图(bar3)
>> A = [82 82 63 81
81 84 61 84
80 79 65 84
>> bar3(A,'grouped');
82
85
86];
100
80
60
40
20
0
1
2
3
61
其它作图(pie3)
pie3([2 4 3 5],[0 1 1 0],...
{'North','South','East','West'})
West
North
East
South
62
surf, contour
例: 山区地貌-- 在某山区测得一些地点的高程如下表。平
面区域为 1200<=x<=4000,1200<=y<=3600)。试作出该山
区的地貌图和等高线图。
X 1200
1600
2000
2400
2800
3200
3600
4000
1200
1130
1250
1280
1230
1040
900
500
700
1600
1320
1450
1420
1400
1300
700
900
850
2000
1390
1500
1500
1400
900
1100
1060
950
2400
1500
1200
1100
1350
1450
1200
1150
1010
2800
1500
1200
1100
1550
1600
1550
1380
1070
3200
1500
1550
1600
1550
1600
1600
1600
1550
3600
1480
1500
1550
1510
1430
1300
1200
980
Y
63
surf, contour
>> A = [ 1130 1250 1280
1230
3500
1320 1450 2000
1420 1400
1390 1500 1500
1500 1400
3000
1500 1200 1000
1100 1350
500
1500 1200 1100
1550
4000
2500
1500 1550 1600
1550
3500
1480 1500 1550 1510
3000
>> x = 1200:400:4000;
2000
>> y = 1200:400:3600; 2500
>> surf(x,y,A);
2000
1500
>> contour(x,y,A,100);
1500
1040
1300
900
1450
1600
1600
1430
900
700
1100
1200
1550
1600
1300
500 700
900 850
1060 950
1150 1010
1380 1070
1600 1550
1200 980];
1500
1400
1300
1200
1100
1000
900
800
700
4000
3500
3000
2500
1500
1000 2000
1000
1500
2500
2000
3000
3500
4000
64
600
matlab工具箱(toolbox)
 矩阵计算(help matfun)
 最优化工具箱(help optim)
 微分方程工具箱(help ode45,help pde)
 多项式计算(help polyfun)
 统计工具箱(help stats)
65
矩阵计算(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.0611
0.0222
-0.0194
0.1889
0.0639
0.1056
-0.1028
66
矩阵计算(\)
>> 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
67
矩阵的秩(rank)
>> A = magic(4)
A =
16
2
3
5
11
10
9
7
6
4
14
15
>> rank(A)
ans = 3
13
8
12
1
68
矩阵的行列式(det)
>> A = pascal(4)
A =
1
1
1
1
2
3
1
3
6
1
4
10
>> det(A)
ans =
1
1
4
10
20
69
矩阵的迹(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
70
矩阵的Cholesky分解(chol)
>> A = pascal(4)
A =
1
1
1
1
2
3
1
3
6
1
4
10
>> B = chol(A)
B =
1
1
1
0
1
2
0
0
1
0
0
0
1
4
10
20
1
3
3
1
71
矩阵的特征值(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.9487
-0.0000
-0.4041
0.9091
0.1010
-0.1391
0.9740
-0.1789
D =
1.0000
0
0
0
2.0000
0
0
0
3.0000
72
最优化工具箱(linprog)
什么是线性规划？
min x1  x2
s.t. 2 x1  x2  10
5x1  8x2  40
0  x1  7
4  x2
>>
>>
>>
>>
>>
>>
>>
>>
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
73
最优化工具箱(fminu)
[x,opts]=fminu(f,x0,options,g,p1,p2,…)
求解
情况
最优解
向量
初值
目标函数
字符串
工具箱
函数
算法
参数
help
foptions
目标函数
其他参数
目标函数
梯度
字符串
74
最优化工具箱(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
75
最优化工具箱(fzero)
[x,fv]=fzero(f,x0,options,p1,p2,…)
初值
函数值
零点向量
目标函数
字符串
工具箱
函数
help
foptions
目标函数
其他参数
算法
参数
76
最优化工具箱(fzero)
>> f = inline('x+exp(x)-3');
>> [x,fv]=fzero(f,2)
x =
0.79205996843068
fv =
0
只能求单变量函数的根！
77
最优化工具箱(fmincon)
min f ( x)  3x1  x2  x3
2
s.t. x1  x2  x3  8
x1  2 x2  6
x  x2  x3  16
2
1
2
2
x2 ( x1  x3 )  12
x1 , x2 , x3  0
x = fmincon(’fun’,x0,...
[1 1 1],8,...
[1 2 0],6,...
[0 0 0],[],...
’cons’...
最优解
);
[6 0 0]’
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 = [];
78
微分方程工具箱(ode45)
[t,y] = ode45(odefun,tspan,y0,options,p1,p2...)
求解区间
求解函数
help
odeset
工具箱
函数
初值
解的向量
可用plot(x,y)画出解
目标函数
其他参数
算法
参数
79
微分方程工具箱(ode45)
y '  y  x  1,
y(0)  1, 0  x  2
真解 y( x)  e  x
x
function dydx =f(x,y)
dydx = y – x + 1;
>> [x,y] = ...
ode45(’f’,[0,2],1);
>> plot(x,y,’r-’);
80
多项式计算
 多项式表示：向量
p = [1 0 2] 表示 x 2  2
 多项式求值
polyval(p,1:5)
 多项式求根
r = roots(p)
多项式相乘
p = [1 0 2];
q = [2 1];
s = conv(p,q);
多项式相除
[q,r] = deconv(s,p);
多项式求导
d = polyder(p)
 用根生成多项式
p = poly(r)
81
多项式计算
 多项式拟合(polyfit)
>> x = linspace(0,pi,8)
x =
0 0.4488 0.8976 1.3464 1.7952 2.2440
>> y = sin(x)
y =
0 0.4339 0.7818 0.9749 0.9749 0.7818
>> 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)
2.6928
3.1416
0.4339
0.0000
1.2
1
0.8
0.6
0.4
0.2
0
-0.2
0
0.5
1
1.5
2
2.5
3
82
3.5
插值(interp1)
yi=interp1(x，y，xi，'method')
插值结果
插值节点
被插值点
插值方法
‘nearest’: 最邻近插值
‘linear’: 线性插值
‘spline’: 三次样条插值
‘cubic’: 立方插值
缺省时: 分段线性插值
83
插值(interp2)
z=interp2(x0,y0,z0,x,y,’method’)
被插值点
的函数值
插值
节点
被插值点
插值方法
‘nearest’ 最邻近插值
‘linear’
双线性插值
‘cubic’
双三次插值
缺省时,
双线性插值
84
统计工具箱
 参数估计
betafit, binofit, mle, …
 概率密度函数
betapdf, binopdf, hygepdf, …
 各种分布的均值、方差
betastat, binostat, chi2stat, …
 线性、非线性的各种回归模型
anova1, anova2, polyfit, regress, …
lsqnonneg, nlinfit, nlintool, …
85
matlab上机操作
 遵守机房各项规定
 保持机房卫生整洁
 matlab系统若有问题，请重新开机，在选择
进入系统前按Ctrl-R
86
matlab练习题
 求下面表达式的值

  x4 
1  2.621  0.36

0.85
x

 2 
 x1  x3 
 
y  174.42 
x6 x7
 x2  x2  x1 
 0.56 3 / 2



 x4 
x 
 2
1.16
其中
x  0.1, 0.3, 0.12, 0.09, 1.57, 17.34, 0.87,   0.6174
87
matlab练习题
 已知矩阵A,B如下, 计算并比较
A*B, B*A, A.*B, B.*A,
A/B, A\B, A./B, A.\B
1
 0
A  3.14 2.5  106
  2
3
2
0
1
6, B  2
3
2  3 102
1
 1
1 
 1
88
matlab练习题
 已知矩阵A,B如下, 给出这两个矩阵的
(1)逆 (2)行列式 (3)迹 (4)特征值和特征向量
(5)最简阶梯形矩阵
1
 0
A  3.14 2.5  106
  2
3
2
0
1
6, B  2
3
2  3 102
1
 1
1 
 1
89
matlab练习题
 已知矩阵A及多项式如下, 给出多项式在矩阵上的值；
比较结果矩阵和原矩阵的特征值和特征向量。
1
 0

4
A  3.14 2.5  10
  2
3
2

6,
1 
f ( x)  x  x  1
2
90
matlab综合训练
 矩阵: 写一个程序，对某一输入的对称矩阵A进行
分解, A=LDL’, 其中L为单位下三角矩阵，D为对角
矩阵。 参考《数值计算》中的Cholesky分解。
 画图: 给出某一曲面方程及一垂直于XOY平面的平
面方程，利用平面把曲面割成两片，画出这两片，
中间的缝隙要逐渐变大。
 算法: 输入平面上的至少四个点，挑选当中的部分
点，并把它们连成一个多边形，要求面积尽量大，
给出其面积之值。
91
matlab命令速查
: ... .* ./ .^ \ > >> >= == ~= % abs all and any asin
asinh bar3 ceil cell char chol comet3 conj contour conv
deconv det diag disp doc eig elfun else elseif end eps
eval exp eye false feval find fix floor fmincon fminu
for format fprintf function fzero help hilb i if imag
inf input interp1 interp2 inv linprog linspace log log10
lookfor j magic matfun max mesh meshgrid mod NaN not
num2str ode45 ones optim or pascal pi pie3 plot plot3
polar poly polyfun polyder polyfit polyval prod rand
rank real roots round sin sinh sort sqrt stats str2num
struct surf switch title trace tril triu true view what
which while who whos xlabel xor ylabel zeros
92
作业: 三维画图
在某海域测得一些点(x,y)处的水深z由下
表给出，船的吃水深度为5英尺，在矩形区
域（75，200）*（-50，150）里的哪些地方
船要避免进入。
x
y
z
129.0 140.0 103.5 88.0 185.5 195.0 105.5
7.5 141.5 23.0 147.0 22.5 137.5 85.5
4
8
6
8
6
8
8
x
y
157.5 107.5
-6.5 -81.0
z
9
9
77.0
3.0
81.0
56.5
8
8
162.0 162.0 117.5
-66.5 84.0 -33.5
9
4
9
93
作业: 三维画图
94
THANK YOU
95