MATLAB 7.0 采用行向量来表示多项式，将多项式的系数按降幂次序存放在行向量中。
实例：设计一个函数 poly2str()，实现把一个行向量表示的多项式转换为常见的字符串表 示的多项式，
%poly2str.m %把多项式的行向量表示转换为字符串表示
function Y=poly2str(X) %X 是表示一个多项式的向量 %Y 多项式的字符串表示 %输入检查，如果 X 不是一个向量则退出
if isvector(X)==0,
     disp('输入错误：输入 X 不是一个向量，请输入一个代表多项式的向量!');
     return;                          %函数返回
end; Y='';                                %输出字符串
n=length(X);
for i=1:n,                            %把多项式的每一次幂转换为字符串
     if(i~=1&&X(i)>0)                 %如果是正系数，必须添加‘+’字符             
        Y=[Y '+'];
    end;                                         %输出系数
    if(X(i)==0),                       %如果该次幂系数为 0，则不输出字符串
        continue;
    elseif(X(i)==1&&i~=n),             %如果该次幂系数为 1，可以不输出系数，只输出 x^n
        Y=Y;
    else
        Y=[Y num2str(X(i))];           %其他情况输出系数
    end;                                      %输出 x^n
    if(i==n-1),                         %1次幂输出字符串’x’
        Y=[Y 'x'];
    elseif(i==n),                       %0次幂不输出 x^n
        Y=Y;
    else
        Y=[Y 'x^' num2str(n-i)];         %其他情况输出 x^n
    end;                                          %如果不是最后一项，输出’+’
end; if(Y(1)=='+')             %修正如果 0 次幂为 0 时，造成字符串末尾有多余的字符串’ + ’
    Y(1)=[];
end;  
实例：，设计一个函数实现从多项式的字符串表示转换为多项式的行向量表示，要求输入 不必降幂排列，并且函数具有实现多项式同类项合并功能，代码设置如下：
%str2poly.m %把多项式的字符串表示转换为行向量表示
function Y=str2poly(X) %X 是字符串形式的多项式 %Y 是行向量形式的多项式 %输入格式检查
if(ischar(X)==0)
      disp('输入错误：输入 X 必须是一个字符串！');
end; %用正则表达式寻找’+’或者’-’的下标位置
index=regexp(X,'\+|\-'); %多项式的项数
L=length(index); %用于存储多项式每一项信息的单元字符串矩阵
term=cell(1,L+1);
term(1)=cellstr(X(1:(index(1)-1)));
for i=1:L-1,
    term(i+1)=cellstr(X(index(i):(index(i+1)-1)));
end;
term(L+1)=cellstr(X(index(L):end)); %如果第一项为空，则删除该项
if(isempty(char(term(1))))
    term(1)=[]; %多项式的项数减一   
    L=L-1;
end;   %多项式系数矩阵
coefficient=[]; %多项式幂次矩阵，它与多项式系数矩阵一一对应
power=[];
for i=1:L+1, %单项多项式字符串表达式
    substring=char(term(i)); %用正则表达式，寻找字符串’x^’，由于’^’是正则表达式中的特殊字符，多以用’\^’来表示
    index2=regexp(substring,'x\^');
    if(isempty(index2)==0), %如果匹配上        
        if(index2(1)==1), %单项多项式字符串为’x^*’形式
            coefficient=[coefficient 1];
            power=[power str2num(substring((index2(1)+2):end))];
        elseif(index2(1)==2),    %单项多项式字符串为’+x^*’或者’-x^*’，’3x^*’形式
           if(substring(1)=='+'),
               coefficient=[coefficient 1];
               power=[power str2num(substring((index2(1)+2):end))];
           elseif(substring(1)=='-'),
               coefficient=[coefficient -1];
               power=[power str2num(substring((index2(1)+2):end))];
           end;
        else %单项多项式字符串为’2.2x^*’
            coefficient=[coefficient str2num(substring(1:(index2(1)-1)))];
            power=[power str2num(substring((index2+2):end))]; 
        end;
    else %单项多项式字符串不含’x^’ %用正则表达式，寻找字符串’x’
        index2=regexp(substring,'x');
        if(isempty(index2)==0), %如果匹配上’x’ 
            if(index2(1)==1), %单项多项式字符串为’x*’形式
                coefficient=[coefficient 1];
                power=[power 1];
            elseif(index2(1)==2),  %单项多项式字符串为’+x*’或者’-x*’，’3x*’形式
               if((substring(1)=='+')==1),
                   coefficient=[coefficient 1];
                   power=[power 1];
               elseif(substring(1)=='-'),
                   coefficient=[coefficient -1];
                   power=[power 1];
               else %单项多项式字符串为’2.2x*’
                   coefficient=[coefficient str2num(substring(1:(index2-1)))];
                    power=[power 1];
               end; 
            else
                coefficient=[coefficient str2num(substring(1:(index2-1)))];
                power=[power 1]; 
            end;
        else %单项多项式字符串不含’x^’,’x’,则断定它是常数项
           coefficient=[coefficient str2num(substring)];
           power=[power 0]; 
        end; 
    end;
end; %合并同类项 
N=max(power)+1;
Y=zeros(1,N);
for i=1:N,
    index3=find(power==(N-i));
    Y(i)=sum(coefficient(index3));
end;
调用代码
X=’-10x^2+8x^2-x^3+4x^5+x^3-10+9-x^7’;
P=str2poly(X)
Str=poly2str(P)

多项式求值
在 MATLAB 7.0 中使用函数 polyval()来计算多项式的值
实例：求多项式的值和矩阵多项式的值，请注意这两者的区别，具体代码如下：
%polyval_example.m %多项式求值
str1='x^2+2x+3';
p1=str2poly(str1);
A=[1 0;
   0 -1]; p1_A=polyval(p1,A)                %求多项式值
 p1_Am=polyvalm(p1,A)             %求矩阵多项式的值

多项式乘法和多项式除法
多项式乘法和除法运算也就是多项式向量的卷积和解卷积运算。在 MATLAB 7.0 中，函 数 conv()和 deconv()用来实现这些运算。这两个函数的调用格式如下。 • w = conv(u,v)，实现向量 u,v 的卷积，在代数上相当于多项式 u 乘上多项式； • [q,r] = deconv(v,u)，实现解卷积运算，他们之间的关系为 v = conv(u,q)+r。在代数上相 当于实现多项式 u 除以 v，得到商 q 和余多项式 r。
实例：求多项式 x2+2x+3 与多项式 2x2+3x+4 的乘积
%poly_conv.m %多项式乘法
str1=’x^2+2x+3’;
str2=’2x^2+3x+4’;
p1=str2poly(str1);
p2=str2poly(str2);
p3=conv(p1,p2);
str3=poly2str(p3)

例如，求上个例子中的多项式乘积结果 2x4+7x3+16x2+17x+12 除以多项式 x2+2x+3 的商，
%poly_deconv.m %多项式除法
str3=’2x^4+7x^3+16x^2+17x+12’;
str1=’x^2+2x+3’;
p3=str2poly(str3);
p1=str2poly(str1);
[p2 r]=deconv(p3,p1);
str2=poly2str(p2)
r

多项式的导数和微分
1．多项式的导数
在 MATLAB 7.0 中使用函数 polyder()来计算多项式的导数，其调用方式如下： • k = polyder(p)，返回多项式 p 的导数； • k = polyder(a,b)，返回多项式 a 与多项式 b 乘积的导数； • [q,d] = polyder(b,a)，返回多项式 a 除以 b 的商的导数，并以 q/d 格式表示
实例：，对多项式求导，代码设置如下：
%polyder_example.m   %多项式求导数的示例
str1=’x^2+2x+3’;
str2=’2x^2+3x+4’;
p1=str2poly(str1);
p2=str2poly(str2); %求多项式’x^2+2x+3’的导数
q1=polyder(p1);
str1_der=poly2str(q1) %求多项式’2x^2+3x+4’的导数
q2=polyder(p2);
str2_der=poly2str(q2) %求多项式’x^2+2x+3’与’2x^2+3x+4’乘积的导数
c=polyder(p1,p2);
str3_der=poly2str(c) %求多项式’x^2+2x+3’除以’2x^2+3x+4’的商的导数
[q d]=polyder(p1,p2);
str_der_q=poly2str(q)
str_der_q=poly2str(d)

．多项式的积分
在 MATLAB 7.0 中使用函数 polyint()来计算多项式的积分，其调用方式如下： • polyint(p,k)，返回多项式 p 的积分，设积分的常数项为 k； • polyint(p)，返回多项式 p 的积分，设积分的常数项为 0
实例：求多项式 3x2+4x+5 的积分，代码设置如下：
%polyint_example.m   %多项式求积分的示例
str1=’3x^2+4x+5’;
p1=str2poly(str1); %求积分并设常数项为 5
p2=polyint(p1,5);
str2=poly2str(p2) %求积分并设常数项为 0
p3=polyint(p1);
str3=poly2str(p3)

多项式的根和由根创建多项式
1．多项式的根
在 MATLAB 7.0 中使用函数 roots()来求多项式的根，其调用方式如下： • r = roots(c)，返回多项式 c 的所有根 r，r 是向量，其长度等于根的个数。
实例：求多项式 2x3-2x2-8x+8 的根，代码设置如下
%roots_example.m %多项式求根
str1=’2x^3-2x^2-8x+8’;
p1=str2poly(str1); r=roots(p1)                      %求多项式根

2.由根创建多项式
与多项式求根相反的过程是由根创建多项式，它由函数 poly()来实现，其调用格式如下： • p = poly(r)，输入 r 是多项式所有根，返回值为代表多项式的行向量形式； • p = poly(A)，输入是 N×N 的方阵，返回值 p 是长度为 N+1 的行向量多项式，它是矩阵 A 的特征多项式，也就是说多项式 p 的根是矩阵 A 的特征值。
实例：，由根[-2 2 1]创建多项式
%poly_example.m %由根创建多项式
r=[-2 2 1];
p=poly(r);
str=poly2str(p)
实例：，求矩阵的特征多项式，代码设置如下
%poly_example2.m %求矩阵 A 的特征多项式
A =[1    2    3;
    4    5    6;
    7    8    0];
p=poly(A);
str=poly2str(p)；
r=roots(p)            %求多项式 p 的根
A_eig=eig(A)          %求矩阵 A 的特征值

多项式部分分式展开
函数 residue()可以将多项式之比用部分分式展开，也可以将一个部分分式表示为多项式 之比。
实例：，求多项式之比的部分分式展开，然后在把部分分式表示为多项式之比的形式，代 码设置如下：
%residue_example.m %求多项式之比的部分分式，然后把部分分式表示为多项式之比的形式
str1=’2x^3+3x^2-4x+1’;
str2=’x^2-3x+2’;
p1=str2poly(str1);
p2=str2poly(str2); [r,p,k]=residue(p1,p2)                    %求多项式之比的部分分式
[p1_res,p2_res]=residue(r,p,k);           %把部分分式表示为多项式之比的形式 str1_res=poly2str(p1_res)               %把多项式显示为字符串形式
str2_res=poly2str(p2_res)

多项式曲线拟合
函数 polyfit()采用最小二乘法对给定数据进行多项式拟合，最后给出多项式的系数。该 函数的调用方式如下： • p = polyfit(x,y,n)，采用 n 次多项式 p 来拟合数据 x 和 y，从而使得 p(x)与 y 最小均方差 最小
实例：下面例子说明函数 polyfit()的用法，并讨论采用不同多项式阶数对拟合结果的影 响
%polyfit_example.m %说明函数 polyfit()的用法，并讨论采用不同多项式阶数对拟合结果的影响
x=0:0.2:10;
y=0.25*x+20*sin(x); %5 阶多项式拟合
p5=polyfit(x,y,5);
y5=polyval(p5,x); %8 阶多项式拟合
p8=polyfit(x,y,8);
y8=polyval(p8,x); %60 阶多项式拟合
p60=polyfit(x,y,60);
y60=polyval(p60,x); %画图
hold on;
plot(x,y,’ro’);
plot(x,y5,’b--’);
plot(x,y8,’b:’);
plot(x,y60,’r-.’);
xlabel(’x’);
ylabel(’y’); legend(’原始数据’,’5 阶多项式拟合’,’8 阶多项式拟合’,’60 阶多项式拟合’);
%得到输出如图 6-1 所示。由图 6-1 可以看出：使用 5 次多项式拟合时，拟合得到的结果 比较差。而使用 8 次多项式拟合时，得到的结果与原始数据符合得很好。但使用 60 次多项式 拟合时，拟合的结果非常差。可见用多项式拟合必须选择适中的阶数，而不是阶数越高精度 越高。 

曲线拟合图形用户接口
为了方便用户的使用，MATLAB 7.0 提供了支持曲线拟合的图形用户接口。它位于 “Figure”窗口的“Tools\Basic Fitting”菜单中。为了使用该工具，首先用待拟合的数据画图， 代码设置如下：
x=0:0.2:10;
y=0.25*x+20*sin(x);
plot(x,y,’ro’);

 插值
插值是在已知数据之间寻找估计值的过程。在信号处理和图像处理中，插值是极其常用 的方法。MATLAB 7.0 提供大量的插值函数。这些函数在获得数据的平滑度、时间复杂度和 空间复杂度方面上有不同的性能。 插值函数保存在 MATLAB 7.0 工具箱的 polyfun 子目录下
函数名 功能描述 pchip 分段三次厄米多项式插值 interp1 一维插值 interp1q 一维快速插值 interpft 一维快速傅立叶插值方法 interp2 二维插值 interp3 三维插值 interpn N 维插值 griddata 栅格数据插值 griddata3 3 维栅格数据插值 griddatan N 维栅格数据插值 spline 三次样条插值 ppval 分段多项式求值

一维插值
一位插值就是对一维函数 y=f(x)进行插值，图 6-6 显示了插值的含义，实心点(x,y)代表的 是已知数据，空心点(xi,yi)的横坐标 xi 代表需要估计值的位置，纵坐标 yi 表示插值后的估计值。在 MATLAB 7.0 中，一维插值有基于多项式的插值和基于快速傅立叶的插值两种类型。
一维多项式插值可以通过函数 interp1()来实现
一维插值可以指定的方法如下： （1）最邻近插值（Nearest neighbor interpolation，method=’nearest’），这种插值方法在已知 数据的最邻近点设置插值点，对插值点的数进行四舍五入。对超出范围的点将返回一个 NaN （Not a Number）。 （2）线性插值（Linear interpolation，method=’linear’）,该方法是未指定插值方法时 MATLAB7.0 默认采用的方法。该方法采用直线连接相邻的两点。对超出范围的点将返回一 个 NaN（Not a Number）。 （3）三次样条插值（Cubic spline interpolation，method=’spline’），这样方法采用三次样条 函数来获得插值点。在已知点为端点情况下，插值函数至少具有相同的一阶和二阶导数。 （4）分段三次厄米多项式插值（Piecewise cubic Hermite interpolation，method=’pchip’）。 （5）三次多项式插值（method=’cubic’），与分段三次厄米多项式插值相同。 （6）MATLAB5 中使用的三次多项式插值（method=’v5cubic’），该方法使用一个 3 次多项 式函数对已知数据进行拟合。 当选择一个插值方法时应该考虑方法的执行速度、占用内存大小和获得数据的平滑度。 以上方法的特点如下。
（1）最邻近插值：最快的插值方法，但是数据平滑方面最差，其得到的数据是不连续的。 （2）线性插值：比邻近插值占用更多的内存，执行速度也稍慢，但其数据平滑方面优于 邻近插值。与邻近插值不同的是，线性插值的数据变化是连续的。 （3）三次样条插值：处理速度最慢，占用内存小于分段三次厄米多项式插值，可以产生 最光滑的结果。但是如果输入数据不均匀或者某些点靠得很近，会出现一些错误。样条插值 是极其有用的插值方法，因此除了提供三次样条插值函数外，MATLAB 7.0 还提供了一个样 条插值工具箱，它位于 toolbox\splines 下。 （4）分段三次厄米多项式插值：处理速度和消耗的内存比线性插值差，但插值得到的数 据和一阶导数都是连续的。 这些方法的相对优劣不仅适用于一维插值，而且适用于二维插值情况甚至高维插值 情况。
下面例子用不同插值方法对一维数据进行插值，并比较其不同，代码设置如下
%interp1_example.m %用不同插值方法对一维数据进行插值，并比较其不同
x = 0:1.2:10; 
y = sin(x); 
xi = 0:0.1:10;  yi_nearest = interp1(x,y,xi,’nearset’); %最邻近插值 yi_linear = interp1(x,y,xi);            %默认插值方法是线性插值 yi_spline = interp1(x,y,xi,’spline ’);  %三次样条插值 yi_cubic = interp1(x,y,xi,’cubic’);     %三次多项式插值 yi_v5cubic = interp1(x,y,xi,’v5cubic’); %MATLAB 5 中使用的三次多项式插值
hold on;
subplot(2,3,1);
plot(x,y,’ro’,xi,yi_nearest,’b-’); title(’最邻近插值’);
subplot(2,3,2);
plot(x,y,’ro’,xi,yi_linear,’b-’); title(’线性插值’);
subplot(2,3,3);
plot(x,y,’ro’,xi,yi_spline,’b-’); title(’三次样条插值’);
subplot(2,3,4);
plot(x,y,’ro’,xi,yi_cubic,’b-’); title(’三次多项式插值’);
subplot(2,3,5);
plot(x,y,’ro’,xi,yi_v5cubic,’b-’); title(’三次多项式插值(MATLAB5)’);

一维快速傅立叶插值
一维快速傅立叶插值通过函数 interpft()来实现，该函数用傅立叶变换把输入数据变换到 频域，然后用更多点的傅立叶逆变换，变换回时域，其结果是对数据进行增采样
实例：利用一维快速傅立叶插值实现数据增采样，代码设置如下：
%interpft_example.m %一维快速傅立叶插值实现数据增采样
x = 0:1.2:10; 
y = sin(x);  n = 2*length(x);             %增采样 1 倍 yi = interpft(y,n);            %一维快速傅立叶插值
xi = 0:0.6:10.4;               
hold on;                       
plot(x,y,’ro’);                %画图  
plot(xi,yi,’b.-’);
title(’一维快速傅立叶插值’);
legend(’原始数据’,’插值结果’);

二维插值
二维插值主要应用于图像处理和数据的可视化，其基本思想与一维插值相同，它是对两 变量的函数 z=f(x,y)进行插值
MATLAB 7.0 中的二维插值函数为 interp2()，其调用格式如下。 • zi = interp2(x,y,z,xi,yi)，原始数据 x,y,z 决定插值函数 z=f(x,y)，返回值 zi 是(xi,yi)在函数 f(x,y)上的值； • zi = interp2(z,xi,yi)，若 z=n×m，则 x=1:n，y=1:m； • zi = interp2(z,ntimes)，在两点之间递归地插值 ntimes 次； • zi = interp2(x,y,z,xi,yi,method)，可采用的插值方法在本小节将详细介绍； • zi = interp2(...,method, extrapval)，当数据超过原始数据范围时，用输入 extrapval 指定 一种外推方法。 二维插值可以采用的插值方法如下： （1）最邻近插值（Nearest neighbor interpolation，method=’nearest’），这种插值方法在已知 数据的最邻近点设置插值点，对插值点的数进行四舍五入。对超出范围的点将返回一个 NaN （Not a Number）。 （2）双线性插值（Bilinear interpolation，method=’linear’）,该方法是未指定插值方法时 MATLAB7.0 默认采用的方法。插值点的值只决定于最邻近的 4 个点的值。 （3）三次样条插值（Cubic spline interpolation，method=’spline’），这样方法采用三次样条
函数来获得插值数据。
实例：）双三次多项式插值（method=’cubic’）。 例如，下面例子比较各种二维插值插值算法的不同，如下：
%interp2_example2 %采用二次插值对三维高斯型分布函数进行插值 [x,y] = meshgrid(-3:0.8:3);                  %原始数据
z = peaks(X,Y); [xi,yi] = meshgrid(-3:0.25:3);                %插值点 zi_nearest = interp2(x,y,z,xi,yi,’nearset’);       %最邻近插值 zi_linear = interp2(x,y,z,xi,yi);               %默认插值方法是线性插值 zi_spline = interp2(x,y,z,xi,yi,’spline ’);        %三次样条插值 zi_cubic = interp2(x,y,z,xi,yi,’cubic’);          %三次多项式插值
hold on;
subplot(2,3,1);
surf(x,y,z); title(’原始数据’);
subplot(2,3,2);
surf(xi,yi,zi_nearest); title(’最邻近插值’);
subplot(2,3,3);
surf(xi,yi,zi_linear); title(’线性插值’);
subplot(2,3,4);
surf(xi,yi,zi_spline); title(’三次样条插值’);
subplot(2,3,5);
surf(xi,yi,zi_cubic); title(’三次多项式插值’); figure;                                  %新开绘图窗口 subplot(2,2,1);                           %画插值结果的等高线
contour(xi,yi,zi_nearest); title(’最邻近插值’);
subplot(2,2,2);
contour(xi,yi,zi_linear); title(’线性插值’);
subplot(2,2,3);
contour(xi,yi,zi_spline); title(’三次样条插值’);
subplot(2,2,4);
contour(xi,yi,zi_cubic);
title(’三次多项式插值’);

数据分析和傅立叶变换
MATLAB7.0 中与数据分析和傅立叶变换相关的函数位于目录\toolboxs\MATLAB\datafun 下，本节将分类介绍这些函数。 MATLAB 7.0 在这些函数中的约定。 • 一维统计数据可以用行向量或者列向量来表示，不管输入数据是行向量或者是列向量， 运算是对整个向量进行的； • 二维统计数据可以采用多个向量表示，也可以采用二维矩阵来表示。对二维矩阵运算 时，运算总是按列进行的。 这两条约定不仅适用于本节提到的函数，而且适用于 MATLAB 7.0 各个工具箱的函数。
基本数据分析函数
MATLAB 7.0 提供的基本数据分析函数的功能和调用格式如表 6-3 所示。
表 6-3 基本数据分析函数 函数名 功能描述 基本调用格式
max 求最大值
C = max(A)，如果 A 是向量，返回向量中的最大值；如果 A 是矩阵 返回一个包含各列最大值的行向量 C = max(A,B)，返回矩阵 A 和 B 之中较大的元素，矩阵 A、B 必须 具有相同的大小 C = max(A,[],dim)，返回 dim 维上的最大值 [C,I] = max(...)，多返回最大值的下标 min 求最小值 与最大值函数 max()调用格式一致
mean 求平均值
M = mean(A)，如果 A 是向量，返回向量 A 的平均值，如果 A 是矩 阵，返回含有各列平均值的行向量 M = mean(A,dim)，返回 dim 维上的平均值 median 求中间值 与求平均值函数 mean()的调用格式一致
std 求标准方差
s = std(A)，如果 A 是向量，返回向量的标准方差；如果 A 是矩阵， 返回含有各列标准方差的行向量 s = std(A,flag)，用 flag 选择标准方差的定义式 s = std(A,flag,dim)，返回 dim 维上的标准方差
var 方差（标准方差的平反）
var(X)，如果 A 是向量，返回向量的方差；如果 A 是矩阵，返回含 有各列方差的行向量 var(X,1)，返回第二种定义的方差 var(X,w)，利用 w 做为权重计算反差 var(X,w,dim)，返回 dim 维上的方差
sort 数据排序
B = sort(A)，如果 A 是向量，升序排列向量；如果 A 是矩阵，升序 排列各个列 B = sort(A,dim)，升序排列矩阵 A 的 dim 维 B = sort(...,mode)，用 mode 选择排序方式：’ascend’为升序，’descend’ 为降序 [B,IX] = sort(...)，多返回数据 B 在原来矩阵中的下标 IX
sortrows 对矩阵的行排序
B = sortrows(A)，升序排序矩阵 A 的行 B = sortrows(A,column)，以 column 列数据作为标准，升序排序矩阵 A 的行 [B,index] = sortrows(A)，多返回数据 B 在原来矩阵 A 中的下标 IX
sum 求元素之和
B = sum(A)，如果 A 是向量，返回向量 A 的各元素之和；如果 A 是 矩阵，返回含有各列元素之和的行向量 B = sum(A, dim)，求 dim 维上的矩阵元素之和 B = sum(A, ’double’)，返回数据类型指定为双精度浮点数 B = sum(A, ’native’)，返回数据类型指定为与矩阵 A 的数据类型相同
prod 求元素的连乘积
B = prod(A)，如果 A 是向量，返回向量 A 的各元素连乘积；如果 A 是矩阵，返回含有各列元素连乘积的行向量 B = prod(A,dim)，返回 dim 维上的矩阵元素连乘积
hist 画直方图
n = hist(Y)，在 10 个等间距的区间统计矩阵 Y 属于该区间的元素个 数 n = hist(Y,x)，在 x 指定的区间统计矩阵 Y 属于该区间的元素个数 n = hist(Y,nbins)，用 nbins 个等间距的区间统计矩阵 Y 属于该区间的 元素个数 hist(...)，直接画出直方图
histc 直方图统计
n = histc(x,edges)，计算在 edges 区间内向量 x 属于该区间的元素个 数 n = histc(x,edges,dim)，在 dim 维上统计 x 出现的次数
trapz 梯形数值积分（等间距）
Z = trapz(Y)，返回 Y 的梯形数值积分 Z = trapz(X,Y)，计算以 X 为自变量时 Y 的梯形数值积分 Z = trapz(...,dim)，在 dim 维上计算梯形数值积分
cumsum 矩阵的累加
B = cumsum(A)，如果 A 是向量，计算向量 A 的累计和；如果 A 是 矩阵，计算矩阵 A 在列方向上的累计和 B = cumsum(A,dim)，在 dim 维上计算矩阵 A 的累计和 cumprod 矩阵的累积 函数调用格式与函数 cumsum()相同 cumtrapz 梯形积分累计 函数调用格式与函数 trapz()相同
实例：．最大值、最小值、平均值、中间值、元素求和
这几个函数的用法都比较简单，下面的例子用来说明这些函数的用法：
%max_example %最大值、最小值、平均值、中间值、元素求和函数使用的例子
x=1:40; y=randn(1,40);                            %正态分布随机数据
hold on;                                 
plot(x,y); [y_max,I_max]=max(y);                     %求向量最大值和相应下标
plot(x(I_max),y_max,’o’); [y_min,I_min]=min(y);                     %求向量最小值和相应下标
plot(x(I_min),y_min,’*’);            y_mean=mean(y);                           %求向量平均值
plot(x,y_mean*ones(1,length(x)),’:’); y_median=median(y);                       %求向量反差
plot(x,y_median*ones(1,length(x)),’-.’); y_sum=sum(y);                             %求向量元素之和
plot(x,y_sum*ones(1,length(x)),’--’);
xlabel(’x’);
ylabel(’y’); legend(’数据’,’最大值’,’最小值’,’平均值’,’中间值’,’向量元素之和’);

MATLAB 7.0 中默认使用（5-1）式来计算数据的标准方差。要使用（5-2）式来计算标 准方差可以使用函数调用格式 std(A,1)。 方差是标准方差的平反。对应标准方差的两种定义，方差也有两种定义。同样 MATLAB 7.0 中默认使用（5-1）式来计算数据的方差。要使用（5-2）式来计算方差可以使用函数调用 格式 var(A,1)。 例如，有一维随机变量 x，要求验证 2*x 的方差是 x 的 4 倍，2x 标准方差是 x 的 2 倍， 代码设置如下：
%std_example.m %验证 2*x 的方差是 x 的 4 倍，2x 标准方差是 x 的 2 倍
x=rand(1,1000);
x_std=std(x);
x_var=var(x);
x2_std=std(2*x);
x2_var=var(2*x); disp([’随机变量 x2 的标准方差与 x1 的标准方差之比 = ’ num2str(x2_std/x_std)]); disp([’随机变量 x2 的标准方差与 x1 的标准方差之比 = ’ num2str(x2_var/x_var)]);

．元素排序
MATLAB 7.0 可以对实数、复数和字符串进行排序。对复数矩阵进行排序时，先按复数 的模进行排序，如果模相等则按其在区间[-π,π]上的相角进行排序。MATLAB7.0 中实现排 序的函数为 sort()。
实例：对复数矩阵排列，代码设置如下：
%sort_complex.m %对复数进行排序
a = [0 -1 -0.5i 1 0.5 i -i -0.5 0.5i];
b = sort(a)
MATLAB 7.0 提供函数 sortrows()用于对矩阵的行进行排列。
实例：%minvar.m %x 是一个向量，k 是小于等于向量 x 长度的正标量 %返回值是 y 是 x 中反差最小的 k 个数据，i 是 y 的元素在 x 向量中的下标
function [y,i]=minvar(x,k)
if(isvector(x)==0),                              %输入数据检查
    disp('输入数据 x 必须是向量!');
    return;
end;
if(isscalar(k)==0)     disp('输入数据 k 必须是标量!');
    return;
end;
if(k<=0||k>length(x)),
    disp('输入数据 k 必须是是小于等于向量 x 长度的正标量!');
    return;
end; k=floor(k);                                     %把 k 截断为整数
[x_sort index]=sort(x);                            %对数据 x 排序
n=length(x);                                    %数据 x 的长度
for i=1:(n-k+1),                                  %对排好序的数据求所有相邻 k 个值
     std_x_sort(i)=std(x_sort(i:(i+k-1)));              %的标准方差
end;
[min_std_x_sort j]  =  min(std_x_sort);             %求最小的相邻 k 个值的标准方差和下标
i=index(j:(j+k-1));                                %求最小标准方差的相邻 k 个值在向量
x                                                %中的方差
y=x(i);                                         %求最小标准方差的 k 个值
调用实例：
%minvar_example.m %验证已知 n 个数据中挑选出方差最小的 k 个数据的函数 minvar()
data = (1:10).^2;              %数列 1 4 9 16 ...
index = randperm(10);          %1:10的随机排列
x = data(index)                %随机扰乱数据 data 的排列
[y i]=minvar(x,4)              %求数据 x 中方差最小的 4 个值

协方差和相关系数矩阵
在数理统计中，协方差和相关函数是极其重要的随机变量的数字特征
%corrcoef_example.m %计算矩阵的相关系数
x = randn(10000,4);                       % 4个完全无关随机变量
y(:,1) = x(:,1);
y(:,2) = x(:,2);             
 y(:,3) = 0.3*x(:,3) + 0.7*x(:,4);              %在第 3 个随机变量和第 4 个随机变量之间
y(:,4) = 0.2*x(:,3) + 0.8*x(:,4);              %引入相关性
[r,p] = corrcoef(y)                          %计算相关矩阵和不相关概率矩阵
[i,j] = find(p<0.05);                       %寻找相关的随机变量对
[i,j]                                    %显示相关的随机变量对的编号

有限差分和梯度
在 MATLAB7.0 中使用函数 diff()来计算差分，其调用格式如下。
实例：利用有限差分计算正弦函数的导数，代码设置如下：
%sin_diff.m %正弦函数的导数
x=0:0.1:10;                        %自变量
y=sin(x);                          %正弦函数
y_der=diff(y)./diff(x);                %导数等于有限差分之比(dy/dx)
hold on;         
x_der=x(1:(end-1));                 %导数的自变量
plot(x,y,'b-');                       %画图
plot(x_der,y_der,'b-.'); 
axis([0 10 -1.2 1.4]);                 %设定坐标轴范围
legend('正弦函数','正弦函数的导数');  %添加图例

在 MATLAB 7.0 中，用函数 gradient()来计算梯度
，计算二维高斯函数的梯度场，代码设置如下
%gradient_example.m %计算二维高斯函数的梯度场
v=-2:0.25:2;                                 
[x,y]=meshgrid(v,v);                    %产生自变量 x,y
z= exp(-(x.^2+y.^2+0.5*x.*y));           %二维高斯函数
[px py]=gradient(z,0.25);                %梯度场
contour(v,v,z,4);                       %画 4 条等高线
hold on;                               
quiver(v,v,px,py);                      %画出梯度场
图中椭圆是二维高斯函数的等高线，箭头的方 向代表梯度的方向，箭头的大小代表梯度的模。由图可见梯度方向总是垂直于等高线的，这 是梯度场的基本性质之一。 
信号滤波和卷积
MATLAB 7.0 中有关信号滤波和卷积的函数，
函数名 功能描述 filter 一维数字滤波器 filter2 二维数字滤波器（将在本书的第 9 章详细介绍） conv 一维卷积和多项式乘法 conv2 二维卷积（将在本书的第 9 章详细介绍） convn N 维卷积 deconv 反卷积和多项式除法 detrend 去除信号中的直流或者线形成分，主要用于 FFT 运算中
一维数字滤波
MATLAB 7.0 的一维数字滤波可以用函数 filter()来实现，它是直接 II 型线性差分方程组 的解，可以用于实现 FIR 和 IIR 滤波
如何得到滤波器的系数 b 和 a 涉及到滤波器设计问题，具体将在本书的第 9 章作出详细 介绍。在有些情况下，可以很容易得到滤波器的系数 b 和 a。例如，一个 5 阶平均值滤波器： y(n)=1/5*x(n) + 1/5*x(n-1) + 1/5*x(n-2) + 1/5*x(n-3) + 1/5*x(n-4)，其系数 a 可以表示为 1，系
数 b 可以表示为[1/5 1/5 /15 1/5 1/5]。 下面通过一个示例来实现对带噪声的正弦信号进行平均值滤波，代码设置如下：
%filter_example.m %对带噪声的正弦信号进行平均值滤波
t=0:0.1:10;                             %时间
n = 6*randn(size(t));                     %高斯白噪声
x = 40*sin(t)+n;                         %在正弦信号中添加噪声
a = 1;                                 %频率值滤波器的系数
b = [1/5 1/5 1/5 1/5 1/5];
y=filter(b,a,x);                          %滤波
plot(t,x,'b-');                            %画原始信号
hold on;
plot(t,y,'r:');                            %画滤波后的信号
axis([0 10 -60 65]);                      %设置坐标轴范围
xlabel('时间/ts(s)');
legend('原始数据','滤波后数据');          %添加图例

．信号卷积
向量 u 和 v 的卷积为 w，则记为w u v = ⊗ 。如果 u，v 的长度分别为 m 和 n，则 w 的长 度为 m+n-1。
从信号中去除直流成分和线性成分，代码设置如下：
%detrend_example.m %从信号中去除直流成分和线性成分
t = 0:0.04:5;                                      %时间
x = 2*t + 0.5*randn(size(t));                         %带线性成分的随机信号
x_no_linear=detrend(x);                            %去除去除线性成分
subplot(2,1,1);                                    %画图
hold on;
plot(t,x,'b-');
plot(t,x_no_linear,'b:'); axis([0 5 -2 14]);                                  %设置坐标轴范围
title('从信号中去除线形成分');                     
legend('原始数据','去除线性成分的数据');            %添加图例
y = 3 + 0.5*randn(size(t));                          %带直流成分的随机信号
y_no_constant = detrend(y,'constant');                 %去除直流成分
subplot(2,1,2);
hold on;
plot(t,y,'b-');
plot(t,y_no_constant,'b:');
axis([0 5 -2 8]);                                   %设置坐标轴范围
title('从信号中去直流成分'); 
legend('原始数据','去除直流成分的数据');             %添加图例

傅立叶变换
对信号进行分析通常可以在时域中进行，也可以在频域中进行，时域分析方法和频域分 析方法各有其优缺点。傅立叶变换是把信号从时域变换到频域，因此它在信号分析中占有极 其重要的地位，如滤波器设计、频谱分析等方面。 傅立叶变换既可以对连续信号进行分析，也可以对离散信号进行分析。连续信号的傅立
叶变换实际上要计算傅立叶积分，这部分内容可以参见本书第 7 章
MATLAB 7.0 提供的有关傅立叶变换的函数如表 6-5 所示。
函数名 功能描述 fft 一维离散快速傅立叶变换 fft2 二维离散快速傅立叶变换 fftn N 维离散快速傅立叶变换 ifft 一维离散快速傅立叶逆变换 ifft2 二维离散快速傅立叶逆变换 ifftn N 维离散快速傅立叶逆变换 fftshift 把频谱的直流成分移到中间 nextpow2 返回大于或等于参数值的二次幂，满足 2p>=abs(A)的最大整数
在 MATLAB 7.0 中一维离散傅立叶变换由函数 fft()来实现
在MATLAB 7.0中一维离散傅立叶逆变换由函数ifft()来实现，其调用方式基本与函数fft() 相同，只是添加一个选项，其调用方式如下。 • y = ifft(..., ’symmetric’)，把频谱 X(k)看成是共轭对称。该选项在 X(k)不是严格共轭对称 时有效； • y = ifft(..., ’nonsymmetric’)，该选项为默认选项
与傅立叶变换函数经常使用的函数是 fftshift()函数。如图 6-19 所示，信号 x(n)经过 FFT 后得到频谱 X(k)，X(k)的头部和尾部代表信号的低频分量大小，X(k)的中间则代表高频部分。 习惯上画频率响应曲线时把直流成分放在曲线的中间，而把高频部分放在曲线的头部和尾部。 因此频谱 X(k)不可以直接用于画频率响应曲线。为了画图时方便，MATLAB 7.0 提供函数 fftshift()用于把 FFT 变换的结果频谱 X(k)转换为适合画图显示的格式。对一维信号操作时， fftshift()函数把信号分为左半部分和右半部分，然后交换其左右部分。 
，计算单位冲激信号经过一个带阻滤波器前后的频谱，代码设置如下：
%fliter_example.m %单位冲击信号通过带阻滤波器
t = 1:40;                                   %信号的时间
x =zeros(size(t));                         
x(1) = 1;                                   %产生单位冲击信号
>> [b,a] = butter(10,[0.3 0.7],'stop');                %设计带阻滤波器
y = filter(b,a,x);                             %滤波 hold on;                                   %画图滤波前后的时域数据
stem(t,x,'marker','o');
stem(t,y,'marker','.');
xlabel('时间/ts(s)');
legend('原始数据','滤波后数据'); 
fx=fft(x);                                  %对单位冲击信号进行傅立叶变换  
fx=fftshift(fx); fy=fft(y);                                  %滤波后的信号进行傅立叶变换
fy=fftshift(fy);
figure; subplot(2,1,1);                             %画图显示滤波前后的数据频谱
f=(t-20)/20; plot(f,abs(fx),'b-',f,abs(fy),'b-.');               %画信号的幅频曲线
xlabel('数字频率/\pi(rad)'); title('幅频曲线');
legend('原始数据的幅频曲线谱','滤波后数据的幅频曲线');
subplot(2,1,2); 
 plot(f,angle(fx),'b-',f,angle(fy),'b-.');           %画信号的相频曲线
xlabel('数字频率/\pi(rad)'); title('相频曲线'); legend('原始数据的相频曲线','滤波后数据的相频曲线');
由图 6-21 可见，单位冲激信号的频谱是一个常数，它通过一个带阻滤波器后的频谱等于 带通滤波器的频率响应曲线。这个结果验证了一个信号通过一个线性系统后的频谱等于信号 频谱乘以线性系统的频率响应的结论。

二维傅立叶变换
在图像处理中经常使用二维傅立叶变换来进行图像滤波等操作。在 MATLAB7.0 中二维 傅立叶变换用函数 fft2()来实现。函数 fft2()可以看成是对信号进行两次一维傅立叶变换，即 fft2(X)=fft(fft(X).’).’。函数 fft2()的调用格式如下：
实例：分析图 6-22 的频谱，代码设置如下
%fft2_example.m %分析图像的频谱
img = imread('coins.png');                            %读图像文件
f_img = fft2(double(img));                            %二维 fft 变换
f_img = fftshift(f_img);                               %把直流成分移动频谱的中间
imshow(img);                                      %显示图像
figure; f_img_abs = abs(f_img);                             %得到频谱幅度
f_img_abs = (f_img_abs-min(min(f_img_abs)))./(max(max(f_img_abs))-min(min(f_img_abs)))*255;  %把频谱幅度变换到[0,255]范围内
imshow(f_img_abs);                                %显示频谱幅度
title('图像的幅频分布');

功能函数 自写
功能函数就是可以用其他函数作为输入变量的函数，例如一个可以画函数图像的函数 fplot()。其调用格式为 fplot(@fun,[-pi,pi])。输入变量@fun 为需要画成图的函数的函数句柄。 MATLAB 7.0 的功能函数都在目录“\toolbox\matlab\funfun”下。
在 MATLAB 7.0 中函数可以用 M 文件、匿名函数和 inline 对象来表示。其中匿名函数是
MATLAB 7.0 的新特性。
%customized_fun.m %用户自己定义的函数
function y=customized_fun(x)
y=2./(1+exp(-x))+3./(1+exp(-2*x));
把 customized_fun.m 文件放在 MATLAB 7.0 的路径下，例如“\work”目录下，就可以使 用该函数了。例如在命令行中输入如下代码：
customized_fun(magic(4))

．一元函数最小值
一元函数在给定区间内的最小值问题可以用函数 fminbnd()来实现
实例：求正弦函数在[0,10]内的最小值，代码设置如下：
%fminbnd_example1.m %求正弦函数在[0 10]内的最小值
x=fminbnd(@sin,0,10)
如果要显示计算函数最小值的过程，则可以使用 options 参数来设定，代码设置如下
%fminbnd_example2.m %求正弦函数在[0 10]内的最小值
x=fminbnd(@sin,0,10,optimset(’Display’,’iter’))

求多元函数的最小值
求多元函数的最小值用函数 fminsearch()来实现。使用该函数时必须指定开始矢量 X0， 此函数返回一个在矢量 X0 附近的局部最小值
实例：求二维函数 2 2 2 2 1 1 f( ) 100( ) (1 ) x x x x = − + − 的局部最小值，代码设置如下
%fminsearch_example.m %寻找二维最小值
banana = @(x)100*(x(2)-x(1)^2)^2+(1-x(1))^2;
[x,fval] = fminsearch(banana,[-1.2, 1])

求一元函数的零点
可以使用函数 fzero()来求一元函数的零点。寻找一元函数零点时，可以指定一个开始点， 或者指定一个区间。当指定一个开始点时，此函数在开始点附近寻找一个使函数值变号的区 间，如果没有找到这样的区间，则函数返回 NaN。如果知道函数零点所在的区间，则可以使 用一个包含两个元素的向量来指定此区间。fzero()函数的调用格式如下。
用函数 fzero()来解一元非线性方程
2 2 1 1 1 ( 4) 1 ( 4) 1 2 x x + = + + − +
。该问题等价
于求函数
2 2 1 1 1 f( ) ( 4) 1 ( 4) 1 2 x x x = + + + − +
的零点，函数 f(x)的曲线图形如图 6-25 所示
%fzero_example.m %求解一元非线性方程
f = @(x) 1./((x+4).^2+1)+1./((x-4).^2+1)-0.5;             %用匿名函数来表示函数
fplot(f,[-10 10]);                                    %画出函数的图
xlabel('x');
ylabel('f(x)');
 x1 = fzero(f,-5)                                     %求某个点附近的零点
x2 = fzero(f,-3)
x3 = fzero(f,3)
x4 = fzero(f,5)
x1_region = fzero(f,[-10,-3])                          %求某个区间内的零点

．一元函数的数值积分
MATLAB 7.0 提供了两个函数 quad()和 quadl()计算一元函数的积分。函数 quad()采用低 阶的自适应递归 Simpson 方法，函数 quadl()则采用高阶的自适应 Lobatto 方法
实例：，求归一化高斯函数的在区间[-1 1]上的定积分，代码设置如下：
%quad_exam.m %求归一化高斯函数的在区间[-1 1]上的定积分,并得到积分过程的中间结点 y=@(x)1/sqrt(pi)*exp(-x.^2);                %归一化高斯函数 quad(y,-1,1,2e-6,1)                        %求定积分，并显示中间迭代过程 fplot(y,[-1 1],’b’);                          %画出函数
hold on; %跟踪数据（运行完上面程序后，可以在命令行中复制这些数据）
trace = [ 9    -1.0000000000    5.43160000e-001     0.1804679399;
      11    -1.0000000000    2.71580000e-001     0.0728222057;
      13    -0.7284200000    2.71580000e-001     0.1076454255;
      15    -0.4568400000    9.13680000e-001     0.4817487615;
      17    -0.4568400000    4.56840000e-001     0.2408826755;
      19    -0.4568400000    2.28420000e-001     0.1142172651;
      21    -0.2284200000    2.28420000e-001     0.1266655031;
      23     0.0000000000    4.56840000e-001     0.2408826755;
      25     0.0000000000    2.28420000e-001     0.1266655031;
      27     0.2284200000    2.28420000e-001     0.1142172651;
      29     0.4568400000    5.43160000e-001     0.1804679399;
      31     0.4568400000    2.71580000e-001     0.1076454255;
      33     0.7284200000    2.71580000e-001     0.0728222057]; x1 = trace(:,2);                           %积分过程的中间结点 y1 = y(x1);                              %中间结点的函数值 plot(x1,y1,’ro’);                           %画图
xlabel(’x’);
ylabel(’y’); legend(’高斯函数’,’求积分过程的中间结点’);
函数 quadl()的调用格式与函数 quad()相同，可以参见上面 quad()函数的用法。

矢量积分
矢量积分相当于多个一元定积分
实例：
%quadv_exam.m %求矢量积分
y=@(x,n)1./(sqrt(2*pi).*(1:n)).*exp(-x.^2./(2*(1:n).^2));     %归一化高斯函数
quadv(@(x)y(x,5),-1,1)

MATLAB 7.0 中使用函数 dblquad()来计算二重积分。该函数先计算内积分值，然后利用 内积分的中间结果来计算二重积分。根据 dxdy 的顺序，称 x 为内积分变量 ， y 为外积分变量。 函数 dblquad()的基本格式如下
计算二维高斯函数在矩形区间[-1 1 -1 1]上的二重积分，代码设置如下:
%quadv_exam.m %求矢量积分
y=@(x,n)1./(sqrt(2*pi).*(1:n)).*exp(-x.^2./(2*(1:n).^2));     %归一化高斯函数
quadv(@(x)y(x,5),-1,1)
函数 dblquad()处理的都是矩形积分区域。若要计算非矩形的积分区间的二重积分，可以 用一个大的矩形积分区域包含非矩形的积分区间，然后在非矩形的积分区间之外的区域上把 二元函数的值取零。
实例：，计算二维高斯函数在圆形区域 2 2 1 x y + < 上的二重积分，代码设置如下
%dblquad_exam2.m %计算二维高斯函数在圆形区域 sqrt(x.^2+y.^2)<1 上的二重积分 %在圆形区域外填充 0 的归一化高斯函数
f=@(x,y)(1/sqrt(pi)*exp(-x.^2)*1/sqrt(pi)*exp(-y.^2)).*(sqrt(x.^2+y.^2)<=1);   
dblquad(f,-2,2,-2,2,1e-6,@quadl)
三重积分可以用函数 triplequad()来实现，其用法与二重积分类似。

实例：求函数f( ) xx e a x b = + * - 的零点：
%fzero_nestedfun.m %用嵌套函数的方法实现处理含参函数
function y = fzero_nestedfun(a,b,x0) %a,b 是含参函数的参数值 %x0 是求函数零点的开始点
options = optimset('Display', 'off');       %关闭显示
y = fzero(@para_fun,x0,options)        %求函数零点    
function y=para_fun(x)            %含参函数
        y=exp(x)+a*x-b;
   end
end
实例：a=1,b=[-10 10]时，含参函数的零点取值随着参数 b 变化的关系，可以通过如下程 序来实现：
%fzero_nestedfun_exam.m %使用含参函数的例子
clear; a=1;                               %参数取值
b=-10:10;
x = zeros(size(b));
for i=1:length(b),            
                 x(i) = fzero_nestedfun(a,b(i),0);     %求函数零点
end; plot(b,x);                           %画图
xlabel(’b’); ylabel(’零点’);

用匿名函数提供函数参数
使用含参函数还可以利通过匿名函数来实现。由于与函数 fzero()一样，一般功能函数只 接受一个自变量，因此函数的参数在使用之前必须先赋值。具体步骤如下
%指数函数与线性函数的交点
function y=exp_linear(x,a,b)
    y=exp(x)+a*x-b;
然后，用如下代码来实现计算函数零点：
%fzero_anonymousfun.m %用匿名函数提供函数参数
clear; a=1;                                      %参数取值
b=-10:10;
x = zeros(size(b));
for i=1:length(b),  
    f = @(x) exp_linear(x,a,b(i));
    x(i) = fzero(f,0);                        %求函数零点
end; plot(b,x);                                  %画图
xlabel(’b’); ylabel(’零点’);

微分方程组数值解
在 MATLAB 7.0 中可以计算 3 类微分方程数值解问题，即常微分方程组的初值问题、延 迟微分方组和常微分方程组的边界问题。
常微分方程组的初值问题
MATLAB 7.0 可以解决 3 种常微分方程组的初值问题，即显式常微分方程组、线性隐式 常微分方程组和完全隐式常微分组

显式常微分方程组的解法
由于显式常微分方程组的多样性，需要根据显式常微分方程组的特点采用相应的解法来 处理。MATLAB7.0给出7种解法，分别由7个不同的函数来实现，它们是函数ode45()、ode23()， ode113()、ode15s()、ode23s()、ode23t()和 ode23tb()。 在介绍这些函数之前，首先介绍一个概念——刚性(stiffness)。定性地讲，对一个常微分 方程组，如果其 Jocabian 矩阵的特征值相差十分悬殊，那么这个方程组就称为刚性方程组。 对于刚性方程组，为了保持解法的稳定，步长选取很困难。有些解法不能用于解刚性方程组， 而有些解法对方程组的刚性要求不严，可以用于解刚性问题。表 6-10 对 MATLAB 7.0 中的 7 个解法进行了对比

表 6-10 常微分方程组解法对比 函数名 采用算法 精度 适用系统 优缺点
ode45()
四阶/五阶龙格-库 塔
中 非刚性方程
属于单步算法（只需要前一步的解即可计算 出当前的解），不需要附加初始值，因而， 计算过程中随意改变步长也不会增加任何 计算量。通常 ode45()对很多问题来说都是 首选的最好方法
ode23()
二阶/三阶龙格-库 塔
低 非刚性方程
属于单步算法，在误差容许范围较宽或者存 在轻微刚度时性能比 ode45()好
ode113()
可 变 阶 Adams PECE 算法
低~高 非刚性方程
属于多步解法（需要前几步的解来计算当前 解）。比 ode45()更适合解决误差允许范围比 较严格的情况
ode15s()
可变阶的数值微分 公式算法（NDFS）
低~中 刚性方程
属于多步解法。如果 ode45()解法速度很慢， 很可能系统是刚性的，可以尝试采用该算法
第 6 章  数据分析
–237–
续表
函数名 采用算法 精度 适用系统 优缺点
ode23s()
基 于 改 进 的 Rosenbrock 公式
低 刚性方程
属于单部解法。比 ode15s()更适用于误差容 许范围较宽的情况。可以解决一些 ode15s() 效果不好的刚性方程
ode23t()
自由内插实现的梯 形规则
低 轻微刚性方程 适用于轻微刚性系统，给出的解无数值衰减
ode23tb
TR-BDF2 方法，即 龙格-库塔公式的第 一级采用梯形规则、 第二级采用 Gear 法
低 刚性方程
对于误差允许范围比较宽的情况，比 ode15s 效果好

解非刚性微分方程实例：
把该一阶常微分方程组用一个 M 文件形式的函数来表述，代码设置如下： 表示一个微分方程的函数必须有时间 t 这个输入变量，然而在函数内部计算时可以不使用 t 这个变量。
%ivpodefun.m %常微分方程
function dydt = ivpodefun(t,y)
dydt    = zeros(2,1);
dydt(1) = y(2); 
dydt(2) = (1-y(1)^2)*y(2)-y(1);
最后，用函数 ode45()来解这微分方程，并画图计算结果，代码设置如下：
%ode45_example.m %解常微分方程
[t,y] = ode45(@ivpodefun,[0 20],[2; 0]);
plot(t,y(:,1),’-’,t,y(:,2),’--’)
title(’常微分方程的解’);
xlabel(’t’);
ylabel(’y’); legend(’y’,’y 的一阶导数’);
如果要处理含参数的微分方程组可以把参数作为表示微分方程组的函数的第 3 个输入变 量。然后再调用函数 ode45()时，把参数值作为最后一个输入变量。其微分方程组可以表示为：
%ivpodefun_para.m %常微分方程
function dydt = ivpodefun_para(t,y,u)
dydt    = zeros(2,1);
dydt(1) = y(2); 
dydt(2) = u*(1-y(1)^2)*y(2)-y(1);
如果要求参数 u=10 的微分方程组，可以用如下代码实现：
[t,y] = ode45(@ivpodefun_para,[0 20],[2; 0],[],10)

设置常微分方程组解法器参数
正如上小节介绍的，在调用常微分方程组函数时可以输入一个 options 结构体，从而设置 解法器的参数。该 options 结构体可以用函数 odeset()来设定
常微分方程组解法器参数如表 6-11 所示。
表 6-11 常微分方程组解法器参数 参数名 可取值 默认值 用途描述 RelTol 正标量 1e-3 用于所有分量的相对误差，解法器的积分估计误差必 须小于相对误差与解的乘积并且小于绝对误差
AbsTol
正 标 量 或 者 向量 1e-6
绝对误差允许范围，如果是标量该绝对误差应用于所 有的分量；如果是向量则单独指定每一个分量的绝对 误差
NormControl ’on’或者’off’ ’off’
如果该值为’on’，解法器采用积分估计误差的模来控 制计算精度；如果该值为’off ’，解法器采用更加严格 的精度控制策略，即严格控制每一分量的计算精度
OutputFcn 函数句柄 @odeplot 或 者[]
每个时间步长计算完后，解法器调用该函数来输出。 如果调用 ode 类函数时，没有输出变量，则默认采用 odeplot 来输出数据；如果有输出变量，则默认值是[]， 即不输出结果。可选的输出函数有： odeplot：一维时域画图 odephas2：二维相位平面画图 odephas3：三维相位平面画图 odeprint：在命令行输出解
OutputSel 正整数向量
所有分量的 下标
OutputSel 向量包含的下标所对应的分量被送给输出 函数 OutputFcn 输出。默认情况下，所有分量都输出
Refine 正整数 1 或 者 4 （ode45）
如果 refine 大于 1，则输出结果被插值，从而提供输 出结果的精度 Stats ’on’或者’off’ ’off’ 如果该值为’on’，输出计算耗费时间，否则不输出
Jacobian
函 数 或 者 常 数矩阵
无 指定常微分方程的 Jacobian 矩阵
JPattern 稀疏矩阵 无
给出微分方程的 Jacobian 矩阵稀疏样式，如果微分方 程的 Jacobian 矩阵的元素不为 0 则 JPattern 相应元素 为 1，否则为 0
Vectorized ’on’或者’off’ ’off’
ode 函数是否被向量化，如果该值等于 1，则 ode 函 数 x 形式为[f(t,y1) f(t,y2) ...]，否则返回格式为 f(t,[y1 y2 ...]) Events 函数 无 定位事件
Mass
常 数 矩 阵 或 者函数
无 指定线性隐式常微分方程组的加权函数 M(t,y)
MStateDependence
’none’ ’weak’ 或者’strong’ ’weak’ 说明加权函数 M(t,y)是否依赖于 y MvPattern 稀疏矩阵 无 ( ( , ) )/ M t y v y抖 的稀疏矩阵样式。用于加权矩阵 M 与 y 强相关的情况
MassSingular
’yes’ ’no’或者 ’maybe’
’maybe’ 加权函数 M(t,y)是否奇异 InitialSlope 向量 无 初始一阶导数值 yp0，满足 M(t0,y0)*yp0 = f(t0,y0) MaxStep 正标量 自动选择 最大步长值 InitialStep 正标量 自动选择 解法器自动选择初始步长
MaxOrder
1，2，3，4， 5
5 ode15s 采用的最大阶数 BDF ’on’或者’off’ ’off’ 在 ode15s 算法中是否采用 BDF 算法
实例：设置解法器的输出函数为二维相位平面画图，代码设置如下：
%odeset_example.m %解常微分方程
option = odeset(’RelTol’,1e-6,’OutputFcn’,’odephas2’);
[t,y] = ode45(@ivpodefun,[0 20],[2; 0],option);
xlabel(’y’); ylabel(’y 的一阶导数’);

线性隐式常微分方程组
线性隐式常微分方程组可以利用 ode 类函数的解法器参数 options 来解决。首先创建 M 文件函数来表示加权函数 M(t,y)，然后用 odeset()函数来把加权函数设置到解法器的 Mass 参 数中，最后用 ode 类函数来解微分方程组。
，求解线性隐式常微分方程组 2 ( 1) ’ 2 2 y y y y + = + - ，首先创建表示微分方程的
函数 f(t,y)和 M(t,y)，表示 f(t,y)的函数，代码设置如下：
%odefun_LinImp.m %线性隐式常微分方程
function dydt =  odefun_LinImp(t,y)
dydt = y.^2 + 2*y -2; 表示 M(t,y)的函数，代码设置如下：
%odefun_LinImp_mass.m %线性隐式常微分方程
function mass =  odefun_LinImp_mass(t,y)
mass = y + 1; 求解微分方程，代码设置如下：
%ode_LinImp_example.m %解常微分方程
option = odeset(’RelTol’,1e-6,’OutputFcn’,’odeplot’,’Mass’,@odefun_LinImp_mass);
第 6 章  数据分析
–241–
[t,y] = ode45(@odefun_LinImp,[0 2],2,option);
xlabel(’t’);
ylabel(’y’);

完全隐式常微分方程组
MATLAB 7.0 的新特性之一是提供了一个完全隐式常微分方程组的解法器 ode15i()。其 调用格式如下：
首先创建一个代表该方程的 M 文件函数，代码设置如下：
%ode_weissfun.m %代表 Weissinger 方程的 M 文件函数
function exp = ode_weissfun(t,y,dydt)
exp = t*y.^2*dydt.^3-y.^3*dydt.^2+t*(t^2+1)*dydt-t^2*y;
解微分方程的程序，代码设置如下：
%ode15i_example.m %解完全隐式微分方程
t0 = 1;                                       %猜想的初值
y0 = sqrt(3/2);                     
yp0 = 0;                                              
[y0,yp0] = decic(@ode_weissfun,t0,y0,1,yp0,0);     %求出自洽初值，并保值 y0 值不变
[t,y] = ode15i(@ode_weissfun,[1 10],y0,yp0);       %求在[1 10]区间内的解
ytrue = sqrt(t.^2 + 0.5);                         %解析解
plot(t,y,t,ytrue,’o’);                             %画图比较数值解与解析解
xlabel(’t’);                                   
ylabel(’y’);    

延迟微分方程组数值解 在 MATLAB 7.0 中使用函数 dde23()来解延迟微分方程。其调用格式如下：
求解延迟微分方程组；
首先，确定延迟向量 lags，在本例中 lags=[1 3]。
其次，创建一个 M 文件形式的函数表示延迟微分方程组，代码设置如下：
%ddefun.m %延迟微分方程
function dydt=ddefun(t,y,Z)
dydt = zeros(2,1);
dydt(1) = Z(1,2).^2 + Z(2,1).^2;
dydt(2) = y(1) + Z(2,1); 
然后，创建一个 M 文件形式的函数表示延迟微分方程组的初始值，代码设置如下：
%ddefun_history.m %延迟微分方程的历史函数
function y=ddefun_history(t)
y = zeros(2,1);
y(1) = 1;
y(2) = t-2;
最后，用 dde23 解延迟微分方程组和用图形显示解，代码设置如下
%dde_example.m %解延迟微分方程的例子 lags = [1 3];                                    %延迟向量 sol = dde23(@ddefun,lags,@ddefun_history,[0,1]);   %解方程
hold on; plot(sol.x,sol.y(1,:),’b-’);                         %画出结果
plot(sol.x,sol.y(2,:),’r-.’); title(’延迟微分方程的解’);                            
xlabel(’t’);
ylabel(’y’);
legend(’y_1’,’y_2’,2);

常微分方程组的边界问题
常微分方程组的边界问题的形式如下：
同时要指定函数 y(x)在某些边界条件下的值，然后在边界条件下求函数 y(x)，其中 x 是 独立变量。MATLAB 7.0 中使用函数 bvp4c()来处理常微分方程组的边界问题。该函数解决第 一类边界条件问题，即边界条件是：
( ( ), ( )) 0 g y a y b = 其中 a 和 b 是求解区间的下界和上界，即要求解出 y 在区间[a b]上的值。 常微分方程组的边界问题与常微分方程组的初值问题的不同之处在于：常微分方程组的 初值问题总是有解的，然而常微分方程组的边界问题有时会出现无解的情况，有时会出现有 限个解，有时会出现无穷多个解。因此在解常微分方程组的边界问题时，不可或缺的一部分
工作是提供猜测解。猜测解决定了解边界问题的算法性能，甚至算法是否成功。 在边界问题中，还经常出现附加的未知参数，其方程形式如下：
dy
f(x,y,p)
dx
=
其边界条件为：
( ( ), ( ), ) 0 g y a y b p = 在这种情况下，边界条件必须充分，从而能够决定未知参数 p。 函数 bvp4c()的调用格式如下
实例：
求 Mathieu 方程的特征值。Mathieu 方程的形式为：
’’ ( 2 cos2 ) 0 y l q x y + - = 其中 q-1 是 Mathieu 方程的阶数，为已知参数。λ是未知参数 Mathieu 方程的特征值。 为在本例中 q=5，即求 4 阶 Mathieu 方程的特征值λ。首先把方程改写成一阶常微分方程的 形式：
1 2
2 1
’ ’ ( 2 cos2 ) y y y q x y l = = - 
指定边界条件为：
’(0) 0 ’( ) 0 (0) 1
y y y
p
= = =
其次，创建 M 文件来描述 Mathieu 方程，代码设置如下：
%bvp_Mathieufun.m %4 阶 Mathieu 方程
function dydx=bvp_Mathieufun(x,y,lambda)
q = 5;
dydx = zeros(2,1);
dydx(1) = y(2);
dydx(2) = -(lambda - 2*q*cos(2*x))*y(1);
然后，创建 M 文件来描述方程的初值，代码设置如下：
%Mathieu_initfun.m %4 阶 Mathieu 方程的初始值
function yinit = Mathieu_initfun(x)
yinit(1) =  cos(4*x); 
yinit(2) =  -4*sin(4*x);
还需要，创建 M 文件来描述方程的边界条件，代码设置如下：
%Mathieu_bcfun.m %4 阶 Mathieu 方程的边界条件
function res = Mathieu_bcfun(ya,yb,lambda)
res = zeros(3,1);
res(1) =  ya(2); 
res(2) =  yb(2); 
res(3) =  ya(1)-1;
最后，用函数 bvp4c 来解方程，并画图显示结果，代码设置如下：
%bvp4c_example.m %Mathieu 方程在边界条件下的解
lambda = 15;                                             %未知参数猜测解
solinit = bvpinit(linspace(0,pi,10),@Mathieu_initfun,lambda);%求初始值
sol = bvp4c(@bvp_Mathieufun,@Mathieu_bcfun,solinit);         %求解方程
xint = linspace(0,pi);                                       %画图
Sxint = deval(sol,xint);
plot(xint,Sxint(1,:))
axis([0 pi -1 1.1]) title(’Mathieu 方程在边界条件下的解’); 
xlabel(’x’)
ylabel(’y’) text(0.1,1,[’4 阶 Mathieu 方程的特征值 = ’ num2str(sol.parameters)]);

/matlab解方程，分析方程方法，函数///