注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

点点滴滴

感悟人生,享受精彩。

 
 
 

日志

 
 

第2章 数值计算与数据分析  

2009-03-30 15:49:06|  分类: 默认分类 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |
第2章  数值计算与数据分析
2.1  基本数学函数
2.1.1  三角函数与双曲函数
函数  sin、sinh
功能  正弦函数与双曲正弦函数
格式  Y = sin(X)   %计算参量X(可以是向量、矩阵,元素可以是复数)中每一个角度分量的正弦值Y,所有分量的角度单位为弧度。
Y = sinh(X)  %计算参量X的双曲正弦值Y
注意:sin(pi)并不是零,而是与浮点精度有关的无穷小量eps,因为pi仅仅是精确值π浮点近似的表示值而已;对于复数Z= x+iy,函数的定义为:sin(x+iy) = sin(x)*cos(y) + i*cos(x)*sin(y), ,
例2-1
x = -pi:0.01:pi; plot(x,sin(x))
x = -5:0.01:5; plot(x,sinh(x))
图形结果为图2-1。
    
图2-1  正弦函数与双曲正弦函数图
函数  asin、asinh
功能  反正弦函数与反双曲正弦函数
格式  Y = asin(X)   %返回参量X(可以是向量、矩阵)中每一个元素的反正弦函数值Y。若X中有的分量处于[-1,1]之间,则Y = asin(X)对应的分量处于[-π/2,π/2]之间,若X中有分量在区间[-1,1]之外,则Y= asin(X)对应的分量为复数。
Y = asinh(X)   %返回参量X中每一个元素的反双曲正弦函数值Y
说明  反正弦函数与反双曲正弦函数的定义为: ,
例2-2
x = -1:.01:1; plot(x,asin(x))
x = -5:.01:5; plot(x,asinh(x))
图形结果为图2-2。
    
图2-2  反正弦函数与反双曲正弦函数图
函数  cos、cosh
功能  余弦函数与双曲余弦函数
格式  Y = cos(X)   %计算参量X(可以是向量、矩阵,元素可以是复数)中每一个角度分量的余弦值Y,所有角度分量的单位为弧度。我们要指出的是,cos(pi/2)并不是精确的零,而是与浮点精度有关的无穷小量eps,因为pi仅仅是精确值π浮点近似的表示值而已。
Y = sinh(X)  %计算参量X的双曲余弦值Y
说明  若X为复数z= x+iy,则函数定义为:cos(x+iy) = cos(x)*cos(y) + i*sin(x)*sin(y), ,
例2-3
x = -pi:0.01:pi; plot(x,cos(x))
x = -5:0.01:5; plot(x,cosh(x))
图形结果为图2-3。
    
图2-3  余弦函数与双曲余弦函数图
函数  acos、acosh
功能  反余弦函数与反双曲余弦函数
格式  Y = acos(X)   %返回参量X(可以是向量、矩阵)中每一个元素的反余弦函数值Y。若X中有的分量处于[-1,1]之间,则Y = acos(X)对应的分量处于[0,π]之间,若X中有分量在区间[-1,1]之外,则Y = acos(X)对应的分量为复数。
Y = asinh(X)   %返回参量X中每一个元素的反双曲余弦函数Y
说明  反余弦函数与反双曲余弦函数定义为: ,
例2-4
x = -1:.01:1; plot(x,acos(x))
x = -5:.01:5; plot(x,acosh(x))
图形结果为图2-4。
    
图2-4  反余弦函数与反双曲余弦函数图
函数  tan、tanh
功能  正切函数与双曲正切函数
格式  Y = tan(X)   %计算参量X(可以是向量、矩阵,元素可以是复数)中每一个角度分量的正切值Y,所有角度分量的单位为弧度。我们要指出的是,tan(pi/2)并不是精确的零,而是与浮点精度有关的无穷小量eps,因为pi仅仅是精确值π浮点近似的表示值而已。
Y = tanh(X)  %返回参量X中每一个元素的双曲正切函数值Y
例2-5
x = (-pi/2)+0.01:0.01:(pi/2)-0.01;  % 稍微缩小定义域
plot(x,tan(x))
x = -5:0.01:5; plot(x,tanh(x))
图形结果为图2-5。
    
图2-5  正切函数与双曲正切函数图
函数  atan、atanh
功能  反正切函数与反双曲正切函数
格式  Y = atan(X)    %返回参量X(可以是向量、矩阵)中每一个元素的反正切函数值Y。若X中有的分量为实数,则Y = atan(X)对应的分量处于[-π/2,π/2]之间。
Y = atanh(X)   %返回参量X中每一个元素的反双曲正切函数值Y。
说明  反正切函数与反双曲正切函数定义为: ,
例2-6
x = -20:0.01:20; plot(x,atan(x))
x = -0.99:0.01:0.99; plot(x,atanh(x))
图形结果为图2-6。
    
图2-6  反正切函数与反双曲正切函数图
函数  cot、coth
功能  余切函数与双曲余切函数
格式  Y = cot(X)    %计算参量X(可以是向量、矩阵,元素可以是复数)中每一个角度分量的余切值Y,所有角度分量的单位为弧度。
Y = coth(X)    %返回参量X中每一个元素的双曲余切函数值Y
例2-7
x1 = -pi+0.01:0.01:-0.01;  %  去掉奇点x = 0
x2 = 0.01:0.01:pi-0.01;  % 做法同上
plot(x1,cot(x1),x2,cot(x2))
plot(x1,coth(x1),x2,coth(x2))
图形结果为图2-7。
    
图2-7  余切函数与双曲余切函数图
函数  acot、acoth
功能  反余切函数与反双曲余切函数
格式  Y = acot(X)   %返回参量X(可以是向量、矩阵)中每一个元素的反余切函数Y
Y = acoth(X)  %返回参量X中每一个元素的反双曲余切函数值Y
例2-8
x1 = -2*pi:pi/30:-0.1; x2 = 0.1:pi/30:2*pi; % 去掉奇异点x = 0
plot(x1,acot(x1),x2,acot(x2))
x1 = -30:0.1:-1.1; x2 = 1.1:0.1:30;
plot(x1,acoth(x1),x2,acoth(x2))
图形结果为图2-8。
    
图2-8  反余切函数与反双曲余切函数图
函数  sec、sech
功能  正割函数与双曲正割函数
格式  Y = sec(X)   %计算参量X(可以是向量、矩阵,元素可以是复数)中每一个角度分量的正割函数值Y,所有角度分量的单位为弧度。我们要指出的是,sec(pi/2)并不是无穷大,而是与浮点精度有关的无穷小量eps的倒数,因为pi仅仅是精确值π浮点近似的表示值而已。
Y = sech(X)  %返回参量X中每一个元素的双曲正割函数值Y
例2-9
x1 = -pi/2+0.01:0.01:pi/2-0.01;  % 去掉奇异点x = pi/2
x2 = pi/2+0.01:0.01:(3*pi/2)-0.01;
plot(x1,sec(x1),x2,sec(x2))
x = -2*pi:0.01:2*pi;
plot(x,sech(x))
图形结果为图2-9。
     
图2-9  正割函数与双曲正割函数图
函数  asec、asech
功能  反正割函数与反双曲正割函数
格式  Y = asec(X)   %返回参量X(可以是向量、矩阵)中每一个元素的反正割函数值Y
Y = asech(X)   %返回参量X中每一个元素的反双曲正割函数值Y
例2-10
x1 = -5:0.01:-1; x2 = 1:0.01:5;
plot(x1,asec(x1),x2,asec(x2))
x = 0.01:0.001:1; plot(x,asech(x))
图形结果为图2-10。
    
图2-10  反正割函数与反双曲正割函数图
函数  csc、csch
功能  余割函数与双曲余割函数
格式  Y = csc(X)    %计算参量X(可以是向量、矩阵,元素可以是复数)中每一个角度分量的余割函数值Y,所有角度分量的单位为弧度。
Y = csch(X)    %返回参量X中每一个元素的双曲余割函数值Y
例2-11
x1 = -pi+0.01:0.01:-0.01; x2 = 0.01:0.01:pi-0.01; % 去掉奇异点x=0
plot(x1,csc(x1),x2,csc(x2))
plot(x1,csch(x1),x2,csch(x2))
图形结果为图2-11。
    
图2-11  余割函数与双曲余割函数图
函数  acsc、acsch
功能  反余割函数与反双曲余割函数。
格式  Y = asec(X)    %返回参量X(可以是向量、矩阵)中每一个元素的反余割函数值Y
Y = asech(X)   %返回参量X中每一个元素的反双曲余割函数值Y
例2-12
x1 = -10:0.01:-1.01; x2 = 1.01:0.01:10; % 去掉奇异点x = 1
plot(x1,acsc(x1),x2,acsc(x2))
x1 = -20:0.01:-1; x2 = 1:0.01:20;
plot(x1,acsch(x1),x2,acsch(x2))
图形结果为图2-12。
    
图2-12  反余割函数与反双曲余割函数图
函数  atan2
功能  四象限的反正切函数
格式  P = atan2(Y,X)   %返回一与参量X和Y同型的、与X和Y元素的实数部分对应的、元素对元素的四象限的反正切函数阵列P,其中X和Y的虚数部分将忽略。阵列P中的元素分布在闭区间[-pi,pi]上。特定的象限将取决于sign(Y)与sign(X)。
例2-13
z=1+2i;
r = abs(z);
theta = atan2(imag(z),real(z))
z = r *exp(i *theta)
feather(z);hold on
t=0:0.1:2*pi;
x=1+sqrt(5)*cos(t);
y=sqrt(5)*sin(t);
plot(x,y);
axis equal; hold off
计算结果为:
theta =
1.1071
z =
1.0000 + 2.0000i
图形结果为图2-13。
图2-13  四象限的反正切函数图
2.1.2  其他常用函数
函数  fix
功能  朝零方向取整
格式  B = fix(A)   %对A的每一个元素朝零的方向取整数部分,返回与A同维的数组。对于复数参量A,则返回一复数,其分量的实数与虚数部分分别取原复数的、朝零方向的整数部分。
例2-14
   >>A = [-1.9, -0.2, 3.1415926, 5.6, 7.0, 2.4+3.6i];
   >>B = fix(A)
计算结果为:
   B =
      Columns 1 through 4
          -1.0000              0             3.0000             5.0000
      Columns 5 through 6
          7.0000             2.0000 + 3.0000i
函数  roud
功能  朝最近的方向取整。
格式  Y = round(X)   %对X的每一个元素朝最近的方向取整数部分,返回与X同维的数组。对于复数参量X,则返回一复数,其分量的实数与虚数部分分别取原复数的、朝最近方向的整数部分。
例2-15
>>A = [-1.9, -0.2, 3.1415926, 5.6, 7.0, 2.4+3.6i];
>>Y = round(A)
计算结果为:
Y =
     Columns 1 through 4
       -2.0000      0         3.0000         6.0000
     Columns 5 through 6
       7.0000     2.0000 + 4.0000i
函数  floor
功能  朝负无穷大方向取整
格式  B = floor(A)   %对A的每一个元素朝负无穷大的方向取整数部分,返回与A同维的数组。对于复数参量A,则返回一复数,其分量的实数与虚数部分分别取原复数的、朝负无穷大方向的整数部分。
例2-16
>>A = [-1.9, -0.2, 3.1415926, 5.6, 7.0, 2.4+3.6i];
>>F = floor(A)
计算结果为:
F =
    Columns 1 through 4
      -2.0000         -1.0000       3.0000         5.0000
    Columns 5 through 6
      7.0000      2.0000 + 3.0000i
函数  rem
功能  求作除法后的剩余数
格式  R = rem(X,Y)   %返回结果X - fix(X./Y).*Y,其中X、Y应为正数。若X、Y为浮点数,由于计算机对浮点数的表示的不精确性,则结果将可能是不可意料的。fix(X./Y)为商数X./Y朝零方向取的整数部分。若X与Y为同符号的,则rem(X,Y)返回的结果与mod(X,Y)相同,不然,若X为正数,则rem(-X,Y) = mod(-X,Y) - Y。该命令返回的结果在区间[0,sign(X)*abs(Y)],若Y中有零分量,则相应地返回NaN。
例2-17
    >>X = [12 23 34 45];
    >>Y = [3 7 2 6];
    >>R = rem(X,Y)
    计算结果为:
    R =
        0     2     0     3
函数  ceil
功能  朝正无穷大方向取整
格式  B = floor(A)   % 对A的每一个元素朝正无穷大的方向取整数部分,返回与A同维的数组。对于复数参量A,则返回一复数,其分量的实数与虚数部分分别取原复数的、朝正无穷大方向的整数部分。
例2-18
>>A = [-1.9, -0.2, 3.1415926, 5.6, 7.0, 2.4+3.6i];
>>B = ceil(A)
计算结果为:
B =
    Columns 1 through 4
     -1.0000        0         4.0000         6.0000
    Columns 5 through 6
     7.0000    3.0000 + 4.0000i
函数  exp
功能  以e为底数的指数函数
格式  Y = exp(X)   %对参量X的每一分量,求以e为底数的指数函数Y。X中的分量可以为复数。对于复数分量如,z = x +i*y,则相应地计算:e^z = e^x*(cos(y) + i*sin(y))。
例2-19
>>A = [-1.9, -0.2, 3.1415926, 5.6, 7.0, 2.4+3.6i];
>>Y = exp(A)
计算结果为:
    Y =
       1.0e+003 *
       Columns 1 through 4
          0.0001        0.0008        0.0231         0.2704
       Columns 5 through 6
          1.0966   -0.0099 - 0.0049i
函数  expm
功能  求矩阵的以e为底数的指数函数
格式  Y = expm(X)   %计算以e为底数、x的每一个元素为指数的指数函数值。若矩阵x有小于等于零的特征值,则返回复数的结果。
说明  该函数为一内建函数,它有三种计算算法:
(1)使用文件expm1.m中的用比例法与二次幂算法得到的Pad近似值;
(2)使用Taylor级数近似展开式计算,这种计算在文件expm2.m中。但这种一般计算方法是不可取的,通常计算是缓慢且不精确的;
(3)在文件expm3.m中,先是将矩阵对角线化,再把函数计算出相应的的特征向量,最后转换过来。但当输入的矩阵没有与矩阵阶数相同的特征向量个数时,就会出现错误。
例2-20
>>A=hilb(4);
>>Y = expm(A)
计算结果为:
Y =
    3.2506    1.2068    0.8355    0.6417
    1.2068    1.7403    0.5417    0.4288
    0.8355    0.5417    1.4100    0.3318
    0.6417    0.4288    0.3318    1.2729
函数  log
功能  自然对数,即以e为底数的对数。
格式  Y = log(X)   %对参量X中的每一个元素计算自然对数。其中X中的元素可以是复数与负数,但由此可能得到意想不到的结果。若z = x + i*y,则log对复数的计算如下:log (z) = log (abs (z)) + i*atan2(y,x)
例2-21  下面的语句可以得到无理数π的近似值:
>>Pi = abs(log(-1))
计算结果为:
Pi =
    3.1416
函数  log10
功能  常用对数,即以10为底数的对数。
格式  Y = log10(X)   %计算X中的每一个元素的常用对数,若X中出现复数,则可能得到意想不到的结果。
例2-22
>>L1 = log10(realmax)  % 由此可得特殊变量realmax的近似值
>>L2 = log10(eps)  % 由此可得特殊变量eps的近似值
>>M = magic(4);
>>L3 = log10(M)
计算结果为:
L1 =
     308.2547
L2 =
     -15.6536
L3 =
     1.2041    0.3010    0.4771    1.1139
     0.6990    1.0414    1.0000    0.9031
     0.9542    0.8451    0.7782    1.0792
     0.6021    1.1461    1.1761         0
函数  sort
功能  把输入参量中的元素按从小到大的方向重新排列
格式  B = sort(A)   %沿着输入参量A的不同维的方向、从小到大重新排列A中的元素。A可以是字符串的、实数的、复数的单元数组。对于A中完全相同的元素,则按它们在A中的先后位置排列在一块;若A为复数的,则按元素幅值的从小到大排列,若有幅值相同的复数元素,则再按它们在区间[-π,π]的幅角从小到大排列;若A中有元素为NaN,则将它们排到最后。若A为向量,则返回从小到大的向量,若A为二维矩阵,则按列的方向进行排列;若A为多维数组,sort(A)把沿着第一非单元集的元素象向量一样进行处理。
B = sort(A,dim)        %沿着矩阵A(向量的、矩阵的或多维的)中指定维数dim方向重新排列A中的元素。
[B,INDEX] = sort(A,…)   %输出参量B的结果如同上面的情形,输出INDEX是一等于size(A)的数组,它的每一列是与A中列向量的元素相对应的置换向量。若A中有重复出现的相同的值,则返回保存原来相对位置的索引。
例2-23
>>A = [-1.9, -0.2, 3.1415926, 5.6, 7.0, 2.4+3.6i];
>>[B1,INDEX] = sort(A)
>>M = magic(4);
>>B2 = sort(M)
计算结果为:
    B1 =
        Columns 1 through 4
          -0.2000      -1.9000       3.1416         2.4000 + 3.6000i
        Columns 5 through 6
           5.6000      7.0000
    INDEX =
          2     1     3     6     4     5
    B2 =
         4     2     3     1
         5     7     6     8
         9     11    10    12
         16    14    15    13
函数  abs
功能  数值的绝对值与复数的幅值
格式  Y = abs(X)   %返回参量X的每一个分量的绝对值;若X为复数的,则返回每一分量的幅值:abs(X) = sqrt(real(X).^2+imag(X).^2)。
例2-24
>>A = [-1.9, -0.2, 3.1415926, 5.6, 7.0, 2.4+3.6i];
>>Y = abs(A)
计算结果为:
Y =
     1.9000    0.2000    3.1416    5.6000    7.0000    4.3267
函数  conj
功能  复数的共轭值
格式  ZC = conj(Z)   %返回参量Z的每一个分量的共轭复数:
conj(Z) = real(Z) - i*imag(Z)
函数  imag
功能  复数的虚数部分
格式  Y = imag(Z)    %返回输入参量Z的每一个分量的虚数部分。
例2-25
>>imag(2+3i)
计算结果为:
ans =
      3
函数  real
功能  复数的实数部分。
格式  Y = real(Z)   %返回输入参量Z的每一个分量的实数部分。
例2-26
>>real(2+3i)
计算结果为:
ans =
     2
函数  angle
功能  复数的相角
格式  P = angle(Z)   %返回输入参量Z的每一复数元素的、单位为弧度的相角,其值在区间[-π,π]上。
说明  angle(z) = imag (log(z)) = atan2 (imag(z),real(z))
例2-27
>>Z =[1-i, 2+i, 3-i, 4+i;
>>1+2i,2-2i,3+2i,4-2i;
>>1-3i,2+3i,3-3i,4+3i;
>>1+4i,2-4i,3+4i,4-4i];
>>P = angle(Z)
计算结果为:
       P =
           -0.7854    0.4636   -0.3218    0.2450
            1.1071   -0.7854    0.5880   -0.4636
           -1.2490    0.9828   -0.7854    0.6435
            1.3258   -1.1071    0.9273   -0.7854
函数  complex
功能  用实数与虚数部分创建复数
格式  c = complex(a,b)   %用两个实数a,b创建复数c=a+bi。输出参量c与a、b同型(同为向量、矩阵、或多维阵列)。该命令比下列形式的复数输入更有用:a + i*b 或a + j*b因为i和j可能被用做其他的变量(不等于sqrt(-1)),或者a和b不是双精度的。
c = complex(a)    %输入参量a作为输出复数c的实部,其虚部为0:c = a+0*i。
例2-28
>>a = uint8([1;2;3;4]);
>>b = uint8([4;3;2;1]);
>>c = complex(a,b)
计算结果为:
c =
   1.0000 + 4.0000i
   2.0000 + 3.0000i
   3.0000 + 2.0000i
   4.0000 + 1.0000i
函数  mod
功能  模数(带符号的除法余数)
用法  M = mod(X,Y)   %输入参量X、Y应为整数,此时返回余数X -Y.*floor(X./Y),若Y≠0,或者是X。若运算数x与y有相同的符号,则mod(X,Y)等于rem(X,Y)。总之,对于整数x,y,有:mod(-x,y) = rem(-x,y)+y。若输入为实数或复数,由于浮点数在计算机上的不精确表示,该操作将导致不可预测的结果。
例2-29
>>M1 = mod(13,5)
>>M2 = mod([1:5],3)
>>M3 = mod(magic(3),3)
计算结果为:
M1 =
     3
M2 =
     1     2     0     1     2
M3 =
     2     1     0
     0     2     1
     1     0     2
函数  nchoosek
功能  二项式系数或所有的组合数。该命令只有对n<15时有用。
函数  C = nchoosek(n,k)   %参量n,k为非负整数,返回n! / ( (n-k)! k!),即一次从n个物体中取出k个的组合数。
C = nchoosek(v,k)   %参量v为n维向量,返回一矩阵,其行向量的分量为一次性从v个物体中取k个物体的组合数。矩阵 C包含 =n! / ( (n-k)! k!)行与k列。
例2-30
>>C = nchoosek(2:2:10,4)
计算结果为:
       C =
           2     4     6     8
           2     4     6    10
           2     4     8    10
           2     6     8    10
           4     6     8    10
函数  rand
功能  生成元素均匀分布于(0,1)上的数值与阵列
用法  Y = rand(n)   %返回n*n阶的方阵Y,其元素均匀分布于区间(0,1)。若n不是一标量,在显示一出错信息。
Y = rand(m,n)、Y = rand([m n])          %返回阶数为m*n的,元素均匀分布于区间(0,1)上矩阵Y。
Y = rand(m,n,p,…)、Y = rand([m n p…])    %生成阶数m*n*p*…的,元素服从均匀分布的多维随机阵列Y。
Y = rand(size(A))   %生成一与阵列A同型的随机均匀阵列Y
rand            %该命令在每次单独使用时,都返回一随机数(服从均匀分布)。
s = rand('state')    %返回一有35元素的列向量s,其中包含均匀分布生成器的当前状态。该改变生成器的当前的状态,见表2-1。
表2-1
命   令 含   义 
Rand(’state’,s) 设置状态为s 
Rand(’state’,0) 设置生成器为初始状态 
Rand(’state’,k) 设置生成器第k个状态(k为整数) 
Rand(’state’,sum(100*clock)) 设置生成器在每次使用时的状态都不同(因为clock每次都不同) 
例:
>>R1 = rand(4,5)
>>a = 10; b = 50;
>>R2 = a + (b-a) * rand(5)  % 生成元素均匀分布于(10,50)上的矩阵
计算结果可能为:
R1 =
    0.6655    0.0563    0.2656    0.5371    0.6797
    0.3278    0.4402    0.9293    0.5457    0.6129
    0.6325    0.4412    0.9343    0.9394    0.3940
    0.5395    0.6501    0.5648    0.7084    0.2206
R2 =
   33.6835   19.8216   36.9436   49.6289   46.4679
   18.5164   34.2597   15.3663   31.0549   49.0377
   19.0026   37.1006   33.6046   39.5361   13.9336
   12.4641   12.9804   35.5420   23.2916   46.8304
   28.5238   48.7418   49.0843   13.0512   10.9265
函数  randn
功能  生成元素服从正态分布(N(0,1))的数值与阵列
格式  Y = randn(n)    %返回n*n阶的方阵Y,其元素服从正态分布N(0,1)。若n不是一标量,则显示一出错信息。
Y = randn(m,n)、Y = randn([m n])   %返回阶数为m*n的,元素均匀分布于区间(0,1)上矩阵Y。
Y = randn(m,n,p,…)、Y = randn([m n p…])    %生成阶数m*n*p*…的,元素服从正态分布的多维随机阵列Y。
Y = randn(size(A))    %生成一与阵列A同型的随机正态阵列Y
randn    %该命令在每次单独使用时,都返回一随机数(服从正态分布)。
s = randn('state')      %返回一有2元素的向量s,其中包含正态分布生成器的当前状态。该改变生成器的当前状态,见表2-2。
表2-2
命    令 含    义 
randn(’state’,s) 设置状态为s 
randn(’state’,0) 设置生成器为初始状态 
rand(’state’,k) 设置生成器第k个状态(k为整数) 
rand(’state’,sum(100*clock)) 设置生成器在每次使用时的状态都不同(因为clock每次都不同) 
例:
>>R1 = rand(4,5)
>>R2 = 0.6 + sqrt(0.1) * randn(5)
计算结果可能为:
R1 =
    0.2778    0.2681    0.5552    0.5167    0.8821
    0.2745    0.3710    0.1916    0.3385    0.5823
    0.9124    0.5129    0.4164    0.2993    0.0550
    0.4125    0.2697    0.1508    0.9370    0.5878
R2 =
    0.4632    0.9766    0.5410    0.6360    0.6931
    0.0733    0.9760    0.8295    0.9373    0.1775
    0.6396    0.5881    0.4140    0.6187    0.8259
    0.6910    0.7035    1.2904    0.5698    1.1134
    0.2375    0.6552    0.5569    0.3368    0.3812
2.2  插值、拟合与查表
插值法是实用的数值方法,是函数逼近的重要方法。在生产和科学实验中,自变量x与因变量y的函数y = f(x)的关系式有时不能直接写出表达式,而只能得到函数在若干个点的函数值或导数值。当要求知道观测点之外的函数值时,需要估计函数值在该点的值。
如何根据观测点的值,构造一个比较简单的函数y=φ(x),使函数在观测点的值等于已知的数值或导数值。用简单函数y=φ(x)在点x处的值来估计未知函数y=f(x)在x点的值。寻找这样的函数φ(x),办法是很多的。φ(x)可以是一个代数多项式,或是三角多项式,也可以是有理分式;φ(x)可以是任意光滑(任意阶导数连续)的函数或是分段函数。函数类的不同,自然地有不同的逼近效果。在许多应用中,通常要用一个解析函数(一、二元函数)来描述观测数据。
根据测量数据的类型:
1.测量值是准确的,没有误差。
2.测量值与真实值有误差。
这时对应地有两种处理观测数据方法:
1.插值或曲线拟合。
2.回归分析(假定数据测量是精确时,一般用插值法,否则用曲线拟合)。
MATLAB中提供了众多的数据处理命令。有插值命令,有拟合命令,有查表命令。
2.2.1  插值命令
命令1  interp1
功能  一维数据插值(表格查找)。该命令对数据点之间计算内插值。它找出一元函数f(x)在中间点的数值。其中函数f(x)由所给数据决定。各个参量之间的关系示意图为图2-14。
 
图2-14  数据点与插值点关系示意图
格式  yi = interp1(x,Y,xi)        %返回插值向量yi,每一元素对应于参量xi,同时由向量x与Y的内插值决定。参量x指定数据Y的点。若Y为一矩阵,则按Y的每列计算。yi是阶数为length(xi)*size(Y,2)的输出矩阵。
yi = interp1(Y,xi)          %假定x=1:N,其中N为向量Y的长度,或者为矩阵Y的行数。
yi = interp1(x,Y,xi,method)   %用指定的算法计算插值:
’nearest’:最近邻点插值,直接完成计算;
’linear’:线性插值(缺省方式),直接完成计算;
’spline’:三次样条函数插值。对于该方法,命令interp1调用函数spline、ppval、mkpp、umkpp。这些命令生成一系列用于分段多项式操作的函数。命令spline用它们执行三次样条函数插值;
’pchip’:分段三次Hermite插值。对于该方法,命令interp1调用函数pchip,用于对向量x与y执行分段三次内插值。该方法保留单调性与数据的外形;
’cubic’:与’pchip’操作相同;
’v5cubic’:在MATLAB 5.0中的三次插值。
对于超出x范围的xi的分量,使用方法’nearest’、’linear’、’v5cubic’的插值算法,相应地将返回NaN。对其他的方法,interp1将对超出的分量执行外插值算法。
yi = interp1(x,Y,xi,method,'extrap')     %对于超出x范围的xi中的分量将执行特殊的外插值法extrap。
yi = interp1(x,Y,xi,method,extrapval)  %确定超出x范围的xi中的分量的外插值extrapval,其值通常取NaN或0。
例2-31
>>x = 0:10; y = x.*sin(x);
>>xx = 0:.25:10; yy = interp1(x,y,xx);
>>plot(x,y,'kd',xx,yy)
插值图形为图2-15。
例2-32
>> year = 1900:10:2010;
>> product = [75.995  91.972  105.711  123.203  131.669  150.697  179.323  203.212  226.505  249.633  256.344  267.893 ];
>>p1995 = interp1(year,product,1995)
>>x = 1900:1:2010;
>>y = interp1(year,product,x,'pchip');
>>plot(year,product,'o',x,y)
插值结果为:
p1995 =
        252.9885
插值图形为图2-16。
    
图2-15  一元函数插值图形                图2-16  离散数据的一维插值图
命令2  interp2
功能  二维数据内插值(表格查找)
格式  ZI = interp2(X,Y,Z,XI,YI)   %返回矩阵ZI,其元素包含对应于参量XI与YI(可以是向量、或同型矩阵)的元素,即Zi(i,j)←[Xi(i,j),yi(i,j)]。用户可以输入行向量和列向量Xi与Yi,此时,输出向量Zi与矩阵meshgrid(xi,yi)是同型的。同时取决于由输入矩阵X、Y与Z确定的二维函数Z=f(X,Y)。参量X与Y必须是单调的,且相同的划分格式,就像由命令meshgrid生成的一样。若Xi与Yi中有在X与Y范围之外的点,则相应地返回nan(Not a Number)。
ZI = interp2(Z,XI,YI)      %缺省地,X=1:n、Y=1:m,其中[m,n]=size(Z)。再按第一种情形进行计算。
ZI = interp2(Z,n)         %作n次递归计算,在Z的每两个元素之间插入它们的二维插值,这样,Z的阶数将不断增加。interp2(Z)等价于interp2(z,1)。
ZI = interp2(X,Y,Z,XI,YI,method)   %用指定的算法method计算二维插值:
’linear’:双线性插值算法(缺省算法);
’nearest’:最临近插值;
’spline’:三次样条插值;
’cubic’:双三次插值。
例2-33:
>>[X,Y] = meshgrid(-3:.25:3);
>>Z = peaks(X,Y);
>>[XI,YI] = meshgrid(-3:.125:3);
>>ZZ = interp2(X,Y,Z,XI,YI);
>>surfl(X,Y,Z);hold on;
>>surfl(XI,YI,ZZ+15)
>>axis([-3 3 -3 3 -5 20]);shading flat
>>hold off
插值图形为图2-17。
例2-34
>>years = 1950:10:1990;
>>service = 10:10:30;
>>wage = [150.697 199.592 187.625
       179.323 195.072 250.287
       203.212 179.092 322.767
       226.505 153.706 426.730
       249.633 120.281 598.243];
>>w = interp2(service,years,wage,15,1975)
插值结果为:
w =
190.6288
命令3  interp3
功能  三维数据插值(查表)
格式  VI = interp3(X,Y,Z,V,XI,YI,ZI)   %找出由参量X,Y,Z决定的三元函数V=V(X,Y,Z)在点(XI,YI,ZI)的值。参量XI,YI,ZI是同型阵列或向量。若向量参量XI,YI,ZI是不同长度,不同方向(行或列)的向量,这时输出参量VI与Y1,Y2,Y3为同型矩阵。其中Y1,Y2,Y3为用命令meshgrid(XI,YI,ZI)生成的同型阵列。若插值点(XI,YI,ZI)中有位于点(X,Y,Z)之外的点,则相应地返回特殊变量值NaN。
VI = interp3(V,XI,YI,ZI)     %缺省地,X=1:N,Y=1:M,Z=1:P,其中,[M,N,P]=size(V),再按上面的情形计算。
VI = interp3(V,n)             %作n次递归计算,在V的每两个元素之间插入它们的三维插值。这样,V的阶数将不断增加。interp3(V)等价于interp3(V,1)。
VI = interp3(…,method)         %用指定的算法method作插值计算:
                    ‘linear’:线性插值(缺省算法);
                    ‘cubic’:三次插值;
                    ‘spline’:三次样条插值;
                    ‘nearest’:最邻近插值。
说明  在所有的算法中,都要求X,Y,Z是单调且有相同的格点形式。当X,Y,Z是等距且单调时,用算法’*linear’,’*cubic’,’*nearest’,可得到快速插值。
例2-35
>>[x,y,z,v] = flow(20);
>>[xx,yy,zz] = meshgrid(.1:.25:10, -3:.25:3, -3:.25:3);
>>vv = interp3(x,y,z,v,xx,yy,zz);
>>slice(xx,yy,zz,vv,[6 9.5],[1 2],[-2 .2]); shading interp;colormap cool
插值图形为图2-18。
图2-18  三维插值图
命令4  interpft
功能  用快速Fourier算法作一维插值
格式  y = interpft(x,n)   %返回包含周期函数x在重采样的n个等距的点的插值y。若length(x)=m,且x有采样间隔dx,则新的y的采样间隔dy=dx*m/n。注意的是必须n≥m。若x为一矩阵,则按x的列进行计算。返回的矩阵y有与x相同的列数,但有n行。
y = interpft(x,n,dim)     %沿着指定的方向dim进行计算
命令5  griddata
功能  数据格点
格式  ZI = griddata(x,y,z,XI,YI)   %用二元函数z=f(x,y)的曲面拟合有不规则的数据向量x,y,z。griddata将返回曲面z在点(XI,YI)处的插值。曲面总是经过这些数据点(x,y,z)的。输入参量(XI,YI)通常是规则的格点(像用命令meshgrid生成的一样)。XI可以是一行向量,这时XI指定一有常数列向量的矩阵。类似地,YI可以是一列向量,它指定一有常数行向量的矩阵。
[XI,YI,ZI] = griddata(x,y,z,xi,yi)   %返回的矩阵ZI含义同上,同时,返回的矩阵XI,YI是由行向量xi与列向量yi用命令meshgrid生成的。
[…] = griddata(…,method)         %用指定的算法method计算:
                  ‘linear’:基于三角形的线性插值(缺省算法);
                  ‘cubic’: 基于三角形的三次插值;
                  ‘nearest’:最邻近插值法;
                  ‘v4’:MATLAB 4中的griddata算法。
命令6  spline
功能  三次样条数据插值
格式  yy = spline(x,y,xx)   %对于给定的离散的测量数据x,y(称为断点),要寻找一个三项多项式 ,以逼近每对数据(x,y)点间的曲线。过两点 和 只能确定一条直线,而通过一点的三次多项式曲线有无穷多条。为使通过中间断点的三次多项式曲线具有唯一性,要增加两个条件(因为三次多项式有4个系数):
1.三次多项式在点 处有:  ;
2.三次多项式在点 处有: ;
3.p(x)在点 处的斜率是连续的(为了使三次多项式具有良好的解析性,加上的条件);
4.p(x)在点 处的曲率是连续的;
对于第一个和最后一个多项式,人为地规定如下条件:
                  ①.
                  ②.
                上述两个条件称为非结点(not-a-knot)条件。综合上述内容,可知对数据拟合的三次样条函数p(x)是一个分段的三次多项式:
                 ,其中每段 都是三次多项式。
                该命令用三次样条插值计算出由向量x与y确定的一元函数y=f(x)在点xx处的值。若参量y是一矩阵,则以y的每一列和x配对,再分别计算由它们确定的函数在点xx处的值。则yy是一阶数为length(xx)*size(y,2)的矩阵。
pp = spline(x,y)   %返回由向量x与y确定的分段样条多项式的系数矩阵pp,它可用于命令ppval、unmkpp的计算。
例2-36
对离散地分布在y=exp(x)sin(x)函数曲线上的数据点进行样条插值计算:
    >>x = [0 2 4 5 8 12 12.8 17.2 19.9 20];  y = exp(x).*sin(x);
    >>xx = 0:.25:20;
    >>yy = spline(x,y,xx);
    >>plot(x,y,'o',xx,yy)
插值图形结果为图2-19。
图2-19  三次样条插值
命令7  interpn
功能  n维数据插值(查表)
格式  VI = interpn(X1,X2,,…,Xn,V,Y1,Y2,…,Yn)   %返回由参量X1,X2,…,Xn,V确定的n元函数V=V(X1,X2,…,Xn)在点(Y1,Y2,…,Yn)处的插值。参量Y1,Y2,…,Yn是同型的矩阵或向量。若Y1,Y2,…,Yn是向量,则可以是不同长度,不同方向(行或列)的向量。它们将通过命令ndgrid生成同型的矩阵,再作计算。若点(Y1,Y2,…,Yn)中有位于点(X1,X2,…,Xn)之外的点,则相应地返回特殊变量NaN。
VI = interpn(V,Y1,Y2,…,Yn)   %缺省地,X1=1:size(V,1),X2=1:size(V,2),…,Xn=1:size(V,n),再按上面的情形计算。
VI = interpn(V,ntimes)   %作ntimes次递归计算,在V的每两个元素之间插入它们的n维插值。这样,V的阶数将不断增加。interpn(V)等价于interpn(V,1)。
VI = interpn(…,method)   %用指定的算法method计算:
                   ‘linear’:线性插值(缺省算法);
                   ‘cubic’:三次插值;
                   ‘spline’:三次样条插值法;
                   ‘nearest’:最邻近插值算法。
命令8  meshgrid
功能  生成用于画三维图形的矩阵数据。
格式  [X,Y] = meshgrid(x,y)   将由向量x,y(可以是不同方向的)指定的区域[min(x),max(x),min(y),max(y)]用直线x=x(i),y=y(j)(i=1,2,…,length(x) ,j=1,2,…,length(y))进行划分。这样,得到了length(x)*length(y)个点,这些点的横坐标用矩阵X表示,X的每个行向量与向量x相同;这些点的纵坐标用矩阵Y表示,Y的每个列向量与向量y相同。其中X,Y可用于计算二元函数z=f(x,y)与三维图形中xy平面矩形定义域的划分或曲面作图。
[X,Y] = meshgrid(x)        %等价于[X,Y]=meshgrid(x,x)。
[X,Y,Z] = meshgrid(x,y,z)   %生成三维阵列X,Y,Z,用于计算三元函数v=f(x,y,z)或三维容积图。
例2-37
    [X,Y] = meshgrid(1:3,10:14)
计算结果为:
    X =
       1     2     3
       1     2     3
       1     2     3
       1     2     3
       1     2     3
   Y =
      10    10    10
      11    11    11
      12    12    12
      13    13    13
      14    14    14
命令9  ndgrid
功能  生成用于多维函数计算或多维插值用的阵列
格式  [X1,X2,…,Xn] = ndgrid(x1,x2,…,xn)  %把通过向量x1,x2,x3…,xn指定的区域转换为数组x1,x2,x3,…,xn。这样,得到了   length(x1)* length(x2)*…*length(xn)个点,这些点的第一维坐标用矩阵X1表示,X1的每个第一维向量与向量x1相同;这些点的第二维坐标用矩阵X2表示,X2的每个第二维向量与向量x2相同;如此等等。其中X1,X2,…,Xn可用于计算多元函数y=f(x1,x2,…,xn)以及多维插值命令用到的阵列。
[X1,X2,…,Xn] = ndgrid(x)   %等价于[X1,X2,…,Xn] = ndgrid(x,x,…,x)
2.2.2  查表命令
命令1  table1
功能  一维查表
格式  Y = table1(TAB,X0)   %返回用表格矩阵TAB中的行线性插值元素,对X0(TAB的第一列查找X0)进行线性插值得到的结果Y。矩阵TAB是第一列包含关键值,而其他列包含数据的矩阵。X0中的每一元素将相应地返回一线性插值行向量。矩阵TAB的第一列必须是单调的。
例2-38 
>>tab = [(1:4)' hilb(4)]
>>y = table1(tab,[1 2.3 3.6 4])
查表结果为:
tab =
      1.0000    1.0000    0.5000    0.3333    0.2500
      2.0000    0.5000    0.3333    0.2500    0.2000
      3.0000    0.3333    0.2500    0.2000    0.1667
      4.0000    0.2500    0.2000    0.1667    0.1429
Warning: TABLE1 is obsolete and will be removed in future versions.  Use INTERP1 or INTERP1Q instead.
> In D:MATLABR12 oolboxmatlabpolyfun able1.m at line 31
y =
    1.0000    0.5000    0.3333    0.2500
    0.4500    0.3083    0.2350    0.1900
    0.2833    0.2200    0.1800    0.1524
    0.2500    0.2000    0.1667    0.1429
由上面结果可知,table1是一将要废弃的命令。
命令2  table2
功能  二维查表
格式  Z = table1(TAB,X0,Y0)   %返回用表格矩阵TAB中的行与列交叉线性线性插值元素,对X0(TAB的第一列查找X0)进行线性插值,对Y0(TAB的第一行查找Y0)进行线性插值,对上述两个数值进行交叉线性插值,得到的结果为Z。矩阵TAB是第一列与第一行列都包含关键值,而其他的元素包含数据的矩阵。TAB(1,1)的关键值将被忽略。[X0,Y0]中的每点将相应地返回一线性插值。矩阵TAB的第一行与第一列必须是单调的。
例2-39
>>tab = [NaN 1:4; (1:4)' magic(4)]
>>y = table2(tab,[2 3 3.7],[1.3 2.3 4])
查表的结果为:
tab =
      NaN     1     2     3     4
        1    16     2     3    13
        2     5    11    10     8
        3     9     7     6    12
        4     4    14    15     1
 Warning: TABLE2 is obsolete and will be removed in future versions.  Use INTERP2 instead.
 > In D:MATLABR12 oolboxmatlabpolyfun able2.m at line 24
 Warning: TABLE1 is obsolete and will be removed in future versions.  Use INTERP1 or INTERP1Q instead.
 > In D:MATLABR12 oolboxmatlabpolyfun able1.m at line 31
 In D:MATLABR12 oolboxmatlabpolyfun able2.m at line 29
 Warning: TABLE1 is obsolete and will be removed in future versions.  Use INTERP1 or INTERP1Q instead.
 > In D:MATLABR12 oolboxmatlabpolyfun able1.m at line 31
 In D:MATLABR12 oolboxmatlabpolyfun able2.m at line 31
 y =
    6.8000   10.7000    8.0000
    8.4000    6.7000   12.0000
    7.4200   12.0200    4.3000
由上面的结果可知,table2是将要废弃的命令。
2.3  数值积分
2.3.1  一元函数的数值积分
函数1  quad、quadl、quad8
功能  数值定积分,自适应Simpleson积分法。
格式  q = quad(fun,a,b)    %近似地从a到b计算函数fun的数值积分,误差为10-6。若给fun输入向量x,应返回向量y,即fun是一单值函数。
q = quad(fun,a,b,tol)   %用指定的绝对误差tol代替缺省误差。tol越大,函数计算的次数越少,速度越快,但结果精度变小。
q = quad(fun,a,b,tol,trace,p1,p2,…)    %将可选参数p1,p2,…等传递给函数fun(x,p1,p2,…),再作数值积分。若tol=[]或trace=[],则用缺省值进行计算。
[q,n] = quad(fun,a,b,…)  %同时返回函数计算的次数n
… = quadl(fun,a,b,…)   %用高精度进行计算,效率可能比quad更好。
… = quad8(fun,a,b,…)   %该命令是将废弃的命令,用quadl代替。
例2-40
>>fun = inline(‘’3*x.^2./(x.^3-2*x.^2+3)’);
>>Q1 = quad(fun,0,2)
>>Q2 = quadl(fun,0,2)
计算结果为:
Q1 =
    3.7224
Q2 =
    3.7224
函数2  trapz
功能  梯形法数值积分
格式  T = trapz(Y)       %用等距梯形法近似计算Y的积分。若Y是一向量,则trapz(Y)为Y的积分;若Y是一矩阵,则trapz(Y)为Y的每一列的积分;若Y是一多维阵列,则trapz(Y)沿着Y的第一个非单元集的方向进行计算。
T = trapz(X,Y)     %用梯形法计算Y在X点上的积分。若X为一列向量,Y为矩阵,且size(Y,1) = length(X),则trapz(X,Y)通过Y的第一个非单元集方向进行计算。
T = trapz(…,dim)   %沿着dim指定的方向对Y进行积分。若参量中包含X,则应有length(X)=size(Y,dim)。
例2-41
>>X = -1:.1:1;
>>Y = 1./(1+25*X.^2);
>>T = trapz(X,Y)
计算结果为:
T =
    0.5492
函数3  rat,rats
功能  有理分式近似。虽然所有的浮点数值都是有理数,有时用简单的有理数字(分子与分母都是较小的整数)近似地表示它们是有必要的。函数rat将试图做到这一点。对于有连续出现的小数的数值,将会用有理式近似表示它们。函数rats调用函数rat,且返回字符串。
格式  [N,D] = rat(X)       %对于缺省的误差1.e-6*norm(X(:),1),返回阵列N与D,使N./D近似为X。
[N,D] = rat(X,tol)    %在指定的误差tol范围内,返回阵列N与D,使N./D近似为X。
rat(X)、rat(X…)      %在没有输出参量时,简单地显示x的连续分数。
S = rats(X,strlen)     %返回一包含简单形式的、X中每一元素的有理近似字符串S,若对于分配的空间中不能显示某一元素,则用星号表示。该元素与X中其他元素进行比较而言较小,但并非是可以忽略。参量strlen为函数rats中返回的字符串元素的长度。缺省值为strlen=13,这允许在78个空格中有6个元素。
S = rats(X)         %返回与用MATLAB命令format rat显示 X相同的结果给S。
例2-42
>>s = 1-1/2+1/3-1/4+1/5-1/6+1/7
>>format rat
>>S1 = rats(s)
>>S2 = rat(s)
>>[n,d] = rat(s)
>>PI1 = rats(pi)
>>PI2 = rat(pi)
计算结果为:
s =
    0.7595
S1 =
    319/420  
S2 =
    1 + 1/(-4 + 1/(-6 + 1/(-3 + 1/(-5))))
n =
    319     
d =
    420     
PI1 =
    355/113  
PI2 =
3 + 1/(7 + 1/(16))
2.3.2  二元函数重积分的数值计算
函数1  dblquad
功能  矩形区域上的二重积分的数值计算
格式  q = dblquad(fun,xmin,xmax,ymin,ymax)    %调用函数quad在区域[xmin,xmax, ymin,ymax]上计算二元函数z=f(x,y)的二重积分。输入向量x,标量y,则f(x,y)必须返回一用于积分的向量。
q = dblquad(fun,xmin,xmax,ymin,ymax,tol)   %用指定的精度tol代替缺省精度10-6,再进行计算。
q = dblquad(fun,xmin,xmax,ymin,ymax,tol,method)   %用指定的算法method代替缺省算法quad。method的取值有@quadl或用户指定的、与命令quad与quadl有相同调用次序的函数句柄。
q = dblquad(fun,xmin,xmax,ymin,ymax,tol,method,p1,p2,…)   %将可选参数p1,p2,..等传递给函数fun(x,y,p1,p2,…)。若tol=[],method=[],则使用缺省精度和算法quad。
例2-43
>>fun = inline(’y./sin(x)+x.*exp(y)’);
>>Q = dblquad(fun,1,3,5,7)
计算结果为:
Q =
  3.8319e+003
函数2  quad2dggen
功能  任意区域上二元函数的数值积分
格式  q = quad2dggen(fun,xlower,xupper,ymin,ymax)  %在由[xlower,xupper, ymin,ymax]指定的区域上计算二元函数z=f(x,y)的二重积分。
q = dblquad(fun,xlower,xupper,ymin,ymax,tol)  %用指定的精度tol代替缺省精度10-6,再进行计算。
q = dblquad(fun,xmin,xmax,ymin,ymax,tol,method)   %用指定的算法method代替缺省算法。method的取值有缺省算法或用户指定的、与缺省命令有相同调用次序的函数句柄。
q=dblquad(fun,xlower,xupper,ymin,ymax,tol,method,p1,p2,…)  %将可选参数p1,p2,..等传递给函数fun(x,y,p1,p2,…)。若tol=[],method=[],则使用缺省精度和算法。
例2-44
计算单位圆域上的积分:
先把二重积分转化为二次积分的形式:
f = inline(’exp(-x.^2/2).*sin(x.^2+y)’,’x’,’y’);
xlower = inline(’-sqrt(1-y.^2)’,’y’); xupper = inline(’sqrt(1-y.^2)’,’y’);
Q = quad2dggen(fun,xlower,xupper,-1,1,1e-4)
计算结果为:
  Q =
      0.5368603818
2.4  常微分方程数值解
函数  ode45、ode23、ode113、ode15s、ode23s、ode23t、ode23tb
功能  常微分方程(ODE)组初值问题的数值解
参数说明:
solver为命令ode45、ode23,ode113,ode15s,ode23s,ode23t,ode23tb之一。
Odefun 为显式常微分方程y’=f(t,y),或为包含一混合矩阵的方程M(t,y)*y’=f(t,y)。命令ode23只能求解常数混合矩阵的问题;命令ode23t与ode15s可以求解奇异矩阵的问题。
Tspan 积分区间(即求解区间)的向量tspan=[t0,tf]。要获得问题在其他指定时间点t0,t1,t2,…上的解,则令tspan=[t0,t1,t2,…,tf](要求是单调的)。
Y0 包含初始条件的向量。
Options 用命令odeset设置的可选积分参数。
P1,p2,… 传递给函数odefun的可选参数。
格式  [T,Y] = solver(odefun,tspan,y0)   %在区间tspan=[t0,tf]上,从t0到tf,用初始条件y0求解显式微分方程y’=f(t,y)。对于标量t与列向量y,函数f=odefun(t,y)必须返回一f(t,y)的列向量f。解矩阵Y中的每一行对应于返回的时间列向量T中的一个时间点。要获得问题在其他指定时间点t0,t1,t2,…上的解,则令tspan=[t0,t1,t2,…,tf](要求是单调的)。
[T,Y] = solver(odefun,tspan,y0,options)   %用参数options(用命令odeset生成)设置的属性(代替了缺省的积分参数),再进行操作。常用的属性包括相对误差值RelTol(缺省值为1e-3)与绝对误差向量AbsTol(缺省值为每一元素为1e-6)。
[T,Y] =solver(odefun,tspan,y0,options,p1,p2…) 将参数p1,p2,p3,..等传递给函数odefun,再进行计算。若没有参数设置,则令options=[]。
1.求解具体ODE的基本过程:
(1)根据问题所属学科中的规律、定律、公式,用微分方程与初始条件进行描述。
F(y,y’,y’’,…,y(n),t) = 0
y(0)=y0,y’(0)=y1,…,y(n-1)(0)=yn-1
而y=[y;y(1);y(2);…,y(m-1)],n与m可以不等
(2)运用数学中的变量替换:yn=y(n-1),yn-1=y(n-2),…,y2=y1=y,把高阶(大于2阶)的方程(组)写成一阶微分方程组: ,
(3)根据(1)与(2)的结果,编写能计算导数的M-函数文件odefile。
(4)将文件odefile与初始条件传递给求解器Solver中的一个,运行后就可得到ODE的、在指定时间区间上的解列向量y(其中包含y及不同阶的导数)。
2.求解器Solver与方程组的关系表见表2-3。
表2-3
函数指令 含  义 函  数 含    义 
求解器 Solver ode23 普通2-3阶法解ODE odefile 包含ODE的文件 
 ode23s 低阶法解刚性ODE 选项 odeset 创建、更改Solver选项 
 ode23t 解适度刚性ODE  odeget 读取Solver的设置值 
 ode23tb 低阶法解刚性ODE 输出 odeplot ODE的时间序列图 
 ode45 普通4-5阶法解ODE  odephas2 ODE的二维相平面图 
 ode15s 变阶法解刚性ODE  odephas3 ODE的三维相平面图 
 ode113 普通变阶法解ODE  odeprint 在命令窗口输出结果 
3.因为没有一种算法可以有效地解决所有的ODE问题,为此,MATLAB提供了多种求解器Solver,对于不同的ODE问题,采用不同的Solver。
表2-4  不同求解器Solver的特点
求解器Solver ODE类型 特点 说明 
ode45 非刚性 一步算法;4,5阶Runge-Kutta方程;累计截断误差达(△x)3 大部分场合的首选算法 
ode23 非刚性 一步算法;2,3阶Runge-Kutta方程;累计截断误差达(△x)3 使用于精度较低的情形 
ode113 非刚性 多步法;Adams算法;高低精度均可到10-3~10-6 计算时间比ode45短 
ode23t 适度刚性 采用梯形算法 适度刚性情形 
ode15s 刚性 多步法;Gear’s反向数值微分;精度中等 若ode45失效时,可尝试使用 
ode23s 刚性 一步法;2阶Rosebrock算法;低精度 当精度较低时,计算时间比ode15s短 
ode23tb 刚性 梯形算法;低精度 当精度较低时,计算时间比ode15s短 
4.在计算过程中,用户可以对求解指令solver中的具体执行参数进行设置(如绝对误差、相对误差、步长等)。
表2-5  Solver中options的属性
属性名 取值 含义 
AbsTol 有效值:正实数或向量 缺省值:1e-6 绝对误差对应于解向量中的所有元素;向量则分别对应于解向量中的每一分量 
RelTol 有效值:正实数 缺省值:1e-3 相对误差对应于解向量中的所有元素。在每步(第k步)计算过程中,误差估计为: e(k)<=max(RelTol*abs(y(k)),AbsTol(k)) 
NormControl 有效值:on、off 缺省值:off 为‘on’时,控制解向量范数的相对误差,使每步计算中,满足:norm(e)<=max(RelTol*norm(y),AbsTol) 
Events 有效值:on、off 为‘on’时,返回相应的事件记录 
OutputFcn 有效值:odeplot、odephas2、odephas3、odeprint 缺省值:odeplot 若无输出参量,则solver将执行下面操作之一: 画出解向量中各元素随时间的变化; 画出解向量中前两个分量构成的相平面图; 画出解向量中前三个分量构成的三维相空间图; 随计算过程,显示解向量 
OutputSel 有效值:正整数向量 缺省值:[] 若不使用缺省设置,则OutputFcn所表现的是那些正整数指定的解向量中的分量的曲线或数据。若为缺省值时,则缺省地按上面情形进行操作 
Refine 有效值:正整数k>1 缺省值:k = 1 若k>1,则增加每个积分步中的数据点记录,使解曲线更加的光滑 
Jacobian 有效值:on、off 缺省值:off 若为‘on’时,返回相应的ode函数的Jacobi矩阵 
Jpattern 有效值:on、off 缺省值:off 为‘on’时,返回相应的ode函数的稀疏Jacobi矩阵 
Mass 有效值:none、M、 M(t)、M(t,y) 缺省值:none M:不随时间变化的常数矩阵 M(t):随时间变化的矩阵 M(t,y):随时间、地点变化的矩阵 
MaxStep 有效值:正实数 缺省值:tspans/10 最大积分步长 
例2-45  求解描述振荡器的经典的Ver der Pol微分方程
y(0)=1,y’(0)=0
令x1=y,x2=dy/dx,则
dx1/dt = x2
dx2/dt = μ(1-x2)-x1
编写函数文件verderpol.m:
function xprime = verderpol(t,x)
global MU
xprime = [x(2);MU*(1-x(1)^2)*x(2)-x(1)];
再在命令窗口中执行:
>>global MU
>>MU = 7;
>>Y0=[1;0]
>>[t,x] = ode45(‘verderpol’,0,40,Y0);
>>x1=x(:,1);x2=x(:,2);
>>plot(t,x1,t,x2)
图形结果为图2-20。
图2-20  Ver der Pol微分方程图
2.5  偏微分方程的数值解
MATLAB提供了一个专门用于求解偏微分方程的工具箱—PDE Toolbox (Paticial Difference Equation )。在本章,我们仅提供一些最简单、经典的偏微分方程,如:椭圆型、双曲型、抛物型等少数的偏微分方程,并给出求解方法。用户可以从中了解其解题基本方法,从而解决相类似的问题。
Matlab能解决的偏微分类型
    其中u = u(x,y),
 ,
c = c(x,y) , ,
2.5.1  单的Poission方程
Poission方程是特殊的椭圆型方程: ,
即 c = 1,a = 0,f = -1
Poission的解析解为: 。在下面计算中,用求得的数值解与精确解进行比较,看误差如何。
方程求解
(1)问题输入:
  c = 1;a = 0;f = 1;%方程的输入。给c,a,f赋值即可
  g = 'circleg'       %  区域G,内部已经定义为 circleg
   b = 'circleb1'      %  u在区域G的边界上的条件,内部已经定义好
(2)对单位圆进行网格化,对求解区域G作剖分,且作的是三角分划:
  [p,e,t] = initmesh(g,'hmax',1)
(3)迭代求解:
  error = [];err = 1;
  while err > 0.001,
  [p,e,t]=refinemesh('circleg',p,e,t);
  u=assempde('circleb1',p,e,t,1,0,1);
  exact=-(p(1,:)^2+p(2,:)^2-1)/4;
  err=norm(u-exact',inf);
  error=[error,err];
  end
(4)误差显示与区域G内的剖分点数:
  Error: 1.292265e-002. Number of nodes: 25
  Error: 4.079923e-003. Number of nodes: 81
  Error: 1.221020e-003. Number of nodes: 289
  Error: 3.547924e-004. Number of nodes: 1089
(5)结果显示:
  subplot(2,2,1),pdemesh(p,e,t)%结果显示
  title('数值解')
  subplot(2,2,2),pdesurf(p,t,u)%精确解显示
  title('精确解')
  subplot(2,2,3),pdesurf(p,t,u-exact')%与精确解的误差
  title('计算误差')
图2-21  Poission方程图
2.5.2  双曲型偏微分方程
1.Matlab能求解的类型:
其中 , , , , ,
2.形传递问题: ,
即:c =1; a = 0; f = 0; d = 1
3.方程求解
(1)问题的输入:
  c = 1; a = 0; f = 0; d = 1;   % 输入方程的系数
  g = 'squareg'    % 输入方形区域G,内部已经定义好
  b = 'sqareb3'    % 输入边界条件,即初始条件
(2)对单位矩形G进行网格化:
  [p,e,t] = initmesh('squareg')
(3)定解条件和求解时间点:
  x = p(1,:)'; y = p(2,:)';
  u0 = atan(cos(pi/2*x));
  ut0 = 3*sin(pi*x).*exp(sin(pi/2*y));
  n = 31;
  tlist = linspace(0,5,n);
(4)求解:
  uu = hyperbolic(uo, ut0,tlist,b,p,e,t,c,a,f,d);
结果显示:计算过程中的时间点和信息
Time:0.166667
Time:0.333333
Time:4.33333
……
Time:4.66667
Time:4.83333
Time:5
428 successful steps
62 failed attempts
982 function evaluations
1 partial derivatives
142 LU decompositions
981 solutions of linear systems
(5)动画显示:
delta=-1:0.1:1;
[uxy,tn,a2,a3]=tri2grid(p,t,uu(:,1),delta,delta);
gp=[tn;a2;a3];
umax=max(max(uu));
umin=min(min(uu));
newplot;M=moviein(n);
for i=1:n,
    pdeplot(p,e,t,'xydata',uu(:,i),'zdata',uu(:,i),…
    'mesh','off','xygrid','on','gridparam',gp,…
    'colorbar','off','zstyle','continuous');
    axis([-1 1 -1 1 umin umax]);
    caxis([umin umax]);
    M(:,i)=getframe;
end
movie(M,5)
图2-22是动画过程中的一个状态。
图2-22  波动方程动画中的一个状态
2.5.3  抛物型偏微分方程
1.Matlab能求解的类型:
其中 , , , , ,
2.热传导方程: ,
即:c = 1; a = 0; f = 1; d = 1;
3.问题计算
(1)问题的输入:
c = 1; a = 0; f = 1; d = 1;  %  输入方程的系数
g = 'squareg';  %  输入方形区域G
b = 'squareb1';  %  输入边界条件
(2)对单位矩形的网格化:
[p,e,t] = initmesh(g);
(3)定解条件和求解的时间点:
u0 = zeros(size(p, 2), 1);
ix = find(sqrt(p(1, :).^2+p(2, :).^2) < 0.4);
u0(ix) = ones(size(ix));
nframes = 20;
tlist=linspace(0,0.1,nframes)  % 在时间[0, 0.1]内20个点上计算,生成20帧
(4)求解方程:
u1 = parabolic(u0, tlist, b, p, e, t, c, a, f, d)
计算结果如下:
Time: 0.00526316
Time: 0.0105263
……
Time: 0.0947368
Time: 0.1
75 successful steps
1 failed attempts
154 function evaluations
1 partial derivatives
17 LU decompositions
153 solutions of linear systems
(5)动画显示:
x = linspace(-1,1,31); y = x;
newplot;
Mv = moviein(nframes);
umax=max(max(u1));
umin=min(min(u1));
for j=1:nframes
   u=tri2grid(p,t,u1(:,j),tn,a2,a3);i=find(isnan(u));u(i)=zeros(size(i));…
   surf(x,y,u);caxis([umin umax]);colormap(cool),…
   axis([-1 1 -1 1 0 1]);…
   Mv(:,j) = getframe;…
end
movie(Mv,10)
图2-23是动画过程中的瞬间状态。
图2-23  热传导方程动画瞬间状态图
  评论这张
 
阅读(480)| 评论(0)
推荐 转载

历史上的今天

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2017