第二讲MATLAB的数值计算
Download
Report
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. 数组乘除(,./,.\)
ab —— 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语言的数值运算
功能,通过学习应该掌握:
如何创建矩阵、修改矩阵
符号的用法
矩阵及数组运算
多项式运算
线性方程组与微分运算