数值计算中 国 民 航 大 学 电 子 信 息 工 程 学 院本讲提纲
与线性代数有关的 matlab指令
解线性方程组
多项式
数据分析函数
泛函指令
信号处理
系统分析中 国 民 航 大 学 电 子 信 息 工 程 学 院与线性代数有关的 matlab指令
三角分解( LU分解) [L,U,P]=lu (A)
说明:矩阵 A的 LU分解,使 LU=PA
行列式 det(A)
逆运算 inv(A)
广义逆 pinv(A)
零空间和值空间的正交基矢量指令,null(A) orth(A)
范数和条件数指令,norm(A,flag) cond(A)
中 国 民 航 大 学 电 子 信 息 工 程 学 院与线性代数有关的 matlab指令
秩 r=rank(A,tol)
奇异值分解 [U,S,V]=svd(A)
特征值分解
d=eig(A) [V,D]=eig(A)
矩阵 A的特征值的条件数
condeig(A)
例 A=randn(2,2)
[v,d]=eig(A)
condeig(A)
r=rank(A)
中 国 民 航 大 学 电 子 信 息 工 程 学 院与线性代数有关的 matlab指令
矩阵函数
A^p,p^A,expm(A),logm(A),sqrtm(A)
funm(A,’FN’)
例,A=(reshape(1:9,3,3))’;
expm(A);
logm(A);
sqrtm(A);
B=funm(A,’sin’)
中 国 民 航 大 学 电 子 信 息 工 程 学 院解线性方程组
线性方程组 Ax=b的解,由左除运算符求解 x=A\b
线性方程组 xA=b的解,由右除运算符求解 x=b/A
线性方程组可以是正定方程,欠定方程,过定方程
和求逆法相比,运算符方法速度快,精度高( A的条件数大)
A奇异时,采用 pinv(A)*b求解
线性方程组 y=Ax+n的解,由
x=inv(A’*A)*A’*y
x=A\y
x=pinv(A)*y
求解中 国 民 航 大 学 电 子 信 息 工 程 学 院多项式
多项式由降幂排列的系数向量表示
多项式运算函数
P=conv(p1,p2)
[q,r]=deconv(p1,p2)
R=root(p)
P=poly(A) 求矩阵 A的特征多项式,或求向量指定根所对应的多项式
P=polyfit(x,y,n)
pA=polyval(p,s)
pM=polyvalm(p,s)
[r,p,k]=residue(b,a) r 留数,p 极点,k直项
poly2str(P,’ s’ )
real(p)
中 国 民 航 大 学 电 子 信 息 工 程 学 院多项式多项式的值和多项式的根
求( s2+2)(s+4)(s+1)/(s3+s+1)的商及余多项式
求 x=5,x=6和 x=[1 2; 3 4]时多项式
x3-6x2-72x-27的值中 国 民 航 大 学 电 子 信 息 工 程 学 院多项式
数据插值数据插值是利用一批已知的测量数据,采用某种算法平滑地估算出测量点之间的数据。 Matlab提供一维插值函数 interp1(x,y,xi,’parameter’)
二维插值函数 interp2(x,y,z,xi,yi,’parameter’)
三维插值函数 interp3
nearest 最近插值,linear 线性插值
spline 三次样条插值 cubic 三次插值例:利用正弦曲线上的 9个样点来说明一维插值函数的使用方法中 国 民 航 大 学 电 子 信 息 工 程 学 院多项式和卷积
数据拟合数据拟合是利用一批已知测量点上的取值,按照某个确定的准则寻找一条可用函数表示的平滑曲线,以使该函数在已知测量点上尽可能接近测量点的取值。
例,用 6阶多项式在区间 [0,2.5]上拟合误差函数
dtexe rf x t 0 22)(?
中 国 民 航 大 学 电 子 信 息 工 程 学 院数据分析函数
Matlab的数据分析函数有三条约定
用矩阵表示数据时,行代表一次观测值,列代表不同的观测变量
若函数的输入量是向量,则运算对整个向量进行
若函数的输入量是矩阵,则运算按列进行
统计分析函数
mean(x),mean(x,dim)
min(x),[y,I]=min(x),[y,I]=min(x,[],dim)
max(x),[y,I]=max (x),[y,I]=max (x,[],dim)
median(x),median(x,dim)
中 国 民 航 大 学 电 子 信 息 工 程 学 院数据分析函数
统计分析函数
sort(x),[y,I]=sort(x),sort(x,dim)
std(x),std(x,1),std(x,flag,dim)
var(x),var(x,1),var(x,w)
A = [-1 1 2 ; -2 3 1 ; 4 0 3]
cov(x),cov(x,y)
corrcoef(x),corrcoef(x,y)
中 国 民 航 大 学 电 子 信 息 工 程 学 院数据分析函数例 逆序输出向量 x中最大的 5个元素;保留 x中最大的五个元素,其他元素全部置换为 0
x=[1;2;5;3;6;8;7;2;9;7;4;0]; x’
[y,j]=sort(x); xm=y(end:-1:end-4);
x(j(1:length(x)-5))=0;x’
中 国 民 航 大 学 电 子 信 息 工 程 学 院数据分析函数
差分和累积指令
diff(x),diff(x,m),diff(x,m,dim) 微分
gradient(F),[Fx,Fy]=gradient(F),梯度
del2(U,hx,hy) 5点 Laplacian
prod(x,n) 乘积
sum(x,n) 和
trapz(x,Y,n) 梯形法求积分
cumprod(x,n) 累积
cumsum(x,n) 累和
cumtrapz(x,y,n) 累计积分中 国 民 航 大 学 电 子 信 息 工 程 学 院
Matlab泛函指令
在 Matlab中,凡以“函数”为输入量的指令,统称为泛函指令
利用泛函指令可以解决下列问题
求函数零点
求函数极值点
数值积分
解常微分方程中 国 民 航 大 学 电 子 信 息 工 程 学 院求函数零点
Matlab 指令
fzero 单变量求零点语法,x=fzero(fun,x0)
x=fzero(fun,x0,options)
x=fzero(fun,x0,options,P1,P2)
[x,fval]=fzero(… )
[x,fval,exitflag]=fzero(… )
中 国 民 航 大 学 电 子 信 息 工 程 学 院求函数零点
参数说明
fun:函数名字符串表达式 例 ‘ y=2x+3’
内联函数 例 inline(‘ y=2x+3’,’ x’ )
M函数文件 例 Hf=@fun
x0:零点的初始猜测,可以为标量或向量。
标量:在 x0 点两侧寻找与 x0最靠近的零点向量:在 [a,b]寻找零点。 [a,b]中最多有 1个零点
Option,用来设置是否显示中间结果,设置以自变量或函数值迭代误差控制的迭代条件,通过 optimset函数来设置参数为,
display
off 不输出
Iter 每一迭代步骤要输出显示
Final 输出显示最终结果中 国 民 航 大 学 电 子 信 息 工 程 学 院求函数零点
tolx,用来给定的 x的收敛允许值
[]:没有选项设置例,fzero(fun,xo,optimset(‘display’,’iter’))
Exitflag 描述函数 fzero的退出条件
exitflag>0:找到零点
exitflag<0:没找到零点
P1,P2:传递给函数的附加参数中 国 民 航 大 学 电 子 信 息 工 程 学 院求函数零点
函数说明
寻找单变量函数值为 0时的自变量的值
根据函数是否穿越横轴来决定零点的
不能确定函数曲线仅触及横轴和不穿越的零点。例 |sin(x)|
中 国 民 航 大 学 电 子 信 息 工 程 学 院求函数零点
解题步骤:
利用 Matlab作图指令获得初步近似解确定零点可能存在的自变量区间
plot指令绘图观察零点位置,或用 ginput指令获得更精确的零点
利用 fzero指令求精确解中 国 民 航 大 学 电 子 信 息 工 程 学 院求函数零点例 1
求函数 f(t)=(sin2t)e-at-b|t|的零点思路,在 matlab中表示函数观察函数的零点分布利用零点指令求零点中 国 民 航 大 学 电 子 信 息 工 程 学 院求函数极值点
fminbnd 求单变量函数的最小值语法,x=fminbnd(fun,x1,x2)
x=fminbnd(fun,x1,x2,options)
x=fminbnd(fun,x1,x2,options,P1,P2)
[x,fval]=fminbnd(?)
[x,fval,exitflag]=fminbnd(?)
[x,fval,exitflag,output]=fminbnd(?)
fminsearch 求多变量函数的最小值语法类似,不同在于,调用 fminbnd时给定求极值的区间,调用 fminsearch时,要给定初始极值向量,
fminsearch函数将计算出所给向量附近的一个极小值中 国 民 航 大 学 电 子 信 息 工 程 学 院求函数极值点例 2 求函数 f(x)=x^3-2x+5在区间 (0,2)上的极小值。
例 3 求函数 f(x,y)=100(y-x^2)^2+(1-x)^2
的极小值点。
中 国 民 航 大 学 电 子 信 息 工 程 学 院数值积分
积分算法开型闭型
Matlab提供闭型数值积分指令
quad 采用递推自适应法计算积分语法,x=quad(fun,x1,x2)
x=quad(fun,x1,x2,tol)
x=quad(fun,x1,x2,tol)
x=quad(fun,x1,x2,tol,trace,P1,P2)
进行反复的积分求解迭代,直到计算得到的相对误差小于 tol时迭代停止。当 trace 取非零时,就会将积分过程显示出来,
允许将系数 p1,p2 直接传递给指定的函数 fun
是否需要计算积分区间端点处的函数值中 国 民 航 大 学 电 子 信 息 工 程 学 院数值积分
参数说明:
X1,x2指定的积分区间
tol 可选项,指定迭代误差
trace 可选项,当不为零时,将会绘制出所采用的积分点的图形例 4 求函数 y=e^(-x^2)在区间 [0,1]上的积分中 国 民 航 大 学 电 子 信 息 工 程 学 院解常微分方程
Matlab提供的解常微分方程指令
微分方程解算指令
被解算指令调用的 ODE文件格式
积分算法参数选项 options处理指令
输出指令常用指令 ode45
指令格式,ode45(fun,tspan,y)
tspan:时间区间
y:初始值中 国 民 航 大 学 电 子 信 息 工 程 学 院解常微分方程利用 ode 指令求解常微分方程步骤,
把微分方程改写成微分方程组例,y’’’ -3y’’ -y’ y=0,y(0)=0,
y’ (0)=1,y’’ (0)=-1
令 y1=y,y2=y’,y3=y’’ 则
y’’’ -3y’’ -y’ y=0 可化为
y3’ =3*y3+y2*y1
写成方程组形式 y1’=y2
y2’=y3
y3’=3*y3+y2*y1
y1(0)=0
y2(0)=1
y3(0)=-1
中 国 民 航 大 学 电 子 信 息 工 程 学 院解常微分方程
把微分方程编写成 m文件
function dy=f(t,y)
在 m文件中,虽然写微分方程时并不同时包含参数 t和 y,但 m文件的第一行必须包含这两个参数中 国 民 航 大 学 电 子 信 息 工 程 学 院解常微分方程
利用 ode45指令求解例 5求微分方程 x’’-u(1-x^2)x’+x=0在初始条件
x(0)=1,x’(0)= 0情况下的解并图示思路,把上述方程转化成方程组的形式令 y1=x,y2=x’ 则有
y1’=y2
y2’=u(1-y1^2)*y2-y1
y1(0)=0
y2(0)=-1
中 国 民 航 大 学 电 子 信 息 工 程 学 院信号处理
快速 fft变换和逆变换
数字滤波中 国 民 航 大 学 电 子 信 息 工 程 学 院快速 FFT和逆 FFT
快速 fft及其逆 fft表达式
教科书定义
Matlab表达式
1
0
( ) D F T [ x ( n ) ] ( ) (0 1 )N knN
n
X k x n W n N?

1
0
1( ) I D F T [ ( ) ] ( ) (0 1 )N kn
Nkx n X k X k W k NN


1
( ) D F T [ x ( n ) ] ( ) ( 1 )N knN
n
X k x n W n N

1
1( ) I D F T [ ( ) ] ( ) ( 1 )N kn
N
k
x n X k X k W k NN?

2j N
NWe

中 国 民 航 大 学 电 子 信 息 工 程 学 院快速 FFT和逆 FFT
Matlab 指令
fft 一维快速傅里叶变换语法,y=fft(x)
返回 x的离散傅里叶变换,若 x为矩阵,返回每列的离散傅里叶变换
y=fft(x,n)
返回 n点傅里叶变换,n用来指定变换的点数若 size(x,1)<N,补零若 size(x,1)>N,截短
Ifft一维逆快速傅里叶变换语法结构同上中 国 民 航 大 学 电 子 信 息 工 程 学 院快速 FFT和逆 FFT
利用 fft和 ifft计算卷积
)(*)()( nhnxny?
FFT
FFT
序列相乘 IFFT
x(n)
h(n)
X(k)
H(k)
X(k)H(k)
y(n)
中 国 民 航 大 学 电 子 信 息 工 程 学 院快速 FFT和逆 FFT
求 a(n)=1(n=3,4,?,12),n=else 时,a(n)=0
b(n)=1(n=2,3,?,9),n=else 时,b(n)=0的卷积
FFT在信号分析中的应用例试从受噪声污染的信号
y(t)=3sin5t-6cos9t+N(t)中鉴别有用信号,
其中,N(t)-N(0,5)
中 国 民 航 大 学 电 子 信 息 工 程 学 院数字滤波
滤波器传递函数
Matlab指令获得滤波器传递函数
butter 巴特沃斯滤波器语法 [B,A]=butter(n,wn)
参数说明,n:滤波器的阶数
wn:滤波器截止频率归一化频率真实的滤波器截止频率 /Nyquest频率
12
1 2 3 1
12
1 2 3 1
()()
()
nb
nb
na
na
b b z b z b zBzFz
A z a a z a z a z




中 国 民 航 大 学 电 子 信 息 工 程 学 院数字滤波
Cheby1 契比雪夫滤波器
[B,A] = cheby(N,R,wn)
R:通带内的波纹大小
Cheby2 契比雪夫滤波器
[B,A] = cheby(N,R,wn)
R:通带内的波纹大小利用 IIR滤波器对数据进行滤波指令
filter利用 IIR滤波器对数据进行滤波语法,y=filter(B,A,x)
以向量 B作为分子,A作为分母的滤波器对数据
x进行滤波中 国 民 航 大 学 电 子 信 息 工 程 学 院数字滤波例 设计一个低通滤波器,从受噪声干扰的多频率混合信号 x(t)中获取 10hz的信号,其中:
x(t)=sin(2*pi*10t)+cos(2*pi*100t)+n(t)
n(t)-N(0,0.2)
中 国 民 航 大 学 电 子 信 息 工 程 学 院系统分析线性时不变 (LTI)对象指令
状态方程
S_ss=ss(A,B,C,D) 利用状态方程四对组生成 LTI
对象
[A,B,C,D]=ssdata(S_ss) 从 lTI对象获取状态方程四对组
零极点增益
,x A x B u y C x D u
12
1 2 3
( ) ( ) ( )()
( ) ( ) ( )
ms z s z s zG s K
s p s p s p


中 国 民 航 大 学 电 子 信 息 工 程 学 院系统分析
s_zpk=zpk(Z,P,K)利用零极点增益获取三对组生成
LTI对象
[Z,P,K]=zpkdata(s_zpk) LTI对象获取零极点增益三对组
传递函数
S_tf=tf(num,den)利用传递函数二对组生成 LTI对象
[num,den]=tfdata(S_tf) LTI对象生成传递函数二对组
1
1 2 1
1
1 2 1
()
mm
mm
nn
nn
N S N S N S NGs
D S D S D S D


中 国 民 航 大 学 电 子 信 息 工 程 学 院系统分析三者关系优先级状态方程 〉 零极点 〉 传递函数
LTI对象传递函数状态方程 零极点增益中 国 民 航 大 学 电 子 信 息 工 程 学 院系统分析求方框图所描述的多环系统的传递函数:
已知此方程图中 1 1 1
1,2,3
10 1 2 4 4
11
4,1,2 2,3 1
62
s
G G G
s s s s
ss
G H H H
ss





G2G1 G4G3
H1
H2/G4
H3
-- +
-
中 国 民 航 大 学 电 子 信 息 工 程 学 院系统分析
根据已知条件创建各模块的 LTI对象
G1=tf(1,[1 10]);G2=tf(1,[1 1]);
G3=tf([1 1],[1 4 4]); G4=tf([1 1],[1 6]);
H1=tf([1 1],[1 2]);H2=2;H3=1;
计算最内环的传递函数
P1=G3*G4/(1-H1*G3*G4);
P1=minreal(P1);
计算次内环的传递函数
P2=minreal(G2*P1/(1+G2*P1*H2/G4));
中 国 民 航 大 学 电 子 信 息 工 程 学 院系统分析
计算总的传递函数
P3=minreal(G1*P2/(1+G1*P2*H3))
LTI对象的优点
直接按方块图算法进行运算
指令直观,简洁
与传统教科书上一致