zoukankan      html  css  js  c++  java
  • MATLAB符号运算(3)

    3.2.5  Taylor级数

    命令1  符号函数的Taylor级数展开式

    函数  taylor

    格式  r = taylor(f,n,v)   %返回符号表达式f中的、指定的符号自变量v(若表达式f中有多个变量时)的n-1阶的Maclaurin多项式(即在零点附近v=0)近似式,其中v可以是字符串或符号变量。

    r = taylor(f)      %返回符号表达式f中的、符号变量v6阶的Maclaurin多项式(即在零点附近v=0)近似式,其中v=findsym(f)。

    r = taylor(f,n,v,a)   %返回符号表达式f中的、指定的符号自变量vn-1阶的Taylor级数(在指定的a点附近v=a)的展开式。其中a可以是一数值、符号、代表一数字值的字符串或未知变量。我们指出的是,用户可以以任意的次序输入参量nva,命令taylor能从它们的位置与类型确定它们的目的。解析函数f(x)在点x=aTaylor级数定义为:

    3-46

    >>syms x y a pi m m1 m2

    >>f = sin(x+pi/3);

    >>T1 = taylor(f)

    >>T2 = taylor(f,9)

    >>T3 = taylor(f,a) 

    >>T4 = taylor(f,m1,m2)

    >>T5 = taylor(f,m,a)

    >>T6 = taylor(f,y)

    >>T7 = taylor(f,y,m)   % 或taylor(f,m,y)

    >>T8 = taylor(f,m,y,a)

    >>T9 = taylor(f,y,a)

    计算结果为:

    T1 =

    1/2*3^(1/2)+1/2*x-1/4*3^(1/2)*x^2-1/12*x^3+1/48*3^(1/2)*x^4+1/240*x^5

    T2 =

    1/2*3^(1/2)+1/2*x-1/4*3^(1/2)*x^2-1/12*x^3+1/48*3^(1/2)*x^4+1/240*x^5-1/1440*3^(1/2)* x^6-1/10080*x^7+1/80640*3^(1/2)*x^8

    T3 =

    sin(a+1/3*pi)+cos(a+1/3*pi)*(x-a)-1/2*sin(a+1/3*pi)*(x-a)^2-1/6*cos(a+1/3*pi)*  (x-a)^3+1/24*sin(a+1/3*pi)*(x-a)^4+1/120*cos(a+1/3*pi)*(x-a)^5

    T4 =

    sin(m2+1/3*pi)+cos(m2+1/3*pi)*(x-m2)-1/2*sin(m2+1/3*pi)*(x-m2)^2-1/6* cos(m2+1/3*pi)*(x-m2)^3+1/24*sin(m2+1/3*pi)*(x-m2)^4+1/120* 

    cos(m2+1/3*pi)*(x-m2)^5

    T5 =

    sin(a+1/3*pi)+cos(a+1/3*pi)*(x-a)-1/2*sin(a+1/3*pi)*(x-a)^2-1/6*cos(a+1/3*pi)*  (x-a)^3+1/24*sin(a+1/3*pi)*(x-a)^4+1/120*cos(a+1/3*pi)*(x-a)^5

    T6 =

    sin(y+1/3*pi)+cos(y+1/3*pi)*(x-y)-1/2*sin(y+1/3*pi)*(x-y)^2-1/6*cos(y+1/3*pi) *(x-y)^3+1/24*sin(y+1/3*pi)*(x-y)^4+1/120*cos(y+1/3*pi)*(x-y)^5

    T7 =

    sin(m+1/3*pi)+cos(m+1/3*pi)*(x-m)-1/2*sin(m+1/3*pi)*(x-m)^2-1/6*cos(m+1/3*pi) *(x-m)^3+1/24*sin(m+1/3*pi)*(x-m)^4+1/120*cos(m+1/3*pi)*(x-m)^5

    T8 =

    sin(a+1/3*pi)+cos(a+1/3*pi)*(x-a)-1/2*sin(a+1/3*pi)*(x-a)^2-1/6*cos(a+1/3*pi)*  (x-a)^3+1/24*sin(a+1/3*pi)*(x-a)^4+1/120*cos(a+1/3*pi)*(x-a)^5

    T9 =

    sin(a+1/3*pi)+cos(a+1/3*pi)*(x-a)-1/2*sin(a+1/3*pi)*(x-a)^2-1/6*cos(a+1/3*pi)*  (x-a)^3+1/24*sin(a+1/3*pi)*(x-a)^4+1/120*cos(a+1/3*pi)*(x-a)^5

    命令2  Taylor级数计算器

    函数  taylortool

    格式  taylortool   %该命令生成一图形用户界面,显示缺省函数f=x*cos(x)在区间[-2*pi2*pi]内的图形,同时显示函数f的前N=7项的Taylor多项式级数和(在a=0附近的)图形,如图1。通过更改f(x)项可得不同的函数图形。

    taylortool('f')  %对指定的函数f,用图形用户界面显示出Taylor展开式。(图3-14

    3-47

    >>taylortool('sin(x*sin(x))')

    再通过改变相关的参量,可得如图3-15

      

             

         

    3-14  Taylor级数计算器                    图3-15  函数sin(x*sin(x))的taylortool界面

    3.2.6  其它

    命令1  Jacobian矩阵

    函数  jacobian

    格式  R = jacobian(w,v)   

    说明  计算wvJacobian矩阵。其中w为符号单值函数表达式或符号列向量,v为一符号行向量。输出参量R=rij)的元素rij 为:,i=1,2,…,size(w),j=1,2,…,length(v)

    3-48

    >>syms x y z u v w

    >>w = [x*y*z; y; x+z];

    >>v = [x,y,z];        

    >>R = jacobian(w,v)

    >>b = jacobian(x+u, v)

    计算结果为:

    R =

        [ y*z,  x*z, x*y]

        [   0,   1,   0]

        [   1,   0,   1]

    b =

        [ 1, 0, 0]

    命令 Jordan标准形

    函数  jordan

    格式  J = jordan(A)   %计算矩阵A的Jordan标准形。其中A为一确切已知的符号或数值矩阵。即它的元素必须是整数或小整数的比值。任何的矩阵输入误差将导致不同的Jordan标准形。即Jordan标准形对数据是敏感的。

    [V,J] = jordan(A)   %返回Jordan标准形矩阵J与相似变换矩阵V,其中V的列向量为矩阵A的广义特征向量。它们满足:V\A*V=J

    3-49

    >>A = [1 -3 -2; -1  1 -1; 2 4 5]

    >> [V,J] = jordan(A)

    >>V = double(V);

    >>Test = all(all(V\A*V == J))

    计算结果为:

    V =

        -1    -1     1

         0    -1     0

         1     2     0

    J =

         3     0     0

         0     2     1

         0     0     2

    Test = 1

    命令 Lamber的W函数

    函数  lambertw

    格式  Y = lambertw(X)   %计算参量X的每一元素x的Lamber的W函数值,其中X为一数值或符号矩阵。Lamber的函数W=W(x)为方程的解:wew = x。

    3-50

    >>W1 = lambertw([ -exp(-1); pi]) 

    >>syms x y

    >>W2 = lambertw([0 x;1 y])

    计算结果为:

    W1 =

         -1.0000 + 0.0000i

          1.0737        

    W2 =

          [         0, lambertw(x)]

          [ lambertw(1), lambertw(y)]

    命令 符号表达式的LaTex的表示式

    函数  latex

    格式  latex(S)   %返回符号表达式S的LaTex格式的表示式。该格式可以使表达式S在图形窗口中进行显示(如命令title、text等)。

    3-51

    >>syms x

    >>f = taylor(sin(1+x)); 

    >>Lat1 = latex(f) 

    >>M = sym(magic(3)); 

    >>Lat2 = latex(M)

    计算结果为:

    Lat1 =

    \sin(1)+\cos(1)\mbox {{\tt `x~`}}-1/2\,\sin(1){\mbox {{\tt `x~`}}}^{2}-1/6\,\cos(1){\mbox {{\tt `x~`}}}^{3}+1/24\,\sin(1){\mbox {{\tt `x~`}}}^{4}+{\frac {1}{120}}\,\cos(1){\mbox {{\tt `x~`}}}^{5}

    Lat2 =

    \left [\begin {array}{ccc} 8&1&6\\\noalign{\medskip}3&5&7\\\noalign{\medskip}4…

    &9&2\end {array}\right ]

    命令 调用Maple内核

    函数  maple

    格式  r = maple('statement')   %将参数命令statement传递给Maple内核,且返回计算结果。在必要时,可以在参量statement后面加上分号(;)。

    r = maple('function',arg1,arg2,  %该命令接受任何的带引号的函数名'function',与相关的输入参量arg1,arg2,…。在必要时,要将输入参量转换成符号表达式。若输入参量为syms,则maple返回一sym,否则返回一类型为char的结果。

    [r, status] = maple()   %有条件地返回警告/错误信息。当语句能顺利执行,则r为计算结果,status0;若语句不能通过执行,r为相应的警告/错误信息,而status为一正整数。

    maple('traceon') maple traceonmaple trace on   %将显示所有的后面的Maple语句与其相应的结果显示于屏幕上

    maple('traceoff') maple traceoffmaple trace off   %将关闭上面的操作特性

    3-52

    >>Pi = maple('evalf(Pi,100)')

    >>syms x

    >>v = [x^2-1;x^2-4]

    >>maple traceon

    >>w = factor(v) 

    计算结果为:

    Pi =

    3.1415926535897932384626433832795028841971693993751058209749445923078164…

                                               06286208998628034825342117068

    v =

        [ x^2-1]

        [ x^2-4]

    statement:

       map(ifactor,array([[x^2-1],[x^2-4]]));

    result:

       Error, (in ifactor) invalid arguments

    statement:

       map(factor,array([[x^2-1],[x^2-4]]));

    result:

       matrix([[(x-1)*(x+1)], [(x-2)*(x+2)]])

    w =

    [ (x-1)*(x+1)]

    [ (x-2)*(x+2)]

    命令 初始化Maple内核

    函数  mapleinit

    格式  mapleinit   该命令用于确定包含Maple库的路径,再装载Maple的线性代数与积分变换包、初始化命令digits、指定几个别名。用户可以编辑mapleinitM-文件,用于改变到Maple包的路径,只需按如下的方法改变变量initstring的值:

    1.若用户已经有一Maple V,Release 5的库在目录C:\Maple\Lib上,在文件mapleinit.m中加入:maplelib = 'C:\MAPLE\LIB'

    2.从MATLAB中删除旧的Maple包版本。

    命令 Maple数学函数的数值计算

    函数  mfun

    格式  Y = mfun('function',par1,par2,par3,par4) 

    说明  计算一指定的Maple软件中已知的数学函数function的数值。每一参量par为该函数相应的具体数值。用户可以输入满4个参量。最后指定的参量可以是矩阵,通常对应于x。其他参量的位数取决于该函数规定的范围。用户可以通过下面的命令获得相关参数的信息:help mfunlist;mhelp function;Maple用16位精度计算函数function。函数function中的任何奇异值将返回NaN。

    3-53

    >>M1 = mfun('dilog',5) 

    >>M2 = mfun('Psi',[3*i 0])

    计算结果为:

    M1 =

         -2.3699

    M2 =

         1.1080 + 1.7375i      NaN

    命令8  列出命令mfun中特定的Maple函数

    函数  mfunlist

    格式  mfunlist 

    1.列出在使用命令mfun中用到的特殊的数学函数。下表中参量的一些约定:x,y:实数参量;z,z1,z2:复数参量;m,n:整数参量

    3-1  mfun特殊函数

    函数名

    定   义

    Mfun名

    参量说明

    Bernoulli数

    与多项式

    生成函数:

    Bernoulli(n)

    Bernoulli(n,t)

    n≥0

    0<|t|<2π

    Bessel函数

    BesselI, BesselJ:第一类Bessel函数

    BesselK, BesselY:第二类Bessel函数

    BesselJ(v,x)

    BesselY(v,x)

    BesselI(v,x)

    BesselK(v,x)

    v为实数

    Beta函数

     

    Beta(x,y)

     

    二项式系数

     

    Binomial(m,n)

     

    完全椭圆积分

    第一、二、三类Legendre完全椭圆积分

    LegendreKc(k)

    LegendreEc(k)

    LegendrePic(a,k)

    a为任意实数

    -<a<∞

    k为任意实数

    0<k<1

    带余模的完全

    椭圆积分

    与余模相关的第一、二、三类Legendre完全椭圆积分

    LegendreKc1(k)

    LegendreEc1(k)

    LegendrePic1(a,k)

    a为任意实数

    -<a<∞

    k为任意实数

    0<k<1

    余差函数

    与它的累积分

    Erfc(z) = 

     

    erfc(n,z) = 

     

    erfc(z)

    erfc(n,z)

    n>0

    Dawson积分

     

    dawson(x)

     

    Ψ-函数

     

    Psi(x)

     

    二重对数积分

     

    dilog(x)

    x>1

    误差函数

            

    erf(z)

     

    Euler数与多项式

    生成Euler数的函数:

     

    euler(n)

    euler(n,z)

    n≥0

    |t|<π/2

    指数积分

     

        

    Ei(n,z)

    Ei(x)

    n≥0

    real(z)>0

    Fresnel正弦

    与余弦积分

     

    FresnelC(x)

    FresnelS(x)

     

    Г-函数

     

    GAMMA(z)

     

    调和函数

     

         =Ψ(n+1) + γ

    harmonic(n)

    n>0

    双曲正弦

    与余弦积分

     

    Chi(z) = γ+ln(z) +

           

    Shi(z)

    Chi(z)

     

    广义超几何函数

    F(n,d,z) = 

    hypergeom(n,d,x)

    其中

    n = [n1,n2,]

    d = [d1,d2,]

    n1,n2,… 为实数

    d1,d2,… 为非负实数

    不完全椭圆积分

    第一、二、三类不完全Legendre完全椭圆积分

    LegendreF(x,k)

    LegendreE(x,k)

    LegendrePi(x,a,k)

    0<x≤∞,a为实数

    -∞<a<∞,k为实数

    0<k<1

    不完全Г-函数

    Г(a,z)=

    GAMMA(z1,z2)

     

    Г-函数的对数

    lnГ(z) = ln(Г(z))

    lnGAMMA(z)

     

    对数积分

     

          = Ei(ln(x))

    Li(x)

    x>1

    Г多项式函数

    其中Ψ(z)为Γ-函数

     

    Psi(n,z)

    N0

     

    移位正弦积分

    Ssi(z)=Si(z) – π/2

    Ssi(z)

     

    对于上面的特殊函数function,用户可以通过下面的命令得到更多的帮助信息:mhelp function

    总的来说,函数的精度跟它的根相比会较低,且当它的参数相对而言较大时,精度也较底。函数的执行时间取决于特定的函数与它的输入参量。总之,其计算将比标准的MATLAB计算慢一些。

    2.正交多项式函数:

    下面的函数需要Maple正交多项式包,它们仅仅对于MATLAB的扩展符号数学工具箱有用。在使用这些函数之前,用户要用下面的命令初始化Maple正交多项式包:maple('with','orthopoly')

    3-2  正交多项式函数

    下表参量的约定:n:非负整数;x:任意实数

    多  项  式

    Maple名

    参量说明

    Gegenbauer多项式

    G(n,a,x)

    a为非有理数代数表达式或者是大于-1/2的有理数

    Hermite多项式

    H(n,x)

     

    Laguerre多项式

    L(n,x)

     

    广义Laguerre多项式

    L(n,a,x)

    a为非有理数代数表达式或者是大于-1的有理数

    Legendre

    P(n,x)

     

    Jacobi

    P(n,a,b,x)

    ab为非有理数代数表达式或者是大于-1的有理数

    第一、二类Chebyshev多项式

    T(n,x)U(n,x)

     

    命令  Maple命令帮助

    函数  mhelp

    格式  mhelp topicmhelp('topic') 

    说明  返回Maple软件中指定的Maple标题topic的在线帮助文档信息。

    命令10  交互式计算Riemann和

    函数  rsums

    格式  rsums(f)   %交互式地通过Riemann和计算函数f(x)的积分。rsums(f)显示函数的图形。用户可以通过拖动图形下方的滑块来调整Riemann和的项数,有效的项数从2128

    3-54

    >>rsums sin(-5*x^2)

    计算图形为图3-16

     

     

    3-16  函数的Riemann和

    命令11  在一符号表达式或矩阵中进行符号替换

    函数  subs

    格式  R = subs(S)   %用从调用的函数中获得的变量值,或MATLAB的工作空间中存在的变量值,替换表达式S中所有出现的相同的变量,同时自动进行化简计算;若是数值表达式,则计算出结果。

    R = subs(S,old,new)   %用新值new替换表达式s中的旧值old,参量old是一符号变量或代表一变量名的字符串,new是一符号/数值变量或表达式。若oldnew为有相同大小的阵列,则用new中相应的元素替换old中的元素;若S与old为标量,而new为阵列或单元阵列,则标量S与old将扩展为与new同型的阵列;若new为数值矩阵的单元阵列,则替换按元素的方向执行。若subs(S,old,new)没有改变S,则subs(S,old,new)被证明是可靠的。这提供了对以前版本的向后兼容性,且不会交换参量的位置。

    3-55

    >>a = 980,C1=3;

    >>y = dsolve('Dy = -a*y')

    >>syms b

    >>subs(y)

    >>subs(a+b,a,4)

    >>subs(cos(a)+sin(b),{a,b},{sym('alpha'),2})

    >>subs(exp(a*t),'a',-magic(2))

    >>subs(x*y,{x,y},{[0 1;-1 0],[1 -1;-2 1]})

    命令12  创建符号数值、变量与对象

    函数  sym

    格式  S = sym(A)   %用输入参量A,构造一类型为‘sym’的对象s。若A为字符串,则S为符号数值或变量;若A为一数值标量或矩阵,则S为代表所给数值的符号表达式。

    x = sym('x')   %创建一名字为‘x’的符号变量,且将结果存于x

    pi = sym('pi')   %创建一符号数值,这可避免了用浮点近似表示π的误差,pi的这种创建方法将暂时地代替了有相同名字、用于生成无理数π的近似值的内建数值函数pi.m。

    x = sym('x','real')   %创建一实符号变量。若x有了具体的值,则命令clear x只能清除x的值,而不能改变x的“属性”。

    x = sym('x','unreal')  %使x变成一纯粹的、没有任何附加属性的符号变量。

    S = sym(A,flag)   %将一数值标量或矩阵转换成符号形式。对浮点数值的转换方法要用第二个参量flag来指定。其中flag可以是'r''d''e''f'

    ’f’:代表“浮点格式”。

    ’r’:代表“有理格式”(该方式为缺省转换格式)。

    ’e’:代表“估计误差”。

    ’d’:代表“十进制格式”。

    命令13  创建多个符号对象的快捷命令

    函数  syms

    格式  syms  arg1  arg2 …     %定义arg1arg2为符号

    syms  arg1 arg2 … real   %该命令是下列命令的简洁形式:

    arg1 = sym('arg1','real');

    arg2 = sym('arg2','real'); 

    syms arg1 arg2 … unreal   %该命令是下列命令的简洁形式:

    arg1 = sym('arg1','unreal');

    arg2 = sym('arg2','unreal'); 

    注:clear x不能清除符号变量x的属性“real”,只能清除变量x。要想清除该属性,要输入:syms x unreal或clear mexclear all。执行后面的两个命令后,Maple内核将重新装载入MATLAB的工作空间(这是不可取的,因为花费时间)。

    3-56

    >>syms x beta real  %符号对象已经生成,执行下面一些操作:

    >>whos

    将显示工作空间中存在变量的详细信息:

          Name       Size         Bytes  Class

           beta       1x1            132  sym object

            x        1x1            126  sym object

          Grand total is 7 elements using 258 bytes

     y = x + i*beta; clear x; y 

    通过上面的操作,我们看到,当x被清除掉后,y的值并没有马上改变:

     y =

         x+i*beta

    命令14  将符号多项式转化为数值多项式

    函数  sym2poly

    格式  c = sym2poly(s)   %返回符号多项式s的数值系数行向量c。多项式自变量次数的系数按降幂排列。即行向量c的第一分量c1为多项式s的最高次数项的系数,c2为第二高次数项的系数,如此类推。

    3-57

    >>syms x u;

    >>c1 = sym2poly(3*x^3 - 2*x^2 – sqrt(5)) 

    >>c2 = sym2poly(u^4 – 3 + 5*u^2)

    计算结果为:

    c1 =

         3.0000   -2.0000    0   -2.2361

    c2 =

         1     0     5     0    -3

    命令15  可变精度算法

    函数  vpa

    格式  R = vpa(A)   %用可变精度算法来计算A中的每一元素,使其成为有d位精确度的十进制数。其中d为命令digits设置的当前位数。R中的每一元素为一符号表达式。

    R = vpa(A,d)R = vpa A d   %用参量d指定的位数(而非命令digits设置的位数)来表示A中的每一元素。R中的每一元素为一符号表达式。

    3-58

    >>digits(25)

    >>q = vpa(sym(sin(pi/6)))

    >>p = vpa(pi)

    >>gold_ratioi = vpa('(sqrt(5)-1)/2')

    >>vpa pi 75 

    >>A = vpa(gallery(5),8)

    >>B = vpa(hilb(3),5)

    计算结果为:

    q =

        .5000000000000000000000000

    p =

        3.141592653589793238462643

    gold_ratioi =

              .6180339887498948482045870

    ans =

         3.14159265358979323846264338327950288419716939937510582097… 494459230781640629

    A =

        [     -9.,     11.,    -21.,     63.,   -252.]

        [     70.,    -69.,    141.,   -421.,   1684.]

        [   -575.,    575.,  -1149.,   3451., -13801.]

        [   3891.,  -3891.,   7782., -23345.,  93365.]

        [   1024.,  -1024.,   2048.,  -6144.,  24572.]

    B =

        [     1., .50000, .33333]

        [ .50000, .33333, .25000]

        [ .33333, .25000, .20000]

    命令16  符号表达式的C语言代码

    函数  ccode

    格式  ccode(s)   %返回C语言的、用于计算符号表达式s的语句段落

    3-59

    >>syms x

    >>s = taylor(exp(x));

    >>ccode(s)

    计算结果为:

    ans =

          t0 = 1.0+x+x*x/2.0+x*x*x/6.0+x*x*x*x/24.0+x*x*x*x*x/120.0;

    注:t0x=0附近的计算公式(Taylor展式)。

    命令17  符号表达式的Fortran语言代码

    函数  fortran

    格式  fortan(s)   %返回一Fortan语言的、用于计算符号表达式s的语句段落

    3-60

    >>syms x

    >>f = taylor(sin(x));

    >>F1 = fortran(f)

    >>H = sym(hilb(4));

    >>F2 = fortran(t*(H))

    计算结果为:

    F1 =

            t0 = x-x**3/6+x**5/120

    F2 =

            T(1,1) = t      T(1,2) = t/2     T(1,3) = t/3     T(1,4) = t/4

            T(2,1) = t/2    T(2,2) = t/3     T(2,3) = t/4     T(2,4) = t/5

            T(3,1) = t/3    T(3,2) = t/4     T(3,3) = t/5     T(3,4) = t/6 

            T(4,1) = t/4    T(4,2) = t/5     T(4,3) = t/6     T(4,4) = t/7

  • 相关阅读:
    20205025模拟
    CDQ分治详解
    点分治详解
    虚树详解
    整体二分详解
    算法学习————线段树合并
    Mvc.ExceptionHandling.AbpExceptionFilter
    小程序前端转换时间格式2021-02-25T12:01:20.254748
    《生命3.0—在亿年的尺度下审视生命的演进》阅读笔记1
    软件杯赛题周总结(1)
  • 原文地址:https://www.cnblogs.com/djcsch2001/p/2333956.html
Copyright © 2011-2022 走看看