MATLABch4数值计算

1、1第4章matlab的数值计算功能 matlab 具有出色的数值计算能力，占据世界上数值计算软件的主导地位 数值微积分 矩阵与代数方程 数据统计分析 多项式与插值24.1 数值微积分有限微分（差分） 与积分相比，数值微分非常困难。一个函数小的变化，容易产生相邻点与积分相比，数值微分非常困难。一个函数小的变化，容易产生相邻点的斜率的大的改变。由于微分这个固有的困难，所以尽可能避免数值微分，特的斜率的大的改变。由于微分这个固有的困难，所以尽可能避免数值微分，特别是对实验获得的数据进行微分。在这种情况下，最好用最小二乘曲线拟合这别是对实验获得的数据进行微分。在这种情况下，最好用最小二乘曲线拟合这种数

2、据，然后对所得到的多项式进行微分。种数据，然后对所得到的多项式进行微分。在在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 -13而而 的的 数值微分则为数值微分则为 dy=diff(y)./diff(x)。 例：要计算下列多项式在例：要计算下列多项式在 -4,

3、 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 ) % 将多项式的微分项绘图将多项式的微分项绘图 ti

4、tle(derivative of fifth-deg. equation)4其他用法：其他用法： 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 相邻元素是否相等相邻元素

5、是否相等 all(diff(x)0) 是否单调递增是否单调递增 all(diff(diff(x)=0) 是否均匀分布是否均匀分布5梯度fx = gradient（f） 求一元（函数）梯度fx,fy= gradient（f） 求二元（函数）梯度计算梯度，绘制引力线图：6 数值求和与数值积分sum（x） 按列求和cumsum （x） 按列求元素累计和trapz(x , y) 给定数据点x和y，计算y=f(x)下的梯形面积积分。cumtrapz(x , y) 给定数据点x和y，计算y=f(x)下的梯形面积累计积分。 求解定积分的数值方法多种多样，如简单的梯形法、辛普生求解定积分的数值方法多种多样，如

6、简单的梯形法、辛普生(simpson) 法、牛顿柯特斯法、牛顿柯特斯(newton-cotes)法等都是经常采法等都是经常采用的方法。它们的基本思想都是将整个积分区间用的方法。它们的基本思想都是将整个积分区间a,b分成分成n个子个子区间区间xi,xi+1，i=1,2,n，其中，其中x1=a，xn=b。这样求定积分问。这样求定积分问题就分解为求和问题。题就分解为求和问题。7matlab提供最简单的积分函数是梯形法提供最简单的积分函数是梯形法trapz(x,y)，其中其中x,y分别代表数目相同分别代表数目相同 的阵列或矩阵，而的阵列或矩阵，而y与与x的关系可以由的关系可以由是一函数型态（如是一函数

7、型态（如y=sin(x)）或是不以函数描述的离散型态）或是不以函数描述的离散型态 。我们看一简单积分式我们看一简单积分式 以下为以下为 matlab 的程式的程式 x=0:pi/100:pi; y=sin(x); k=trapz(x,y) k = 1.99988 精度可控的数值积分s1=quad（fun,a,b,tol) 采用递推自适应采用递推自适应simpson法计算积分。法计算积分。s1=quadl（fun,a,b,tol) 采用递推自适应采用递推自适应lobatto法求数值积分。法求数值积分。s2=dblqual（fun,xmin,xmax,ymin,ymax,tol) 二重二重(闭型闭

8、型)数值数值积分指令。积分指令。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.90089函数极值的数值求解函数极值的数值求解 matlab提供了基于单纯形算法求解函数极值的函数提供了基于单纯形

9、算法求解函数极值的函数fminbnd和和fminsearch， 它们分别用于求单变量函数和多变量它们分别用于求单变量函数和多变量函数的最小值，其调用格式为：函数的最小值，其调用格式为： x,fval,exitflag,output=fmin(fun,x1,x2,options) x,fval,exitflag,output=fmins(fun,x0,options)这两个函数的调用格式相似。这两个函数的调用格式相似。其中其中fmin函数用于求单变量函数的最小值点。函数用于求单变量函数的最小值点。fun是被最小化是被最小化的目标函数名，的目标函数名，x1和和x2限定自变量的取值范围。限定自变量的

10、取值范围。fmins函数用于求多变量函数的最小值点，函数用于求多变量函数的最小值点，x0是求解的初始值是求解的初始值向量。向量。 matlab没有专门提供求函数最大值的函数，但只要注意到没有专门提供求函数最大值的函数，但只要注意到-f(x)在区间在区间(a,b)上的最小值就是上的最小值就是f(x)在在(a,b)的最大值，的最大值，10常微分方程的数值求解常微分方程的数值求解基于龙格库塔法，基于龙格库塔法，matlab提供了求常微分方程数值解的函数，提供了求常微分方程数值解的函数，一般调用格式为：一般调用格式为：t,y=ode45(odefun,tspan,y0)其中其中odefun是定义是定义

11、f(t,y)的函数文件名，该函数文件必须返回一个的函数文件名，该函数文件必须返回一个列向量。列向量。tspan形式为形式为t0,tf，表示求解区间。，表示求解区间。y0是初始状态列向是初始状态列向量。量。t和和y分别给出时间向量和相应的状态向量。分别给出时间向量和相应的状态向量。11 单个高阶常微分方程处理方法单个高阶常微分方程处理方法124.2 矩阵与线性代数-矩阵运算加(+)、减(-)乘：(*） 矩阵之间的乘、向量与矩阵相乘、标量与矩阵相乘除：右除（/）、左除（） 数学上没有矩阵除法的定义。幂 转置 其他： 指数运算 y = expm(a) 还有:expm1、expm2、expm3 开方

12、y = sqrtm(a) 对数 y = logm(a)131矩阵的秩矩阵的秩矩阵线性无关的行数与矩阵线性无关的行数与列数称为矩阵的秩。在列数称为矩阵的秩。在matlab中，求矩阵秩中，求矩阵秩的函数是的函数是rank(a)。2矩阵的迹矩阵的迹矩阵的迹等于矩阵的对矩阵的迹等于矩阵的对角线元素之和，也等于角线元素之和，也等于矩阵的特征值之和。在矩阵的特征值之和。在matlab中，求矩阵的中，求矩阵的迹的函数是迹的函数是trace(a)。 a=2 1; 4 3; rank(a) ans = 2 % 表示a秩数为2且等于矩阵的列数 trace(a)ans = 5矩阵的秩与迹矩阵的秩与迹14行列式和逆行

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 6d = det(a) 求行列式的值d = 1c = inv(a) 求逆矩阵c = 3 -3 1 -3 5 -2 1 -2 115矩阵特征值矩阵特征值特征方程：ax=x 由|a-i|=0 求特征值i 由ax= i x 求特征向量x ieig求解ax=x特征值问题 d = eig(a) 只返回特征值向量d v,d = eig(a) 返回特征向量阵v和特

14、征值对角阵d v,d = eig(a,nobalance)矩阵中有与截断误差（eps)相当的元素时,计算精度更高。 v,d = eig(a,b) 计算广义特征向量阵v和广义特征值阵d，使av=bvd成立。16【例例】求下列对称矩阵】求下列对称矩阵a的特征值的特征值/特征向量特征向量a =1.0000 1.0000 0.5000a =1.0000 1.0000 0.5000 1.0000 1.0000 0.2500 1.0000 1.0000 0.2500 0.5000 0.2500 2.0000 0.5000 0.2500 2.0000 eig(a)%返回特征值ans = -0.0166 1.

15、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.7103d = -0.0166 0 0 0 1.4801 0 0 0 2.536517求解线性方程组矩阵方程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.2xzero= 1.

16、2995yzero=humps(xzero , 1.2)% evaluate at xzeroyzero= 3.5527e-15214.3 概率分布和统计数据分析数据集单变量的统计数据存储在1n或n1的向量里。多变量的数据用矩阵来表示，不同的变量放在不同的列。22数据预处理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 = fin

17、d(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),:) = ;23删除越界的数据点:load count.dat;%读数据mu=mean(count);%计算平均值mu = 32.0000 46.5417 65.5833sigma=std(count);%计算偏差sigma =25.3703 41.4057 68.0281n,p=size(count);outl

18、iers=abs(count-mu(ones(n,1),:)3*sigma(ones(n,1),:)%判断偏差过大的数据所在位置nout=sum(outliers); %计算越界数据点数 nout = 1 0 0count(any(outliers),:)=;%删除越界数据24概率分布函数概率分布函数 二项分布二项分布(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

19、,mu,sigma)r=normrnd(mu,sigma,m,n)各种概率分布的交互式观察界面各种概率分布的交互式观察界面disttool25样本分布的频数直方图描述样本分布的频数直方图描述hist指令的使用示例。指令的使用示例。randn(state,1),rand(state,31)%初始化初始化x=randn(100,1);y=rand(100,1);%生成正态和均匀分布实验样本生成正态和均匀分布实验样本%观察正态数据组的频数直方图在不同区间分段数时的变化观察正态数据组的频数直方图在不同区间分段数时的变化subplot(1,2,1),histfit(x,20)%20区间情况区间情况sub

20、plot(1,2,2),hist(y,20)%20区间情况（带正态拟合线）区间情况（带正态拟合线） -4-202402468101200.5101234567891026基本数据分析函数max(x) :求各列最大值min(x) :求各列最小值mean(x):求各列平均值median(x) :各列求中值std(x) :求各列标准差var (x) :求各列方差cov(x) :求x各列之间的协方差阵corrcoef(x):求x各列之间的相关系数271、求取最大值、求取最大值(max)y,i= max (x)：返回矩阵x的各列中的最大元素值及其该元素的位置赋予行向量y与i；当x为向量时，则y与i为单变

21、量。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 = 328 例例分别查找下面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) %查出二

22、维数组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 29【例例】试取下面两个23的二维数组x、y所有同一位置上的元素值大者构成一个新矩阵p。 x=4 5 6;1 4 8 % 产生二维数组x x = 4 5 6

23、 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 830min函数用来查取数据序列的最小值。它的用法与命令格式与max函数完全一样，所不同的是执行的结果是最小值。2、求取最小值、求取最小值(min)3、求和、求和y=sum(x)：sum(x)将返回矩阵x各列元素之和赋予行向量y；若x为二维数组，y为一个向量；若x为向量，则y为单变量。y=sum(x,dim)：按数组x的第dim维的方向的元素求其和赋予y。若dim=1，

24、为按列操作；若dim=2，为按行操作。31 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.33334、求平均值、求平均值32所谓中值，是指在数据序列中其值的大小恰好在中间。例如，数据序列9,-2,

25、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)33【例例】试分别求下面数列x1与x2的中值。x1=9 -2 5 7 12; % 奇数个元素y1=median(x)y1 = 7x2=9 -2 5 6 7 12; % 偶数个元素y2=median(

26、x)y2 = 6.500034【例例】对下面二维数组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.500035y= prod(x)：prod(x)将返回矩阵x各列元素之积赋予行向量y。若x

27、为二维数组，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 48y2= prod(x,2) y2 = 120 32367、协方差与相关系数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.

28、0000 -4.1667 2.3333 -1.5000 3.0000 -1.5000 1.0000s=corrcoef(x)相关系数 r = corrcoef(a) 1.0000 -0.8486 0.9333 -0.8486 1.0000 -0.9820 0.9333 -0.9820 1.0000374.4 多项式与卷积-多项式的表示一般多项式形式： a1xn+a2xn-1+anx+an+1在matlab中，用行向量来表示多项式。 例如： 多项式 4x4+2x2+x+1 表示成 ： 4 0 2 1 138多项式的操作(1)-多项式的相乘和相除多项式相乘（conv） w = conv(u,v)

29、返回u和v两向量的卷积多项式相除（deconv） q,r = deconv(v,u) v是被除多项式， q是商，r 是余项 v = conv(u,q)+r.39【例例】求多项式与多项式的乘积/相除a=1 8 0 0 -10a= 1 8 0 0 -10b=2 -1 3b = 2 -1 3c=conv(a,b)c = 2 15 -5 24 -20 10 -30本例的运行结果是求得一个6次多项式：2x6+15x5-5x4+24x3-20 x2+10 x-30a=1 8 0 0 -10;b=2 -1 3;p,r=deconv(a,b) p = 0.5000 4.2500 1.3750 r = 0 0

30、0 -11.3750 -14.1250商多项式p为: 0.5x2+4.25x+1.375，余多项式r为: -11.375x-14.125。40多项式的操作(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

31、 0 0 -10a = 1 8 0 0 -10 x=roots(a)41多项式的操作(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

32、 = 377 179 439 111 81 136 490 253 63942【例】对前例的4次多项式，分别取x=1.2和下面矩阵的23个元素为自变量计算该多项式的值。a=1 8 0 0 -10; % 4次多项式系数x=1.2; % 取自变量为一数值y1=polyval(a,x) y1 = -97.3043x=-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.321643多项式的操作(4)-多项式的建立(poly)可以用poly函数求矩阵的特征

33、多项式 若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与root

34、s互为逆函数；44操作(4)-多项式曲线拟合（polyfit）对已知的离散数据，采用多项式模型构造光滑曲线进行描述。p = polyfit(x,y,n) 利用已知的向量x和y所确定的数据点，采用最小二乘法构造n阶多项式 p是n阶多项式的系数向量，n是多项式阶数p,s = polyfit(x,y,n) 返回估计或预测误差的结构s45例：x=-2.0 -1.6 -1.2 -0.8 -0.4 0 0.4 0.8 1.2 1.6 2.0y=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

35、=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(原始数据,二阶多项式,十阶多项式)46最小二乘问题实际中常常需要找到一个函数来描述观测量之间的关系。可归结为确定该函数的系数，可通过求解线性方程组得到。多项式回归 用多项式拟合：y=a0+a1t+a2t2 解线性方程组: a=xy线性参数回归 用线性参数的函数拟合 如指数函数：y=a0+a1e-t+a2te-t多变量回归47基本数据拟合

36、界面拟合数据并绘制拟合曲线求拟合的残余向量对拟合的结果求插值将结果保存4849t = 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 onhold onx = ones(size(t) t t.2%二次多项式回归a = xyt1 = (0:0.1:2.5);y1 = ones(size(t1) t1 t1.2*a;plot(t1,y1,-,linewidth,2), grid onx = ones(size(t) exp(-t) t.*exp(-t);%线性回归a = xyt2 = (0:0

37、.1:2.5);y2 = ones(size(t1) exp(-t2) t2.*exp(-t2)*a;plot(t2,y2,r-,linewidth,2), grid onlegend(数据点,多项式拟合,线性拟合,2);hold off00.511.522.50.50.60.70.80.911.11.21.31.450操作(5)-多项式求导（ployder）k = polyder(p) 返回多项式p的导数k，例如： p = 1 0 -2 -5 %多项式 x3-2x-5 q = polyder(p) q = 3 0 -2 %多项式 3x2-2k = polyder(a,b) 返回多项式a、b乘

38、积的导数q,d = polyder(a,b) 返回a/b商的导数，格式为分子q，分母d51插值的概念插值在认为“测量数据”完全正确的情况下，研究如何平滑地估算出“测量数据”之间的其他点的值。当不能很快地求出所需中间点的函数时，插值是一个非常有价值的工具。可以利用已知的数据构造插值表。简单地，可以先根据测量数据绘图，然后估计所需点处的值。matlab提供了一维、二维、 三次样条等许多插值选择52插值原理内插值和查表一样53插值-一元插值interp1yi = interp1(x,y,xi) x、y是基准数据，xi是插值点的自变量 如果y是矩阵，按列对y进行插值yi = interp1(x,y,xi,method) method是插值方法 nearest：最近领域插