Transcript 第二讲MATLAB的数值计算

第二讲 MATLAB的数值计算
—— matlab 具有出色的数值计
算能力，占据世界上数值计算软
件的主导地位
数值运算的功能
创建矩阵
矩阵运算
多项式运算
线性方程组
数值统计
线性插值
函数优化
微分方程的数值解
一、命令行的基本操作
1. 创建矩阵的方法
直接输入法
规则：
 矩阵元素必须用[ ]括住
 矩阵元素必须用逗号或空格分隔
 在[ ]内矩阵的行与行之间必须
用分号分隔
》a=1; b=2; c=3;
》x=[5 b c; a*b a+c c/b]
x=
5.000 2.000 3.000
2.000 4.000 1.500
》y=[2, 4, 5;
3 6 8]
y=
2 4 5
3 6 8
矩阵元素
矩阵元素可以是任何matlab表达式 ，可以
是实数 ，也可以是复数，复数可用特殊函
数I，j 输入。大的矩阵可以用分行输入，
回车键代表分号。
a=[1 2 3;4 5 6]
x=[2 pi/2; sqrt(3) 3+5i]
符号的作用
逗号和分号的作用
逗号和分号可作为指令间
的分隔符，matlab允许多条语句
在同一行出现。
分号如果出现在指令后，
屏幕上将不显示结果。
注意：只要是赋过值的变量，不管是
否在屏幕上显示过，都存储在工作空
间中，以后可随时显示或调用。变量
名尽可能不要重复，否则会覆盖 。
当一个指令或矩阵太长时，可用•••
续行
冒号的作用
用于生成等间隔的向量，默认
间隔为1。
用于选出矩阵指定行、列及元
素。
循环语句
2.用matlab函数创建矩阵
空阵 [ ] — matlab允许输入空阵，当一
项操作无结果时，返回空阵。
rand —— 随机矩阵
eye —— 单位矩阵
zeros ——全部元素都为0的矩阵
ones ——全部元素都为1的矩阵
diag ——产生对角矩阵
例 》eye(2，3)
ans=
1 0 0
0 1 0
》ones(2，3)
ans=
1 1 1
1 1 1
》V=[5 7 2]；A=diag(V)
A=
5 0 0
0 7 0
0 0 2
》zeros(2，3)
ans=
0 0 0
0 0 0
例 》eye(2)
ans=
1 0
0 1
》zeros(2)
ans=
0 0
0 0
》ones(2)
ans=
1 1
1 1
例 在区间[20,50]内均匀分布的5阶随机
矩阵。
命令如下：
x=20+(50-20)*rand(5)
此外，常用的函数还有reshape(A,m,n)，
它在矩阵总元素保持不变的前提下，将
矩阵A重新排成m×n的二维矩阵。
也可用linspace函数产生行向量。其调用
格式为： linspace(a, b, n)
其中a和b是生成向量的第一个和最后一
个元素，n是元素总数。
例 》a=linspace(1 , 10 , 10)
a=
1 2 3 4 5 6 7 8 9 10
还有伴随矩阵、稀疏矩阵、魔方矩阵
(magic)、对角矩阵、范德蒙等矩阵的
创建，就不一一介绍了。
注意：matlab严格区分大小写字母，因
此a与A是两个不同的变量。
matlab函数名必须小写。
3.用m文件创建矩阵
对于比较大且比较复杂的矩阵，可
以为它专门建立一个M文件。下面通
过一个简单例子来说明如何利用M文
件创建矩阵。
例 利用M文件建立MYMAT矩阵。
(1) 启动有关编辑程序或Matlab文本编辑器，
并输入待建矩阵。
(2) 把输入的内容以纯文本方式存盘(设文
件名为mymatrix.m)。
(3) 在Matlab命令窗口中输入mymatrix，即
运行该M文件，就会自动建立一个名为
MYMAT的矩阵，可供以后使用。
4.用冒号表达式创建矩阵
利用冒号表达式可以线性等间距地建立一个
向量来创建矩阵
一般格式是：e1:e2:e3
其中e1为初始值，e2为步长，e3为终止值。
或者为：（start: step: end）
例 》a=[1:2:10]
a=
1 3 5 7 9
5. 矩阵的修改
 直接修改
可用键找到所要修改的矩阵，用键
移动到要修改的矩阵元素上即可修改。
 指令修改
可以用A(,)=  来修改。
例如
a=[1 2 0;3 0 5;7 8 9]
a =1 2 0
3 0 5
还可以用函数subs
7 8 9
a(3,3)=0
修改，matlab6.0还
a =1 2 0
可用find函数修改。
3 0 5
7 8 0
二、数据的保存与获取
把Matlab工作空间中一些有用的数
据长久保存下来的方法是生成mat数
据文件。
save —— 将工作空间中所有的变
量存到matlab.mat文件中。
默认文件名
save data——将工作空间中所有
的变量存到data.mat文件中。
save data a b ——将工作空间中
a和b变量存到data.mat文件中。
下次运行Matlab时即可用load指
令调用已生成的mat文件。
load ——
load data ——
load data a b ——
即可恢复保
存过的所有
变量
mat文件是标准的二进制文件，
还可以ASCII码形式保存。
三、矩阵运算
1. 矩阵加、减（＋,－）运算
规则：
 相加、减的两矩阵必须有相同的行和
列两矩阵对应元素相加减。
 允许参与运算的两矩阵之一是标量。
标量与矩阵的所有元素分别进行加
减操作。
2. 矩阵乘（）运算
规则：
A矩阵的列数必须等于B矩阵的行数
标量可与任何矩阵相乘
a=[1 2 3;4 5 6;7 8 0];b=[1;2;3];c=a*b
c =14
32
23
d=[-1;0;2]; f=pi*d
f = -3.1416
0
6.2832
矩阵除的运算在线性代数中没有，
有矩阵逆的运算，在matlab中有两种
矩阵除运算。
两种除法：\和/，分别表示左除和右除。
如果A矩阵是非奇异方阵，则A\B和B/A
运算可以实现。
A\B等效于A的逆左乘B矩阵，而B/A等效
于A矩阵的逆右乘B矩阵。
对于矩阵来说，左除和右除表示两种不
同的除数矩阵和被除数矩阵的关系。对
于矩阵运算，一般A\B≠B/A。
3. 矩阵乘方—— a^n,a^p,p^a
a ^ p —— a 自乘p次幂
方阵
>1的整数
对于p的其它值,计算将涉及特征值
和特征向量，如果p是矩阵，a是标量
a^p使用特征值和特征向量自乘到p次
幂；如a,p都是矩阵，a^p则无意义。
a=[1,2,3;4,5,6;7,8,9];a^2
ans =30
36
42
66
81
96
102 126 150
※当一个方阵有复数特征值或负实
特征值时，非整数幂是复数阵。
a^0.5
ans =
0.4498 + 0.7623i 0.5526 + 0.2068i 0.6555 -0.3487i
1.0185 + 0.0842i 1.2515 + 0.0228i 1.4844 - 0.0385i
1.5873 - 0.5940i
1.9503 - 0.1611i 2.3134 + 0.2717i
4. 矩阵的其它运算
inv —— 矩阵求逆
det —— 行列式的值
eig —— 矩阵的特征值
diag —— 对角矩阵
’ —— 矩阵转置
sqrt —— 矩阵开方
5. 矩阵的范数
矩阵范数的函数为：
(1) norm(V)或norm(V,2)：计算矩阵V的
2—范数。
(2) norm(V,1)：计算矩阵V的1—范数。
(3) norm(V,inf)：计算矩阵V的∞—范数。
6.矩阵的一些特殊操作
矩阵的变维
a=[1:12];b=reshape(a,3,4)
c=zeros(3,4);c(:)=a(:)
矩阵的变向
rot90:旋转; fliplr:左右翻; flipud:上下翻
矩阵的抽取
diag:抽取主对角线;(对于非方阵的情况?)
tril: 抽取主下三角;
triu:抽取主上三角，然后其余补零元素
矩阵的扩展
关系运算
关系符号
意义
<
<=
>
>=
==
~=
小于
小于或等于
大于
大于或等于
等于
不等于
关系运算符的运算法则
(1) 当两个比较量是标量时，直接比较两数的
大小。若关系成立，关系表达式结果为1，
否则为0。
(2) 当参与比较的量是两个维数相同的矩阵时，
比较是对两矩阵相同位置的元素按标量关
系运算规则逐个进行，并给出元素比较结
果。最终的关系运算的结果是一个维数与
原矩阵相同的矩阵，它的元素由0或1组成。
(3) 当参与比较的一个是标量，而另一个是矩阵
时，则把标量与矩阵的每一个元素按标量关
系运算规则逐个比较，并给出元素比较结果。
最终的关系运算的结果是一个维数与原矩阵
相同的矩阵，它的元素由0或1组成。
注意：其书写方法与数学中的不等式符号不尽
相同。
7. 矩阵的数组运算
数组运算指元素对元素的算术运算，
与通常意义上的由符号表示的线性代数
矩阵运算不同。
1.
数组加减(.+,.-)
a.+b
a.- b
对应元素相加减（与矩阵加
减等效）
2. 数组乘除(，./，.\）
ab —— a，b两数组必须有相同的行
和列两数组相应元素相乘。
a=[1 2 3;4 5 6;7 8 9];
b=[2 4 6;1 3 5;7 9 10];
a.*b
ans =
2
8
18
4
15
30
49
72
90
a=[1 2 3;4 5 6;7 8 9];
b=[2 4 6;1 3 5;7 9 10];
a*b
ans =
25
55
85
37
85
133
46
109
172
a./b=b.\a
—— 给出a,b对应元素间的商.
a.\b=b./a
a./b=b.\a — 都是a的元素被b的对应元
素除， “/”是斜杠
a.\b=b./a — 都是b的元素被a的对应元
素除， “\”是反斜杠
例: a=[1 2 3];b=[4 5 6]; c1=a.\b; c2=b./a
c1 = 4.0000 2.5000 2.0000
c2 = 4.0000 2.5000 2.0000
3. 数组乘方(.^) — 元素对元素的幂
例:
a=[1 2 3];b=[4 5 6];
z=a.^2
z=
1.00
4.00
z=a.^b
z=
1.00
32.00
(1 .^4
2 .^5
9.00
729.00
3 .^6)
四、 多项式运算
matlab语言把多项式表达成一个行向量，
该向量中的元素是按多项式降幂排列的。
f(x)=a0xn+a1xn-1+……an-1x+an
可用行向量 p=[a0 a1 …… an-1 an ]表示
1. poly —— 产生特征多项式系数向量
特征多项式一定是n+1维的
特征多项式第一个元素一定是1
例:a=[1 2 3;4 5 6;7 8 0];
p=poly(a)
p =1.00 -6.00 -72.00 -27.00
p是多项式p(x)=x3-6x2-72x-27的系数
matlab描述方法，我们可用：
p1=poly2str(p,‘x’) — 函数文件，显示
数学多项式的形式
p1 =x^3 - 6 x^2 - 72 x – 27
注意：多项式中缺少的幂次用‘0’补齐。
2.roots —— 求多项式的根
a=[1 2 3;4 5 6;7 8 0];p=poly(a)
p=
1.00 -6.00 -72.00 -27.00
r=roots(p)---------求由p构成的多项式的根
r = 12.12
-5.73 ——显然 r是矩阵a的特征值
-0.39
当然我们可用poly令其返回多项式
形式(这是poly的第二个功能)
p2=poly(r)
p2 =
1.00 -6.00 -72.00 -27.00
matlab规定多项式系数向量用行
向量表示，一组根用列向量表示。
P=poly(r),输入r是多项式所有根，返回值
为代表多项式的行向量形式。
P=poly(A)，输入是N*N的方阵，返回值p
是长度为N+1的行向量多项式，它是矩
阵A的特征多项式，也就是说多项式p的
根是矩阵A的特征值。
求根的另一种方法
str1='x^3-6x^2-72x-27';
p1=str2poly(str1);
r=roots(p1);
注：str2poly 实现把一个字符串表示的多项式转换为
一个行向量表示的多项式。
poly2str 同理。
3.conv多项式乘运算(向量卷积)
例:a(x)=x2+2x+3; b(x)=4x2+5x+6;
c = (x2+2x+3)(4x2+5x+6)
a=[1 2 3];b=[4 5 6];
c=conv(a,b)或c=conv([1 2 3],[4 5 6])
c = 4.00 13.00 28.00 27.00 18.00
p=poly2str(c,‘x’) 其中x表示自变量
p = 4 x^4 + 13 x^3 + 28 x^2 + 27 x + 18
4.deconv多项式除运算(解卷积)
a=[1 2 3];
c = [4.00 13.00 28.00 27.00 18.00]
d=deconv(c,a)
d =4.00
5.00
6.00
[d,r]=deconv(c,a)
余数
c除a后的整数
它们之间的关系为: c = conv(a,d)+r
5.多项式导数或微分
matlab提供polyder函数计算多项式的导数。
命令格式：
polyder(p): 求p的导数
polyder(a,b): 求多项式a,b乘积的导数
[p,q]=polyder(a,b): 求多项式a除以b的商的
导数，并以p/q的格式表示。
例：a=[1 2 3 4 5]; poly2str(a,'x')
ans = x^4 + 2 x^3 + 3 x^2 + 4 x + 5
b=polyder(a)
b=4 6 6 4
poly2str(b,'x')
ans =4 x^3 + 6 x^2 + 6 x + 4
6.多项式的积分
matlab提供polyint函数计算多项式的积分。
命令格式：
polyint(p,k): 求多项式p的积分，设积分的
常数项为k， polyint(p) 默认k=0
例：a=[1 2 3 4 5]; poly2str(a,'x')
ans = x^4 + 2 x^3 + 3 x^2 + 4 x + 5
b=polyint(a,8)
b = 0.2 0.5 1.0 2.0 5.0 8.0
poly2str(b,'x')
ans =0.2 x^5 + 0.5 x^4 + x^3 + 2 x^2 + 5 x + 8
五、代数方程组求解
matlab中有两种除运算左除和右除。
对于方程ax=b，a 为an×m矩阵，有三种情况：
 当n=m时，此方程成为“恰定”方程
 当n>m时，此方程成为“超定”方程
 当n<m时，此方程成为“欠定”方程
matlab定义的除运算可以很方便地解上
述三种方程
1.恰定方程组的解
方程ax=b (a为非奇异)
x=a-1 b
矩阵逆
两种解:
x=inv(a)b — 采用求逆运算解方程
x=a\b — 采用左除运算解方程
注：若a为奇异的，则Matlab适当给出
警告信息或者给出结果为inf。
例:
x1+2x2=8
2x1+3x2=13
1 2 x1
8
=
2 3 x2
13
方程ax=b
a=[1 2;2 3];b=[8;13];
x=inv(a)*b
x=
2.00
3.00
a
x = b
 x=a\b
x=
2.00
3.00
2.超定方程组的解
方程的个数大于未知量个数时，方程一般无解。
方程解 (a ' a)x=a ' b
x=(a‘ a)-1 a ’ b —— 求逆法 (也用到了最小二
乘解的原理)
x=a\b —— matlab用最小二乘法找一
个准确地基本解。
例:
x1+2x2=1
1 2
2x1+3x2=2
2 3
3x1+4x2=3
3 4
a
x1
x2
1
=2
3
x = b
a=[1 2;2 3;3 4];b=[1;2;3];
解1 x=a\b
解2 x=inv(a'a)  a'  b
x=
x=
1.00
1.00
0
0.00
3.欠定方程组的解
当方程数少于未知量个数时,即不定
情况,有无穷多个解存在。
matlab可求出两个解：
用除法求的解x是具有最多零元素的
解
是具有最小长度或范数的解，这个
解是基于伪逆pinv求得的。
x1
1 2 3
1
x2 =
2 3 4
2
x3
x1+2x2+3x3=1
2x1+3x2+4x3=2
a=[1 2 3;2 3 4];b=[1;2];
a
x=a\b
x=pinv(a)b
x=
x=
1.00
0.83
0
0.33
0
-0.17
x = b
六、微分方程求解
微分方程求解的仿真算法有多种，常用
的有Euler（欧拉法）、Runge Kutta(龙
格-库塔法。
Euler法称一步法，用于一阶微分方程
dy
 f ( x, y ), y ( x0 )  y0
dx
当给定仿真步长时：
yn 1  yn
dy yn 1  yn


dx xn 1  xn
h
所以
yn+1 = yn + h·f (xn,yn)
y(x0)=y0
n=0,1,2…
Runge Kutta法
yn 1
1
1
 yn 
k1 
k2
2
2
k1=hf(xn,yn)
k2=hf(xn+h,yn+k1)
龙格-库塔法：实际上取两点斜率
的平均 斜率来计算的，其精度高
于欧拉算法 。
龙格-库塔法：ode23 ode45
··
·
2
例：x+(x -1)x+x=0
·
为方便令x1=x，x2=x分别对x1,x2求一
阶导数，整理后写成一阶微分方程组
形式
·
x1=x2
·
x2=x2(1-x12)-x1
1. 建立m文件
2. 解微分方程
建立m文件
function xdot=wf(t,x)
xdot=zeros(2,1)
xdot(1)=x(2)
xdot(2)=x(2)*(1-x(1)^2)-x(1)
给定区间、初始值;求解微分方程
t0=0; tf=20; x0=[0 0.25]';
[t,x]=ode23('wf', t0, tf, x0)
plot(t,x), figure(2),plot(x(:,1),x(:,2))
命令格式:
[T,Y] = ODE23(ODEFUN,TSPAN,Y0)
建立m文件
function dxdt=wf(t,x)
dxdt=[x(2);x(2)*(1-x(1)^2)-x(1)];
求解微分方程
[t,x]=ode23(@wf,[0 30],[0 0.25]);
plot(t,x);
figure(2)
plot(x(:,1),x(:,2))
3
2
1
0
-1
-2
-3
0
5
10
15
20
25
30
3
2
1
0
-1
-2
-3
-2.5
-2
-1.5
-1
-0.5
0
0.5
1
1.5
2
2.5
七、函数优化
寻优函数：
fmin —— 单变量函数
fmins —— 多变量函数
constr —— 有约束条件
无约束条件
例1：f(x)=‘x2+3x+2’在[-5 5]区间的最
小值
f=fmin('x^2+3*x+2',-5,5)
例2：f(x)=100(x2-x12)2+(a-x1)2在x1=a,
x2=a2处有最小值
function f=xun(x,a)
f=100*(x(2)-x(1).^2).^2+(a-x(1)).^2;
x=fmins('xun',[0,0],[ ],[ ],sqrt(2))
八、数据分析
数据分析相关的函数位于目录：
toolboxs\matlab\datafun下
Matlab对矩阵操作的规定：如果是向量，则对数据整
体操作；如果是矩阵，则对矩阵的列操作。
max —— 各列最大值
mean —— 各列平均值
sum —— 各列求和
std —— 各列标准差
var —— 各列方差
sort —— 各列递增排序
cumsum —— 元素累计和
cumprod —— 元素累计积
八、数据分析
数据分析相关的函数位于目录：
toolboxs\matlab\datafun下。
Matlab对矩阵操作的规定：如果是向量，则对数据
整体操作；如果是矩阵，则对矩阵的列操作。
max —— 各列最大值
mean —— 各列平均值
sum —— 各列求和
std —— 各列标准差
var —— 各列方差
sort —— 各列递增排序
cumsum —— 元素累计和
cumprod —— 元素累计积
例：x=[1 3 2 4],y=[1 2 3 8;5 6 7 4]
sort(x),sort(y),max(y)
>>1 2 3 4
>>1 2 3 4
5 6 7 8
>>5 6 7 8
九、拟合与插值
1. 多项式拟合
采用最小二乘法对给定的数据进行多项式
拟合，最后给出多项式的系数。
p=polyfit(x,y,n)，采用n次多项式p来拟合
数据x和y，从而使得y与p(x)最小均方
差最小。
x0=0:0.1:1;
y0=[-.447 1.978 3.11 5.25 5.02 4.66 4.01 4.58
3.45 5.35 9.22];
p=polyfit(x0,y0,3)
p = 56.6915 -87.1174 40.0070 -0.9043
xx=0:0.01:1;yy=polyval(p,xx);
plot(xx,yy,'-b',x0,y0,'or')
曲线拟合图形用户接口
Matlab7.0提供了支持曲线拟合的图形用户接口。
在“Figure”窗口“Tools\Basic Fitting”菜单
中。
为了使用该工具，先用待拟合的数据画图。
x=0:0.2:10;
y=0.25*x+20*sin(x);
plot(x,y,'ro');
在复选框“Plot fits”中选择“cubic”。
2.插值
插值的定义——是对某些集合给定的数据点
之间函数的估值方法。
当不能很快地求出所需中间点的函数时，插
值是一个非常有价值的工具，它可以在已知
数据之间寻找估计值，常用到信号处理和图
像处理中。
Matlab提供了一维、二维、 三次样条等许多
插值选择。
interp1——一维插值
interp2——二维插值
interp3——三维插值
spline——三次样条插值
griddata ——栅格数据插值
 利用已知点确定未知点
 粗糙—— 精确
 集合大的—— 简化的
一维插值就是对一维函数y=f(x)进行插值。
yi=interp1(x,y,xi,method),x必须是向量，y可以
是向量也可以是矩阵。(x,y)代表的是已知数据。
这时，xi代表需要估计值的位置，yi表示插值
后的估计值。method用于指定插值的方法：
1.method='nearest'，在已知数据的最临近点设置插
值点，对插值点的数进行四舍五入。对超出范围的点
返回一个NaN。此方法是最快的插值方法，但数据平
滑方面最差，其得到的数据是不连续的。
2.method='linear'，采用直线连接相邻的两
点，即线性插值，是此函数的缺省默认方
法。执行速度比最临近插值稍慢，数据平
滑要由于临近插值，且数据是连续的。
3. method='spline'，采用三次样条函数来
获得插值点。处理速度最慢，可以产生最
光滑的结果。Matlab提供了一个样条插值
工具箱，位于 toolbox\splines下。
4. method='pchip'，采用分段三次埃尔米
特多项式插值。
例：x=0:2*pi;y=sin(x);xi=0:0.1:8;
yi1=interp1(x,y,xi, 'linear')
yi2=interp1(x,y,xi, 'nearest')
yi3=interp1(x,y,xi, 'spline')
yi4=interp1(x,y,xi, 'cubic')
p=polyfit(x,y,3);yy=polyval(p,xi);
subplot(3,2,1);plot(x,y,'o');
subplot(3,2,2);plot(x,y,'o',xi,yy);
subplot(3,2,3);plot(x,y,'o',xi,yi1);
subplot(3,2,4);plot(x,y,'o',xi,yi2);
subplot(3,2,5);plot(x,y,'o',xi,yi3);
subplot(3,2,6);plot(x,y,'o',xi,yi4);
二维插值
二维插值主要应用于图像处理和数据的可视
化，对双变量的函数z=f(x,y)进行插值。
zi=interp2(x,y,z,xi,yi,method),原始数据x,y,z
决定插值函数z=f(x,y),返回值zi是(xi,yi)在函
数f(x,y)上的值。
method同样可以采用最临近插值、双线性
插值、三次样条插值。
小
结
本节介绍了matlab语言的数值运算
功能，通过学习应该掌握：
 如何创建矩阵、修改矩阵
 符号的用法
 矩阵及数组运算
 多项式运算
 线性方程组与微分运算