MATLAB第4章数值运算基础.ppt

物理与电气工程学院 1 第4章数值运算基础 MATLAB的计算包括数值计算和符号计算两类 本章将带大家学习数值计算部分 其中将主要学习与我们专业密切相关的多项式 方程组求解 数据分析和数字信号处理的快速傅里叶变换 物理与电气工程学院 2 第1节多项式polynomialMATLAB用行向量表示多项式 将多项式的系数按降幂次序存放在行向量中 如 的系数行向量为 注意 多项式中缺少的幂次系数一定要用 0 补齐 物理与电气工程学院 3 一 创建多项式1 系数矢量直接输入法在命令窗口直接输入多式的系数向量 例4 1 输入系数矢量 创建多项式x 3 4 x 2 3 x 2 p 1 432 poly2sym p 将矢量P表示为多项式的手写形式 Polynomialcoefficientvectortosymbolicpolynominal 4 2 方阵特征多项式输入法p poly A 若A为n n的矩阵 则返回值P将是一个含有n 1个元素的行向量 也就是该矩阵特征多项式的系数 例4 2 求矩阵 123 456 780 的特征多项式系数 并转换为多项式形式 a 123 456 780 p poly a poly2sym p 将矢量P表示为多项式的手写形式d1 roots p 由特征多项式求得的特征值d2 eig a 由特征值函数求得的特征值 5 3 根矢量创建p poly A A为待求多项式的根矢量 则返回值将是对应多项式的系数行矢量 该多项式的根为矢量A 此时p poly A 与A roots p 互逆 系统定义P0 1 A x1x2x3 p 1p1p2p3 6 例4 3 根据根矢量 0 5 0 3 0 4i 0 3 0 4i 创建多项式r 0 5 0 3 0 4i 0 3 0 4i p poly r pr real p ppr poly2sym pr 二 多项式运算1 求多项式的值MATLAB的多项式求值方式有两种 按数组运算规则和按矩阵运算规则计算多项式的值 y polyval p x 按数组规则运算 计算多项式在x处的值 p是多项式的系数矢量 x是指定自变量的值 可以是标量 向量或矩阵 如果x是向量或矩阵 则函数的返回值是与x对应的同维向量或矩阵 物理与电气工程学院 7 例4 4 求多项式3x 2 2x 1在5 7和9处的值 p 321 polyval p 579 y polyvalm p x 将矩阵整体 而不是矩阵元素 作为自变量进行计算 p是多项式的系数向量 相当于用矩阵x代替多项式的变量对矩阵而不是对数组进行计算 x必须是方阵 例4 5 求多项式3x 2 2x 1对于矩阵 25 79 的值p 321 polyvalm p 25 79 物理与电气工程学院 8 2 求多项式的根格式 C roots p p为多项式的系数矢量 C为函数返回多项式的根矢量如果C为复数 则必成对出现 例4 6 分别用两种方法求多项式x 5 5x 4 3x 3 6x 2 4x 10的根 a 1 53 64 10 r roots a 物理与电气工程学院 9 3 多项式的乘除运算多项式的乘法conv格式 c conv a b 多项式的乘法运算 也是矢量的卷积运算向量a长度为m 向量b长度为n a和b的卷积定义为 运算结果矢量c为长度k m n 1 例4 7 计算两多项式x 4 5x 3 3x 2 4x 2和x 3 2x 2 5x 3的乘法a 1 53 42 b 12 53 c conv a b poly2sym c 10 多项式的除法deconv格式 q r deconv c a 多项式的除法运算 也是矢量的解卷积运算过程 向量a对向量c进行解卷积得到 商 向量q和 余 向量r q和r分别是商多项式和余多项式 c和a分别是被除多项式和除多项式 使得 c conv a q r 例4 8 计算例4 7中求得的乘积被x 3 2x 2 5x 3除所得结果c 1 3 1230 3633 226 b 12 53 q r deconv c b 物理与电气工程学院 11 4 多项式微积分polyder p 返回多项式系数向量p的导数 例4 9 计算多项式3x 4 5x 3 2x 2 6x 10的微分 p 3 52 610 dp polyder p poly2sym dp polyint p 返回多项式系数p的积分 例4 10 计算多项式12x 3 15x 2 4x 6的积分 p 12 154 6 ip polyint p poly2sym ip 物理与电气工程学院 12 5 多项式的部分分式展开MATLAB提供了residue命令来执行部分分式展开或多项式系数之间的转换 该命令通常用于信号与控制领域中 格式如下 r p k residue b a 该命令是求多项式之比b s a s 的部分分式展开 返回留数r 极点p和直项向量k a和b分别是分母和分子多项式的系数向量 b a residue r p k 上一条语句的逆过程 只是分母多项式系数以归一化形式 最高次项系数a 1 为1 13 如果分母多项式a s 不含重根 则两个多项式可写成如下形式 其中 pi称为极点 ri称为留数 k s 是直项 如果b的次数低于a的次数 则直项为空 如果分母多项式a s 含m重根p j p j m 1 则这m项应展开为 物理与电气工程学院 14 b x 5x 3 3x 2 2x 7 例3 11 两多项式的比为 a x 4x 3 8x 3求部分分式展开 a 4083 b 53 27 r p k residue b a b1 a1 residue r p k 分母最高次项归1 r2 p2 k2 residue 11 1 21 出现重根 笔算结果 物理与电气工程学院 15 6 多项式拟合对于给定的一组数据 xi yi i 1 2 n 如果要采用多项式模型对数据组进行描述 形成如多项式y x f x p p1xn p2xn 1 p3xn 2 pn 1的形式 求取参数p使得量值 2 p 的值最小的过程 称为对数据组进行多项式拟合 其中 MATLAB系统设计polyfit函数采用最小二乘法原理对给定的数组 xi yi i 1 2 n 进行多项式的曲线拟合 最后给出拟合的多项式系数 物理与电气工程学院 16 p polyfit x y n 其中 x y分别表示横 纵坐标向量 n是给定的拟合多项式的最高阶数 返回一个多项式系数向量p 如n 3 若p 10 512 则y 1 x 3 0 5 x 2 1 x 1 2 x 0 p S polyfit x y n 返回n阶多项式系数p S为误差估计结构 p S mu polyfit x y n 对归一化以后的数据进行多项式拟合 用x x u1 u2替代x 其中mu u1u2 u1 mean x u2 std x 物理与电气工程学院 17 例3 13 求误差函数的6阶拟合多项式 x 0 0 1 2 5 生成0至2 5间隔为0 1的自变量y erf x 计算误差函数p polyfit x y 6 求6阶拟合多项式x 0 0 1 5 生成0至5间隔为0 1的自变量y1 erf x 计算误差函数f polyval p x 计算拟合函数的值plot x y1 o x f 绘图函数 p0 s0 mu0 polyfit x y 6 x x mean x std x p1 s1 mu1 polyfit x mu 1 mu 2 y 6 可以看出 拟合区间 02 5 内拟合曲线拟近原曲线 而在区间以外的曲线误差较大 第2节线性代数给定两个矩阵A和B 求X的解 使得 AX BXA B 在MATLAB中 求解线性方程组时 主要采用前面章节介绍的除法运算符 和 求解 X A B是AX B的解X B A是XA B的解对于方程AX B 要求矩阵A和B有相同的行数 X和B有相同的列数 且它们的行数等于A的列数 物理与电气工程学院 19 根据矩阵A的结构 m n 可以将方程分为以下3类 m n方阵系统 可偿试求精确解m n超定系统 可偿试求最小二乘解m n欠定系统 可偿试求少m个非零解 反斜线运算因子 对于不同结构矩阵A 将采用不同的算法处理 MATLAB会自动检测参数特性 选择运用不同的算法解方程组 20 一 方阵系统由n个未知数构成n个方程 方程有唯一的一组解 其一般形式写成如下形式 Ax b其中A是方阵 b是列向量格式 x A b A为非奇异 且条件数不大 x精确 A为非奇异 且条件数很大 注意解的可靠性 A为奇异 解不存在或存在但不唯一 A为接近奇阵或病态矩阵Matrixisclosetosingularorbadlyscaled Resultsmaybeinaccurate 尽量使用格式 x A b 而不用格式 x inv A b 物理与电气工程学院 21 例4 14 a为方阵 b为矢量 求方阵系统的根 a fix 15 rand 3 3 b fix 15 rand 3 1 det a cond a x a b 例4 15 a和b均为方阵 求方阵系统的根 a fix 15 rand 3 3 b fix 15 rand 3 3 x a b 物理与电气工程学院 22 二 超定系统方程组Ax b A为m n矩阵 如果A列满秩 且m n 则方程没有精确解 此时方程组为超定方程组 一般采用最小二乘法 左除法 x A b建立在奇异值分解基础之上 广义逆法 x pinv A b速度较快 可靠性较差一些 实验中数据的曲线拟合就是就可以解超定方程组的方法来解 一般情况下需要将非线性问题转换为线性问题来解决 物理与电气工程学院 23 一组实验数据 时间t和测量数据y 如下表所示 认为x和y有上式的关系式 则由6组实验数据就可形成超定方程组 就可求解出C1和C2 使得将t代入函数得到的值与实际y值之间差值的最小平方和最小 如何以矩阵的形式表示 物理与电气工程学院 24 例4 16 求表中数据的最小二乘解 t 00 30 81 11 62 3 y 0 820 720 630 600 550 50 e ones size t exp t 6组实验数据变换得到的c e ycp pinv e yt1 0 0 1 2 5 y1 ones size t1 exp t1 c yp ones size t1 exp t1 cp plot t y t1 y1 ro t1 yp b 物理与电气工程学院 25 三 欠定系统程组Ax b A为m n矩阵 如果m n 即未知数的个数多于方程的个数 理论上有无穷个解 MATLAB寻求一个基本特解 左除法 x A b最少元素解广义逆法 x pinv A b最小范数解 例3 17 求解欠定系统 a fix 15 rand 2 3 b fix 15 rand 2 1 p a b 物理与电气工程学院 26 例4 18 使用两种方法求解欠定系统 a 111 11 1 b 106 p a bnorm p a p bq pinv a bnorm q a q b 27 第3节数据分析MATLAB对数据分析函数有两条约定 输入参数 向量 则不论是行向量还是列向量 运算都是对向量整体进行的输入参数 矩阵 则运算是按列进行的 一 基本统计命令 mean x 均值或期望 median x 中位数 排序后中间的那个数 或中间两个数的均值 var x 方差 std x 标准差 物理与电气工程学院 28 sum x 求和 max x 最大值 min x 最小值 sum x 求和 prod x 求积 sx cumsum x 逐元素累和 sx结构同x px cumprod x 逐元素累积 px结构同x xs xn sort x 升序排序 xs为重新排序后数据 xn为xs各元素在x中的原始位置 y sortrows x n 矩阵x以第n列为基准 按递增顺序排列 物理与电气工程学院 29 例4 19 基本统计函数的应用 A fix 10 randn 4 4 Amed median A Amean mean A Avar var A Astd std A 等同于std A 0 Astd1 std A 1 std A 2是对A方差的最佳无偏估计Amax max A Amin min A Asum sum A Aprod prod A Acumsum cumsum A Acumprod cumprod A As An sort A 物理与电气工程学院 30 有关正态分布的几个常用函数fp normpdf x mu sigma 正态分布密度函数在x处的值 mu mean X sigma std X P normcdf x0 mu sigma 正态分布随机变量小于等于x0的累积概率值 x0 norminv P mu sigma 正态分布逆累积函数 p P X x0 与normcdf是互逆的 probabilitydensity概率密度cumulativedistribution累计分布inverse相反的 物理与电气工程学院 31 附加 norm 1x 30 0 1 30 ynorm normpdf x 0 1 标准正态分布ynorm1 normpdf x 3 2 中心点由期望值确定 方差越大分布越分散plot x ynorm x ynorm1 物理与电气工程学院 32 附加例norm 3 设X 3 2 2 求P X 3 P 2 X 5 p1 normcdf 3 3 2 p2 normcdf 5 3 2 normcdf 2 3 2 物理与电气工程学院 33 附加例norm 3 设X 3 2 2 确定c 使得P X c 0 612c norminv 0 612 3 2 p normcdf c 3 2 用于验证 34 二 协方差矩阵和相关阵C cov x 如果x是向量 则返回C是该向量的方差 C是一个标量 C cov x 如果x是矩阵 则每一列相当于一个变量 返回C是该矩阵的列与列之间的协方差矩阵 C是一个列数与x相同的对称方阵 主对角线元素是x对应列矢量的方差c i j mean x i mean x i x j mean x j C cov x y x和y是同维向量或矩阵 相当于cov x y 既先将x和y按展开成列向量 再求这两个列向量的协方差 diag cov x 如果x是矩阵 则C就是该矩阵每个列向量的方差 物理与电气工程学院 35 S corrcoef x 计算x的相关系数矩阵S 将x的每一列视为一个变量 S i j 是i列与j列之间的相关系数 对于输入矩阵x 相关系数矩阵S和协方差矩阵C之间的关系为 S corrcoef x y 相当于corrcoef x y 即先将x和y按展开成列向量 再求这两个列向量的相关系 S为对称矩阵 主对角线元素值为1 物理与电气工程学院 36 例4 20 计算协方差和相关系数矩阵 x rand 10 3 y rand 10 3 z rand 1 10 cx cov x cy cov y cz cov z vz var z cxy cov x y 物理与电气工程学院 37 三 微分 差分与梯度1 微分及差分计算数值向量或矩阵的数值差分 对于一个数值向量或矩阵F diff F 就是计算F 2 m F 1 m 1 其调用格式为 Y diff F n dim 其中 F是向量或数组 n是差分阶数 dim是指定沿着数组的哪一维进行差分 difference微分 差分derivative导数 物理与电气工程学院 38 例4 21 计算三维数组的差分 a rand 3 3 2 diff a 1阶差分 列内diff a 2 第2阶差分 列内diff a 3 第3阶差分 行内diff a 4 第4阶差分 行内diff a 5 第5阶差分 页间diff a 6 第6阶差分 空 物理与电气工程学院 39 2 近似梯度gradient函数用来进行数值梯度的计算 一般情况下 函数F x y z 的梯度定义为 FX FY gradient F h1 h2 FX FY FZ gradient F h1 h2 h3 FX FY FZ分别按矩阵的列 行 页方向 h1 h2 h3指定点间隔 物理与电气工程学院 40 例4 22 计算表达式xe x 2 y 2 的梯度 v 2 0 2 2 生成 2到2间隔为0 2的自变量 x y meshgrid v 产生数据网格z x exp x 2 y 2 计算z px py gradient z 2 2 求二维梯度contour x y z 绘制等高线holdon 保持绘图quiver x y px py 绘制矢量图 抖动holdofffiguremesh x y z 41 第4节插值插值 Interpolation 是在原始数据点之间按照一定的关系插入新的数据点 以便能较准确地分析数据的变化规律 一 一维插值一维插值就是对一维函数y f x 的数据进行插值 yi interp1 x y xi method 它是在 基准 数据的基础上 研究如何 平滑 地估算出 基准 数据间其他点的函数值 原始数据点 x y x为横坐标向量 y为纵坐标向量 xi待插值点的横坐标 yi为插值函数计算出的待插值点的纵坐标 物理与电气工程学院 42 X的数据必须按单调方式排列 如果y为矩阵 则插值将按照y的列向量进行 xi为指定的拟插值点的横坐标 yi是在xi指定位置计算出的插值结果 如果xi的某元素xi i 超出x的定义范围 那么线性插值和最近邻插值方法将相应的yi i 将取值为Nan 物理与电气工程学院 43 method用于指定插值的方法 缺省时默认采用线性插值方法 物理与电气工程学院 44 例4 23 一维插值函数插值方法的对比 x 0 10 y sin x xi 0 25 10 将插值方法定义为单元数组strmod nearest linear spline cubic 将X轴标识定义为单元数组strlb a method nearest b method linear c method spline d method cubic fori 1 4yi interp1 x y xi strmod i subplot 2 2 i plot x y ro xi yi b xlabel strlb i end 物理与电气工程学院 45 spline三次样条插值函数yi spline x y xi 参数使用方法与interp1相似 并指定用三次样条方法 例4 24 三次样条插值 x0 0 10 y0 sin x0 x 0 25 10 y spline x0 y0 x plot x0 y0 or x y k 46 二 二维插值二维函数插值是基于一维函数插值的基本思想 它是对两个变量的函数z f x y 进行插值 z为矩阵 由x y确定的点上的值 zi interp2 x y z xi yi method 输入参数为原始数据点 x y z x y为两个独立的矩阵 它们也必须是按照单调方式排列的 z为函数矩阵 是由x y确定的点上的值 xi为指定插值点横坐标的数值数组 yi是指定插值点纵坐标的数值数组 method参数同一维插值 物理与电气工程学院 47 例4 25 二维插值四种方法的对比 x y z peaks 7 生成双峰函数值 x y数据间隔为1mesh x y z 绘制网格图 xi yi meshgrid 3 0 2 3 3 0 2 3 生成待插值的数据网格z1 interp2 x y z xi yi nearest z2 interp2 x y z xi yi linear z3 interp2 x y z xi yi spline z4 interp2 x y z xi yi cubic figure mesh xi yi z1 figure mesh xi yi z2 figure mesh xi yi z3 figure mesh xi yi z4 物理与电气工程学院 48 三 多维插值多维插值包括interp3和interpn函数 使用方法与interp2类似 vi interp3 x y z v xi yi zi method 例4 26 三维插值实例 x y z v flow 10 slice x y z v 69 5 2 2 2 xi yi zi meshgrid 0 25 10 3 25 3 3 25 3 vi interp3 x y z v xi yi zi figure slice xi yi zi vi 69 5 2 2 2 shadingflat 物理与电气工程学院 49 第5节快速傅利叶变换离散序列的傅立叶变换 DiscreteFourierTransform DFT 和傅立叶逆变换 InvenseDiscreteFourierTransform IDFT 在数学上定义和MATLAB软件中的定义基本相同 不同在于 MATLAB中的序列或向量元素下标从1开始 而不是数学上的从0开始 物理与电气工程学院 50 MATALB提供了fft 内置函数 iff fft2 fftn ifftn fftshift ifftshift等函数 用来计算矩阵的离散快速傅立叶变换 在数据长度是2的幂次时 采用基 2算法进行计算 计算速度会显著增加 因此在使用时 尽量使用数据长度为2的幂次或者用零来填补数据 快速傅里叶变换的结果为复数 物理与电气工程学院 51 一 函数fft和ifftY fft X X为矢量傅里叶变换序列 变换序列长度与X等长 Y fft X n 指定变换序列长度n 如果X的长度小于n 则用0来补足 如果X的长度大于n 则去掉X多出的部分 n可缺省时 视变换序列长度与X等长 Y fft X n dim 对指定维进行离散傅里叶变换 Y fft X dim n可缺省时 视变换序列长度与X等长 若X为矩阵则计算每列的傅立叶变换序列 物理与电气工程学院 52 Y fft X N 点域 N 2 1N1N 2 频域 Fs 2 Fs N0Fs 2 Fs N X的采样频率Fs 点频率 Fs N 物理与电气工程学院 53 例fft test fft实验fs 102 采样频率fs 1 ts 采样频率要大于信号最高频率的2倍ts 1 fs 采样间隔t 0 ts 5 x 100 sin 2 pi 50 t randn size t 加入高斯噪声正胘信号 频率为50HzN 512 变换序列点数 最好是2的幂y fft x N Pf fs N 点频率 fs Nf 0 N 2 1 Pf plot f abs y 1 N 2 幅频特性 物理与电气工程学院 54 55 例4 27 产生一个正弦衰减曲线 进行快速傅里叶变换 并画出幅值图 相位图 实部图和虚部图 tp 0 1 6 2048 时域采样时间间隔ts 1 6fs 6 抽样频率 根据奈奎斯特定理 抽样频率大于信号最高频率的2倍yt sin 4 pi tp exp tp 80 信号最高频率2Hzyf fft yt 快速傅里叶变换序列点数与yt长度相同pf fs length yt 变换序列点频率lee length yt 2 频谱取段的长度f 0 lee 1 pf 频域轴坐标ya abs yf 1 lee 幅值yp angle yf 1 lee 180 pi 相位 角度 plot tp yt 绘制正弦衰减曲线figure plot f ya 绘制FFT幅值曲线figure plot f yp 绘制FFT相位曲线 56 物理与电气工程学院 57 二 快速傅里叶变换的长度与运算速度fff函数的变换点数n决定变换的速度 n为2的整数次幂 运算速度最快 n为合数 fft采用质因数分解的算法 速度取决于质因数的大小 n为质数 运算速度最慢 几个函数解释isprime x 判断x是否质数factor x 对x进行因式分解cputime计算机cpu在此时的时刻 物理与电气工程学院 58 例3 28 比较快速傅里叶变换的长度与运算速度的关系x rand 70000 1 isprime 65539 判断65539是否为质数2 16 计算2的16次方factor 66000 计算66000的因数分解factor 65535 t cputime y fft x 65539 e cputime tt cputime y fft x 2 16 e cputime tt cputime y fft x 66000 e cputime tt cputime y fft x 65535 e cputime t 物理与电气工程学院 59 三 fftshift和ifftshift函数fftshift用于把傅立叶变换结果 频域数据 中的直流分量 频率为0Hz处的值 移到中间位置 Y fftshift X X是傅立叶变换结果如果X是向量 则交换X的左右半边如果X是矩阵 则交换其一 三象限和二 四象限ifftshift相当于fftshift函数的操作逆转 用法相同 物理与电气工程学院 60 例fft shift 变换结果频谱直流分量移位fs 102 ts 1 fs 采样间隔 t 0 ts 5 x 100 sin 2 pi 50 t randn size t N 1024 最好是2的幂y fft x N f 0 N 2 1 fs N 点频率 fs Nfigure 1 plot f abs y 1 N 2 figure 2 F N 2 N 2 1 fs N 频率轴坐标值plot F abs y 请判断频谱的正误figure 3 Y1 fftshift y 变换结果频谱直流分量移位F1 N 2 N 2 1 fs N plot F1 abs Y1 物理与电气工程学院 62 期中练习1 1 生成5行6列的矩阵 矩阵元素为1到61之间的偶数 要求元素数值按单下标顺序递增 2 找出矩阵中大于9的元素的行列标号 并有逻辑数组将其中元素值为6的换为8 3 将矩阵A B C组成一个新的矩阵 三个矩阵分别为新矩阵的第1 2 3行 物理与电气工程学院 63 期中练习2 4 绘制一个周期为3 的三角波 5 绘制分离饼图 其中5 6对应的扇区分离 6 用阶梯图表现函数 7 在一个图形窗口的两个子窗口内分别绘制网格图和表面图 物理与电气工程学院 64 期中练习3 8 绘制曲线 并用 标出y最大值的点 9 用多项式函数求在x 0 9时的值 10 x和y的实验测试数据如下表 用多项式似合x和y的二次关系式 并求出拟合值与原始值之间的差值平方均值 物理与电气工程学院 65 期中练习4 11 求下列方程组的解