MATLAB---ch4数值计算.ppt

第4章MATLAB的数值计算 功能 matlab 具有出色的数值计 算能力，占据世界上数值计算 软件的主导地位 数值微积分 矩阵与代数方程 数据统计分析 多项式与插值 1 4.1 数值微积分 有限微分（差 分） 与积分相比，数值微分非常困难。一个函数小的变化，容易产生相邻点的 斜率的大的改变。由于微分这个固有的困难，所以尽可能避免数值微分，特别 是对实验获得的数据进行微分。在这种情况下，最好用最小二乘曲线拟合这种 数据，然后对所得到的多项式进行微分。 在MATLAB中，没有直接提供求数值导数的函数，只有计算向前 差分的函数diff，其调用格式为：Y = diff(X) 计算向量中连续两个元素的差Y=X(2)-X(1) X(3)-X(2) . X(n)-X(n-1) 例如：A = 9 -2 3 0 1 5 4; diff(A) ans = -11 5 -3 1 4 -1 2 而 的 数值微分则为 dy=diff(y)./diff(x)。 例：要计算下列多项式在 -4, 5 区间的微分 x=linspace(-4,5); % 产生100个x的离散点 p=1 -3 -11 27 10 -24; f=polyval(p,x); plot(x,f) % 将多项式函数绘图 title(Fifth-deg. equation) dfb=diff(f)./diff(x); % 注意要分别计算diff(f)和diff(x) xd=x(2:length(x); % 注意只有99个df值，对应x2,x3,.,x100的点 plot(xd,dfb ) % 将多项式的微分项绘图 title(Derivative of fifth-deg. equation) 3 其他用法： DX=diff(X)：计算向量X的向前差分， DX(i)=X(i+1)-X(i)，i=1,2,n-1。 DX=diff(X,n)：计算X的n阶向前差分。例 如，diff(X,2)=diff(diff(X)。 DX=diff(A,n,dim)：计算矩阵A的n阶差分 ，dim=1时(缺省状态)，按列计算差分；dim=2 ，按行计算差分。 几个应用： diff(x)=0 相邻元素是否相等 all(diff(x)0) 是否单调递增 all(diff(diff(x)=0) 是否均匀分布 4 梯度 FX = gradient（F） 求一元（函数）梯度 FX,FY= gradient（F） 求二元（函数）梯度 计算梯度，绘制引力线图： 5 数值求和与数值积分 sum（X） 按列求和 cumsum （X） 按列求元素累计和 trapz(x , y) 给定数据点x和y，计算y=f(x)下 的梯形面积积分。 cumtrapz(x , y) 给定数据点x和y，计算y=f(x) 下的梯形面积累计积分。 求解定积分的数值方法多种多样，如简单的梯形法、辛普生 (Simpson)法、牛顿柯特斯(Newton-Cotes)法等都是经常采用 的方法。它们的基本思想都是将整个积分区间a,b分成n个子区 间xi,xi+1，i=1,2,n，其中x1=a，xn=b。这样求定积分问题就 分解为求和问题。 6 MATLAB提供最简单的积分函数是梯形法trapz(x,y)， 其中x,y分别代表数目相同 的阵列或矩阵，而y与x的关系可以由是 一函数型态（如y=sin(x)）或是不以函数描述的离散型态 。 我们看一简单积分式 以下为 MATLAB 的程式 x=0:pi/100:pi; y=sin(x); k=trapz(x,y) k = 1.9998 7 精度可控的数值积分 S1=quad（fun,a,b,tol) 采用递推自适应Simpson法计算积分。 S1=quadl（fun,a,b,tol) 采用递推自适应Lobatto法求数值积分 。 S2=dblqual（fun,xmin,xmax,ymin,ymax,tol) 二重(闭型)数值 积分指令。 S3=triplequal （fun,xmin,xmax,ymin,ymax,zmin,zmax,tol)三重( 闭型)数值积分指令。 例: 求定积分。 (1) 建立被积函数文件fesin.m。 function f=fesin(x) f=exp(-0.5*x).*sin(x+pi/6); (2) 调用数值积分函数quad求定积分。 S=quad(fesin,0,3*pi) S = 0.90088 函数极值的数值求解 MATLAB提供了基于单纯形算法求解函数极值的函数 fminbnd和fminsearch，它们分别用于求单变量函数和多变量函 数的最小值，其调用格式为： x,fval,exitflag,output=fmin(fun,x1,x2,options) x,fval,exitflag,output=fmins(fun,x0,options) 这两个函数的调用格式相似。 其中fmin函数用于求单变量函数的最小值点。fun是被最小化的 目标函数名，x1和x2限定自变量的取值范围。 fmins函数用于求多变量函数的最小值点，x0是求解的初始值向 量。 MATLAB没有专门提供求函数最大值的函数，但只要注意到- f(x)在区间(a,b)上的最小值就是f(x)在(a,b)的最大值， 9 常微分方程的数值求解 基于龙格库塔法，MATLAB提供了求常微分方程数值解的函数 ，一般调用格式为： t,y=ode45(odefun,tspan,y0) 其中odefun是定义f(t,y)的函数文件名，该函数文件必须返回一个 列向量。tspan形式为t0,tf，表示求解区间。y0是初始状态列向量 。t和y分别给出时间向量和相应的状态向量。 10 单个高阶常微分方程处理方法 11 4.2 矩阵与线性代数-矩阵运 算 加(+)、减(-) 乘：(*） 矩阵之间的乘、向量与矩阵相乘、标量 与矩阵相乘 除：右除（/）、左除（） 数学上没有矩阵除法的定义。 幂 转置 其他： 指数运算 Y = expm(A) 还有:expm1、 expm2、expm3 开方 Y = sqrtm(A) 对数 Y = logm(A) 12 1矩阵的秩 矩阵线性无关的行数与 列数称为矩阵的秩。在 MATLAB中，求矩阵秩 的函数是rank(A)。 2矩阵的迹 矩阵的迹等于矩阵的对 角线元素之和，也等于 矩阵的特征值之和。在 MATLAB中，求矩阵的 迹的函数是trace(A)。 A=2 1; 4 3; rank(A) ans = 2 % 表示A秩数为2且 等于矩阵的列数 trace(A) ans = 5 矩阵的秩与迹 13 行列式和逆 行列式：d = det(A) 计算方阵的行 列式值 矩阵的逆：C= inv(A) AC=CA=I ，则 ：C=A-1 其中A是非奇 异方阵 实际数值计算 中较少使用 A = pascal(3) 产生33的 对称矩阵 A = 1 1 1 1 2 3 1 3 6 d = det(A) 求行列式的值 d = 1 C = inv(A) 求逆矩阵 C = 3 -3 1 -3 5 -2 1 -2 1 14 矩阵特征值 特征方程：Ax=x 由|A-I|=0 求特征值i 由Ax= i x 求特征向量x i eig求解Ax=x特征值问题 d = eig(A) 只返回特征值向量d V,D = eig(A) 返回特征向量阵V和特 征值对角阵D V,D = eig(A,nobalance)矩阵中有与 截断误差（eps)相当的元素时,计算精度更高 。 V,D = eig(A,B) 计算广义特征向量阵 V和广义特征值阵D，使AV=BVD成立。 15 【例】求下列对称矩阵A的特征值/特征 向量 A =1.0000 1.0000 0.5000 1.0000 1.0000 0.2500 0.5000 0.2500 2.0000 eig(A)%返回特征值 ans = -0.0166 1.4801 2.5365 V,D = eig(A) %返回特征值及特征向量 V = 0.7212 0.4443 0.5315 -0.6863 0.5621 0.4615 -0.0937 -0.6976 0.7103 D = -0.0166 0 0 0 1.4801 0 0 0 2.5365 16 求解线性方程组 矩阵方程AX=b的解 X=Ab 矩阵方程XA=b的解 X=A/b 三种方程组情况的求解（对mn矩阵A ） m=n, 恰定方程组，可以精确求解 mn, 超定方程组，可以求得最小二乘 解 mxzero=fzero( humps , 1.2)% look for a zero near 1.2 xzero= 1.2995 yzero=humps(xzero , 1.2)% evaluate at xzero yzero= 3.5527e-15 20 4.3 概率分布和统计数据分析 数据集 单变量的统计数据存储在1n或n1 的向量里。 多变量的数据用矩阵来表示，不同 的变量放在不同的列。 21 数据预处理 MATLAB用NaN表示非数 用NaN处理缺少的数据为数据分析带来方便 TF = isnan(A) 返回与A同维的矩阵TF，A中为 NaN的得到逻辑1，否则为逻辑0 a = -2 -1 0 1 2 isnan(0./a) Warning: Divide by zero. ans = 0 0 1 0 0 删除NaN i = find(isnan(x);x = x(i) x = x(find(isnan(x) x = x(isnan(x); x(isnan(x) = ; X(any(isnan(X),:) = ; 删除任意含有NaN的 行 function X = excise(X) X(any(isnan(X),:) = ; 22 删除越界的数据点: load count.dat;%读数据 mu=mean(count);%计算平均值mu = 32.0000 46.5417 65.5833 sigma=std(count);%计算偏差sigma =25.3703 41.4057 68.0281 n,p=size(count); outliers=abs(count-mu(ones(n,1),:)3*sigma(ones(n,1),:) %判断偏差过大的数据所在位置 nout=sum(outliers); %计算越界数据点数 nout = 1 0 0 count(any(outliers),:)=;%删除越界数据 23 概率分布函数 二项分布(Binomial distribution) Pk=binopdf(k,N,p) Fk=binocdf(k,N,p) R=binordf(N,p,m,n) 正态分布（Normal distribution） Px=normpdf(x,Mu,Sigma) Fx=normcdf(x,Mu,Sigma) R=normrnd(Mu,Sigma,m,n) 各种概率分布的交互式观察界面 disttool 24 样本分布的频数直方图描述 hist指令的使用示例。 randn(state,1),rand(state,31)%初始化 x=randn(100,1);y=rand(100,1);%生成正态和均匀分布实验样本 %观察正态数据组的频数直方图在不同区间分段数时的变化 subplot(1,2,1),histfit(x,20)%20区间情况 subplot(1,2,2),hist(y,20)%20区间情况（带正态拟合线） 25 基本数据分析函数 max(X) :求各列最大值 min(X) :求各列最小值 mean(X):求各列平均值 median(X) :各列求中值 std(X) :求各列标准差 var (X) :求各列方差 Cov(X) :求X各列之间的协方差阵 Corrcoef(x):求X各列之间的相关系数 26 1、求取最大值(max) Y,I= max (X)：返回矩阵X的各列中的最大元素 值及其该元素的位置赋予行向量Y与I；当X为向量时，则Y与I 为单变量。 Y,I=max(X,DIM)：按数组X的第DIM维的方 向查取其最大的元素值及其该元素的位置赋予向量Y与I。 【例】查找下面数列x的最大值。 x=3 5 9 6 1 8 % 产生数列x x = 3 5 9 6 1 8 y=max(x) % 查出数列x中的最大值赋予 y y = 9 y,l=max(x) % 查出数列x中的最大值及其该 元素的位置赋予y,l y = 9 l = 3 27 例分别查找下面34的二维数组x中各列和各行元素中的最大值 。 x=1 8 4 2;9 6 2 5;3 6 7 1 % 产生二维数组x x = 1 8 4 2 9 6 2 5 3 6 7 1 y=max(x) %查出二维数组x中各列元素的最大值产生赋予 行向量y y = 9 8 7 5 y,l=max(x) %查出二维数组x中各列元素的最大值及其这些元素的行 下标赋予y,l y = 9 8 7 5 l = 2 1 3 2 y,l=max(x, ,1) % 本命令的执行结果与上面命令完全相 同 y = 9 8 7 5 l = 2 1 3 2 y,l=max(x, ,2) % 由于本命令中DIM=2,故查找操作在各 行中进行 y = 8 l = 2 9 1 7 3 28 【例】试取下面两个23的二维数组x、 y所有同一位置上的元素值大者构成一个 新矩阵p。 x=4 5 6;1 4 8 % 产生二维数组x x = 4 5 6 1 4 8 y=1 7 5;4 5 7 % 产生二维数组y y = 1 7 5 4 5 7 p=max(x,y)% 在x,y同一位置上的两个元素中 查找出最大值 % 赋予与x,y同样大小的二 维数组p p = 4 7 6 4 5 8 29 min函数用来查取数据序列的最小值。 它的用法与命令格式与MAX函数完全 一样，所不同的是执行的结果是最小值。 2、求取最小值(min) 3、求和 Y=sum(X)：sum(X)将返回矩阵X各列元素之和 赋予行向量Y；若X为二维数组，Y为一个向量；若X 为向量，则Y为单变量。 Y=sum(X,DIM)：按数组X的第DIM维的方向 的元素求其和赋予Y。若DIM=1，为按列操作；若 DIM=2，为按行操作。 30 Y= mean(X)： mean (X)将返回矩阵X各列元素的 平均值赋予行向量Y。若X为二维数组，Y为一个向量 ；若X为向量，则Y为单变量。 Y= mean(X,DIM)：按数组X的第DIM维的方向的 元素求其平均值赋予向量Y。若DIM=1，为按列操作； 若DIM=2，为按行操作。例如： x=4 5 6;1 4 8; y1= mean(x,1) y1 = 2.5000 4.5000 7.0000 y2= mean(x,2) y2 = 5.0000 4.3333 4、求平均值 31 所谓中值，是指在数据序列中其值的大 小恰好在中间。例如，数据序列9,-2,5, 7,12的中值为7 。 如果数据为偶数个时，则中值等于中间 的两项之平均值。 median函数调用的命令格式有： Y=median(X)：返回矩阵X各列元素的中 值赋予行向量Y。若X为二维数组，Y为一个向量；若 X为向量，则Y为单变量。 Y=median(X,DIM)：按数组X的第DIM 维方向的元素求其中值赋予向量Y。若DIM=1，为按列 操作；若DIM=2，为按行操作。 5、求中值(median) 32 【例】试分别求下面数列x1与x2的中 值。 x1=9 -2 5 7 12; % 奇数个元素 y1=median(x) y1 = 7 x2=9 -2 5 6 7 12; % 偶数个元素 y2=median(x) y2 = 6.5000 33 【例】对下面二维数组x，试从不同维方 向求出其中值。 x=1 8 4 2;9 6 2 5;3 6 7 1 % 产生一个二维数组x x = 1 8 4 2 9 6 2 5 3 6 7 1 y0=median(x) % 按列操作 y0 = 3 6 4 2 y1=median(x,1) % 此时DIM=1,故按列操作,结果y1为行向 量 y1 = 3 6 4 2 y2=median(x,2)% 此时DIM=2,故按行操作, 结果y2为列向 量 y2 = 3.0000 5.5000 4.5000 34 Y= prod(X)：prod(X)将返回矩阵X各列 元素之积赋予行向量Y。若X为二维数组，Y 为一个向量；若X为向量，则Y为单变量。 Y= prod(X,DIM)：按数组X的第DIM维 的方向的元素求其积赋予向量Y。若DIM=1， 为按列操作；若DIM=2，为按行操作。 6、求积 x=4 5 6;1 4 8; y1= prod(x,1) y1 = 4 20 48 y2= prod(x,2) y2 = 120 32 35 7、协方差与相关系数 C=cov(X)矩阵X的协方差 定义：cov(x1,x2)=E(x1-1)(x2- 2) 其中： =Exi，是x期望值 协方差矩阵： A = -1 1 2 ; -2 3 1 ; 4 0 3. C =cov(A) 10.3333 -4.1667 3.0000 -4.1667 2.3333 - 1.5000 3.0000 -1.5000 1.0000 s=corrcoef(X)相关系数 R = corrcoef(A) 1.0000 -0.8486 0.9333 -0.8486 1.0000 - 0.9820 0.9333 -0.9820 1.0000 36 4.4 多项式与卷积 -多项式的表示 一般多项式形式： a1xn+a2xn-1+anx+an+1 在MATLAB中，用行向量来表示多项 式。 例如： 多项式 4x4+2x2+x+1 表示成 ： 4 0 2 1 1 37 多项式的操作(1)-多项式的相乘和 相除 多项式相乘（conv） w = conv(u,v) 返回u和v两向量的 卷积 多项式相除（deconv） q,r = deconv(v,u) V是被除多项式， q是商，r 是余项 v = conv(u,q)+r. 38 【例】求多项式与多项式的乘积/ 相除 A=1 8 0 0 -10 A= 1 8 0 0 -10 B=2 -1 3 B = 2 -1 3 C=conv(A,B) C = 2 15 -5 24 -20 10 -30 本例的运行结果是求得一个6次 多项式： 2x6+15x5-5x4+24x3-20x2+10x-30 A=1 8 0 0 -10; B=2 -1 3; P,r=deconv(A,B) P = 0.5000 4.2500 1.3750 r = 0 0 0 -11.3750 -14.1250 商多项式P为: 0.5x2+4.25x+1.375， 余多项式r为: -11.375x-14.125。 39 多项式的操作(2)求根 (roots) 命令格式：x=roots(A) 这里A为多项式的系数A(1),A(2),A(N),A(N+1) 的行向量； 以列向量的形式返回多项式的根X。 【例】试用ROOTS函数求多项式x4+8x3-10的 根 这是一个4次多项式，它的五个系数依次 为：1,8,0,0,-10。下面先产生多项式系数的向量 A，然后求根： x = -8.0194 -0.5075 + 0.9736i -0.5075 - 0.9736i 1.0344 A=1 8 0 0 -10 A = 1 8 0 0 -10 x=roots(A) 40 多项式的操作(3)-多项式求值（ polyval） polyval函数用来求代数多项式的值，调用的命 令格式为： y = polyval(p,x) 返回向量p所代表的多项式在x处的值。返回的值 赋值给y。 若x为一数值，则Y也为一数值。 若x为向量或矩阵，则对向量或矩阵中的每个元 素求其多项式的值。 Y = polyvalm(p,X) 在矩阵意义下，返回向量p在X处的值 P=1 0 -2 -5; %多项式 x3-2x-5 X = 2 4 5; -1 0 3; 7 1 5; Y = polyvalm(p,X) %求：X3-2X-5I Y = 377 179 439 111 81 136 490 253 639 41 【例】对前例的4次多项式，分别取 x=1.2和下面矩阵的23个元素为自变量计 算该多项式的值。 A=1 8 0 0 -10; % 4次多项式系数 x=1.2; % 取自变量为一 数值 y1=polyval(A,x) y1 = -97.3043 x=-1 1.2 -1.4; 2 -1.8 1.6; % 给 出一个矩阵x y2=polyval(A,x) y2 = -17.0000 5.8976 -28.1104 70.0000 -46.1584 29.3216 42 多项式的操作(4)-多项式的建立 (poly) 可以用poly函数求矩阵的特征多项式 若x为NN的矩阵x，则poly(x)返回一个向量 赋值给A，该向量的元素为矩阵x的特征多项式之系数 ：A(1),A(2),A(N),A(N+1)。例如： A = 1.2 3 -0.9; 5 1.75 6; 9 0 1; p=poly(A) p = 1.0000 -3.9500 -1.8500 -163.2750 roots(p) %求得矩阵特征值，可以用eig直接求 若已知多项式的全部根，则可以用poly函数 建立起该多项式。 调用它的命令格式是：A=poly(x) 若x为具有N个元素的向量，则poly(x)建立以x 为其根的多项式，且将该多项式的系数赋值给向量A 。在此种情况下，poly与roots互为逆函数；43 操作(4)-多项式曲线拟合（ polyfit） 对已知的离散数据，采用多项式模型 构造光滑曲线进行描述。 p = polyfit(x,y,n) 利用已知的向量x和y所确定的数据点 ，采用最小二乘法构造n阶多项式 p是n阶多项式的系数向量，n是多项式 阶数 p,S = polyfit(x,y,n) 返回估计或预测误差的结构s 44 例： x=-2.0 -1.6 -1.2 -0.8 -0.4 0 0.4 0.8 1.2 1.6 2.0 y=2.8 2.96 2.54 3.44 3.56 5.4 6.0 8.4 9.5 13.3 15. %原始数据 p1=polyfit(x,y,2) %二阶拟合多项式 p2=polyfit(x,y,10) %十阶拟合多项式 x1=linspace(-2,2,100) y1=polyval(p1,x1) %二阶拟合数据 y2=polyval(p2,x1) %十阶拟合数据 plot(x,y,o,x1,y1,r:,x1,y2,k) legend(原始数据,二阶多项式,十阶多项式) 45 最小二乘问题 实际中常常需要找到一个 函数来描述观测量之间的关 系。 可归结为确定该函数的系 数，可通过求解线性方程组 得到。 多项式回归 用多项式拟合： Y=a0+a1t+a2t2 解线性方程组: a=Xy 线性参数回归 用线性参数的函数拟 合 如指数函数： Y=a0+a1e-t+a2te-t 多变量回归 46 基本数据拟合界面 拟合数据并绘制拟合曲 线 求拟合的残余向量 对拟合的结果求插值 将结果保存 47 48 t = 0 .3 .8 1.1 1.6 2.3;%原始数据 y = 0.5 0.82 1.14 1.25 1.35 1.40; plot(t,y,o), grid on hold on X = ones(size(t) t t.2 %二次多项式回归 a = Xy T1 = (0:0.1:2.5); Y1 = ones(size(T1) T1 T1.2*a; plot(T1,Y1,-,LineWidth,2), grid on X = ones(size(t) exp(-t) t.*exp(-t); %线性回归 a = Xy T2 = (0:0.1:2.5); Y2 = ones(size(T1) exp(-T2) T2.*exp(-T2)*a; plot(T2,Y2,r-,LineWidth,2), grid on legend(数据点,多项式拟合,线性拟合,2);hold off 49 操作(5)-多项式求导（ ployder） k = polyder(p) 返回多项式p的导数k，例如： p = 1 0 -2 -5 %多项式 x3-2x-5 q = polyder(p) q = 3 0 -2 %多项式 3x2-2 k = polyder(a,b) 返回多项式a、b乘积的导数 q,d = polyder(a,b) 返回a/b商的导数，格式为分子q， 分母d 50 插值的概念 插值在认为“测量数据”完全正确的情况下 ，研究如何平滑地估算出“测量数据”之间的其 他点的值。 当不能很快地求出所需中间点的函数时，插 值是一个非常有价值的工具。 可以利用已知的数据构造插值表。 简单地，可以先根据测量数据绘图，然后估 计所需点处的值。 Matlab提供了一维、二维、 三次样条等许多 插值选择 51 插值原理 内插值和查表一样 52 插值-一元插值interp1 yi = interp1(x,y,xi) x、y是基准数据，xi是插值点的自变量 如果y是矩阵，按列对y进行插值 yi = interp1(x,y,xi,method) Method是插值方法 nearest：最近领域插值 linear：线性插值(缺省方法) spline 三次样条插值 pchip：分段三次hermite插值 cubic：同pchip v5cubic：MATLAB用的三次插值 yi = interp1(x,Y,xi,method,extrap