计算机仿真技术与CAD
上QQ阅读APP看本书,新人免费读10天
设备和账号都新为新人

1.6 MATLAB的符号运算

MATLAB的优点不仅在于其强大的数值运算功能,而且也在于其强大的符号运算功能。MATLAB的符号运算是通过集成在MATLAB中的符号数学工具箱(Symbolic MathToolbox)来实现的。它可完成几乎所有符号表达式的运算功能,如符号表达式的生成、复合和化简,符号矩阵的求解,符号微积分的求解,符号函数的画图,符号代数方程的求解,以及符号微分方程的求解等。

1.6.1 符号表达式的生成

在MATLAB中的符号数学工具箱中,符号表达式是指代表数字、函数和变量的MATLAB字符串或字符串数组,它不要求变量有预先确定的值。

符号表达式可以是符号函数或符号方程。其中,符号函数没有等号,而符号方程必须有等号。MATLAB在内部把符号表达式表示成字符串,以便与数字相区别。符号表达式可由以下三种方法生成。

1.用单引号生成符号表达式

在MATLAB中,符号表达式如同字符串一样也可利用单引号来直接设定。例如:

      >>fun=′sin(x)′

结果显示:

      fun=
          sin(x)
      >>fun=′a*x^2+b*x+c=0′

结果显示:

      fun=
          a*x^2+b*x+c=0

2.用sym()函数生成符号表达式

在MATLAB可自动确定变量类型的情况下,不用sym()函数来显式地生成符号表达式。但在某些情况下,特别是在建立符号数组时,必须要用sym()函数来将字符串转换成符号表达式。例如:

      >>A=sym(′[sin(x)b;c d]′)

结果显示:

      A=
        [ sin(x),b]
        [c,d]

其结果表示A为一个2 ×2维的符号矩阵。但命令

      >>A=′[sin(x)b;c d]′

的结果为

      A=
        [sin(x)b;c d]

其结果表示A为一字符串。

3.用命令syms生成符号表达式

在MATLAB中,利用命令syms只能生成符号函数,而不能生成符号方程。例如

      >>syms K t T;fun=K*(exp(-t/T))

结果显示:

      fun=
          K*exp(-t/T)

另外,在MATLAB中,利用symvar()函数可知道符号表达式中哪些变量为符号变量。同时MATLAB会自动把i,j,pi,inf,nan,eps等特殊字母不当成符号变量。例如,

      >>symvar(′5*pi+j*K*(exp(-t/T)′)

结果显示:

            ′K′
            ′T′
            ′t′

1.6.2 符号表达式的基本运算

在MATLAB的符号工具箱中,利用相关函数,可对符号表达式进行分子/分母的提取、基本代数运算、相互转换、化简和替换等基本运算。

1.符号表达式的提取分子/分母运算

在MATLAB中,如果符号表达式为有理分式的形式或可展开为有理分式的形式,则可通过numden()函数来提取符号表达式中的分子与分母。其调用格式为

[num,den] =numden(f)

其中,f表示所求符号表达式;num和den表示返回所得的分子与分母。例如:

      >>f=sym(′(x+d)/(a*x^2+b*x+c)′);[num,den] =numden(f)

运行结果:

      num=
          x+d
      den=
          a*x^2+b*x+c

2.符号表达式的基本代数运算

在MATLAB中,符号表达式的加、减、乘、除四则运算及幂运算等基本的代数运算,分别由symadd(),symsub(),symmul(),symdiv()及sympow()函数来实现。其中求和函数symadd()的调用格式为

h=symadd(f,g)

式中,f,g表示待运算的符号表达式;h表示结果符号表达式。其中,当f,g为符号矩阵时,以上四则运算及幂运算的命令仍然成立。以上其他函数的调用格式均与求和函数symadd()的调用格式一致。

3.符号表达式与数值表达式的相互转换

在MATLAB中,利用numeric()函数(仅适用于MATLAB 6.5及以前的版本)或eval()函数可将符号表达式转换成数值表达式。反之,sym()函数可将数值表达式转换成符号表达式。例如:

      >>f=′abs(-1)+sqrt(1)/2′,p=eval(f),n=sym(p)

运行结果:

      f=
        abs(-1)+sqrt(1)/2
      p=
        1.5000
      n=
        3/2

若已知数值多项式系数向量,则可以通过符号运算工具箱提供的函数poly2sym()将其转换成多项式表达式。若已知多项式表达式,则可以由函数sym2poly()将其转换成系数向量形式。它们调用格式为

f=poly2sym(p)和p=sym2poly(f)

其中,p为由多项式系数按降幂排列构成的系数向量;f为多项式表达式。

      >>syms x;p=sym2poly(x^2+3*x+2),f=poly2sym(p)

运行结果:

      p=
        1    3    2
      f=
        x^2+3*x+2

4.符号表达式的化简

在MATLAB中,simple()函数可按有关数学规则把符号表达式化简成最简形式,其调用格式为

y=simple(f)

其中,f表示化简前的符号表达式;y表示按有关规则化简后的符号表达式。例如:

      >>f=sym(′2*sin(x)*cos(x)′),y=simple(f)

运行结果:

      f=
        2*sin(x)*cos(x)
      y=
        sin(2*x)

另外,在MATLAB的符号数学工具箱中提供的符号表达式化简函数还有:

      pretty()   %将符号表达式转换成与公式编辑器显示的符号表达式相类似的形式;
      collect()  %将符号表达式的同类项进行合并;
      horner()   %将一般的符号表达式转换成嵌套形式的符号表达式;
      factor()   %对符号表达式进行因式分解;
      expand()   %对符号表达式进行展开;
      simplify() %利用各种类型的代数恒等式对符号表达式进行化简。

5.符号表达式的替换

在MATLAB的符号数学工具箱中,subexpr()函数和subs()函数可以进行符号表达式的替换。其中subexpr()函数用于把复杂表达式中所含的多个相同子表达式用一个符号代替,使其表达简洁,其调用格式为

g=subexpr(f,′S′)

其中,f,g分别表示置换前后的符号表达式;S表示置换复杂表达式中子表达式的符号变量。复杂表达式中被置换的子表达式是自动寻找的,只有比较长的子表达式才被置换,至于比较短的子表达式,即便多次重复出现,也不被置换。例如:

      >>f=solve(′x^3+a*x+1′);r=subexpr(f,′ss′)
      >>f=solve(′a*x^2+b*x+c′);r=subexpr(f,′ss′)

subs()函数除具有与subexpr()函数一样的可以用一个符号变量替换复杂表达式中所含的多个相同子表达式的作用外,还可以求解被替换的复杂符号表达式的值,其调用格式为

      g=subs(f,old,new)  %表示用new置换符号表达式f中的old后产生g;
      g=subs(f,new)     %表示用new置换符号表达式f中的自由变量后产生g。

例如,对于符号表达式fx)=asin(x)+5,可进行下列替代运算:

      >>syms a x;f=a*sin(x)+5;g1=subs(f,′sin(x)′,′y′),g2=subs(f,a,9)
      >>g3=subs(f,{a,x},{2,sym(pi/3)}),g4=subs(f,{a,x},{2,pi/3})
      >>g5=subs(subs(f,a,2),x,0:pi/6:pi),g6=subs(f,{a,x},{0:6,0:pi/6:pi})

1.6.3 符号表达式的微积分

MATLAB的符号工具箱中,符号表达式的微积分包括符号序列求和、符号极值、符号微分和符号积分等运算。

1.符号序列求和

对于求和问题,在MATLAB中可利用符号序列求和函数symsum()来实现,其调用格式为

      y=symsum(f,′x′,a,b) %求符号表达式f在指定变量x取遍[a,b]中所有整数和y;
      y=symsum(f,′x′)       %求符号表达式f在指定变量x取遍[0,x-1]中所有整数和y;
      y=symsum(f,a,b)        %求符号表达式f对独立变量从a到b的所有整数和y。

例1-23】 求的值。

MATLAB命令如下

      >>syms t k;f1=[t,k^3];f2=[1/(2*k-1)^2,(-1)^k/k];
      >>y1=simple(symsum(f1,′t′)),y2=simple(symsum(f2,1,inf))

运行结果:

      y1=
            [1/2*t*(t-1),k^3*t]
      y2=
            [1/8*pi^2,-log(2)]

2.符号极值

在MATLAB中,符号极值由函数limit()来实现,其调用格式为:

      y=limit(f,x,a)            %求符号表达式f对变量x趋于a时的极值y;
      y=limit(f,a)               %求符号表达式f对独立变量趋于a时的极值y;
      y=limit(f)                  %求符号表达式f对独立变量趋于0时的极值y;
      y=limit(f,x,a,′right′) %求符号表达式f对变量x从右边趋于a时的极值y;
      y=limit(f,x,a,′left′)  %求符号表达式f对变量x从左边趋于a时的极值y。

例1-24】 求符号表达式fx)=(a-2x2 x,x趋于a/6的极值。

MATLAB命令如下

      >>syms a x;f=(a-2*x)^2*x;y=limit(f,1/6*a)

运行结果:

      y=
        2/27*a^3

例1-25】 求二元函数的极限值。

解MATLAB命令如下

      >>syms a x y;f=exp(-1/(y^2+x^2))*sin(x)^2/x^2*(1+1/y^2)^(x+a^2*y^2);
      >>L=limit(limit(f,x,1/sqrt(y)),y,inf)

运行结果:

      L=
        exp(a^2)

3.符号微分

在MATLAB中,符号微分由diff()函数来实现。diff()函数可同时计算数值微分与符号微分,MATLAB能根据其输入参数的类型(数值或符号字符串),自动对其进行数值微分或符号微分。其调用格式为

      y=diff(f)          %求符号表达式f对独立变量的微分y;
      y=diff(f,n)       %求符号表达式f对独立变量的n次微分y;
      y=diff(f,′x′)   %求符号表达式f对变量x的微分y;

y=diff(f,′x′,n)%求符号表达式f对变量x的n次微分y。

例1-26】 求fx)=3ax-x3的一阶微分和二阶微分。

MATLAB命令如下

      >>f=′3*a*x-x^3′;df=diff(f,′x′),ddf=diff(f,′x′,2)

运行结果:

      df=
          3*a-3*x^2
      ddf=
          -6*x

对于多元函数的偏导数也可以嵌套使用函数diff()来求取。

例1-27】 已知二元函数fx,y)=(x2 -2x)e -x2-y2-xy,试求∂f/x,f/y和∂y/x

MATLAB命令如下

      >>syms x y;f=(x^2-2*x)*exp(-x^2-y^2-x*y);
      >>dfdx=simple(diff(f,x)),dfdy=diff(f,y),dydx=simple(diff(f,x))/diff(f,y)

运行结果:

      dfdx=
            -exp(-x^2-y^2-x*y)*(-2*x+2+2*x^3+x^2*y-4*x^2-2*x*y)
      dfdy=
         (x^2-2*x)*(-2*y-x)*exp(-x^2-y^2-x*y)
      dydx=
            -(-2*x+2+2*x^3+x^2*y-4*x^2-2*x*y)/(x^2-2*x)/(-2*y-x)

除了利用以上函数求解偏微分方程外,MATLAB还提供了偏微分方程工具箱,可以比较规范地求解各种常见的二阶偏微分方程。在MATLAB环境下键入pdetool,将启动偏微分方程求解界面。另外也可以利用Simulink求解微分方程。由于篇幅原因,这里不再介绍,具体内容可参考有关文献。

4.符号积分

在MATLAB中,符号积分由int()函数来实现。因为积分比微分要复杂得多,故在很多情况下,积分不一定能成功。当MATLAB进行符号积分而找不到原函数时,它将返回未经计算的函数。int()函数的调用格式为

      y=int(f)          %求符号表达式f对独立变量的不定积分y;
      y=int(f,′x′)   %求符号表达式f对变量x的不定积分y;
      y=int(f,a,b)    %求符号表达式f对独立变量从a到b的定积分y;
      y=int(f,′x′,a,b)%求符号表达式f对变量x从a到b的定积分y。

例1-28】 试求以下函数的定积分:

通过下面的MATLAB语句可求出所需函数的定积分:

      >>syms a;f=′a/sqrt(2*pi)*exp(-x^2/2)′;y=int(f,′x′,-inf,inf)
      >>syms a;y=int(′a/sqrt(2*pi)*exp(-x^2/2)′,′x′,-inf,inf)

结果显示:

      y=
        a

例1-29】 求积分方程

MATLAB命令如下

      >>syms x y z;
      >>f=int(int(int(′x^2+y^2+z^2′,′z′,sqrt(x*y),x^2*y),′y′,sqrt(x),x^2),′x′,1,2)
      >>p=eval(f)

结果显示:

      f=
        1610027357/6563700-6072064/348075*2^(1/2)+14912/4641*2^(1/4)+64/225*2^(3/4)
      p=
        224.9215

1.6.4 符号表达式的变换

MATLAB的符号工具箱中,可以进行以下几种符号表达式的变换。

1.Laplace变换及其反变换

在MATLAB中,给出了求解Laplace变换及其反变换的函数laplace()和ilaplace(),其调用格式分别为

F=laplace(f,t,s)和f=ilaplace(F,s,t)

其中,f表示时域函数f(t);t表示时间变量;F表示频域函数F(s);s表示频域变量。

      >>syms k t s;f=k*t^0;F=laplace(f,t,s),f1=ilaplace(F)

结果显示:

      F=
          k/s
      f1=
          k
      >>syms k T t s;f=k*exp(-t/T);F=laplace(f,t,s),f1=ilaplace(F)

结果显示:

      F=
          k/(s+1/T)
      f1=
          k*exp(-t/T)

2.Z变换及其反变换

在MATLAB中,给出了求解Z变换及其反变换的函数ztrans()和iztrans(),其调用格式分别为

F=ztrans(f,n,z)和f=iztrans(F,z,n)

其中,f表示时域序列f(n)或时间函数f(t);n表示时间序列;F表示Z域函数F(z);z表示Z域变量。

      >>syms k t z;f=k*t^1;F=ztrans(f,t,z),f1=iztrans(F)

结果显示:

      F=
          k*z/(z-1)^2
      f1=
          k*n

3.Fourier变换及其反变换

在MATLAB中,给出了求解Fourier变换及其反变换的函数fourier()和ifourier(),其调用格式分别为

F=fourier(f,t,ω)和f=ifourier(F,ω,t)

其中,f表示时间函数f(t);t表示时间变量;F表示频域函数F(ω);ω表示频域变量。

1.6.5 符号表达式的求解

MATLAB的符号工具箱中,符号表达式的求解包括符号代数方程的求解和符号微分方程的求解等。

1.符号代数方程求解

在MATLAB中,符号代数线性方程、符号代数非线性方程,以及符号超越方程均可利用solve()函数对其进行求解。solve()函数的调用格式为

                [x,y,z,…] =solve(′eq1′,′eq2′,′eq3′,…,′a′,′b′,′c′,…)

                [x,y,z,…] =solve(′eq1,eq2,eq3,…′,′a,b,c,…′)
                [x,y,z,…] =solve(exp1,exp2,exp3,…,′a′,′b′,′c′,…)

其中,eq1,eq2,eq3,…表示符号方程,或不含“等号”的符号表达式(或称为函数)(此时函数是对eq1=0,eq2=0,eq3=0,…求解);exp1,exp2,exp3,…仅表示符号表达式;a,b,c,…是符号方程的求解变量名;x,y,z,…是符号方程的解赋值的变量名。

例1-30】 求方程3x2 -3a2 =0的解。

MATLAB命令如下:

      >>x=solve(′3*x^2-3*a^2=0′,′x′)
或
      >>x=solve(′3*x^2-3*a^2′,′x′)
      >>f=′3*x^2-3*a^2=0′;x=solve(f,′x′)
      >>f=′3*x^2-3*a^2′;x=solve(f,′x′)

结果显示:

      x=
          -a
          a

例1-31】 求方程组的解。

MATLAB命令如下:

      >>[x,y] =solve(′3*x+y=a′,′x-y=a′,′x′,′y′)
或
      >>f=′3*x+y=a′;g=′x-y=a′;[x,y] =solve(f,g,′x′,′y′)

结果显示:

      x=
        1/2*a
      y=
        -1/2*a

例1-32】 求非线性方程组的解。

可利用以下命令:

      >>f=′sin(x)+y^2+ln(z)-7=0′;g=′3*x+2^y-z^3+1=0′;h=′x+y+z-5=0′;
      >>[x1,y1,z1] =solve(f,g,h,′x′,′y′,′z′)

结果显示:

      x1=
          0.5991
      y1=
          2.3959
      z1=
          2.0050

2.符号微分方程求解

在MATLAB中,可利用dsolve()函数对符号微分方程进行求解,其调用格式为

                [y1,y2,…] =dsolve(′eq1′,′eq2′,…,′cond1′,′cond2′,…,′x′)

                [y1,y2,…] =dsolve(′eq1,eq2,…′,′cond1,cond2,…′,′x′)
                [y1,y2,…] =dsolve(exp1,exp2,…,′cond1,cond2,…′,′x′)

其中,eq1,eq2,…表示所求符号微分方程,或不含“等号”的符号微分表达式(此时函数是对eq1=0,eq2=0,…求解的);exp1,exp2,…仅表示符号微分表达式;cond1,cond2,…表示初始条件或边界条件;x表示独立变量,当x省略时,表示独立变量为t;y1,y2,…表示输出量。微分方程的记述规定:当y是因变量时,用Dny表示y对x的n阶导数。例如,Dny表示形如的导数。

例1-33】 求微分方程,在边界条件x(0)=0,x(1)=1的解。

MATLAB命令如下

      >>f=′-2*D2x+12*a*t=0′;x=dsolve(f,′x(0)=0,x(1)=1′,′t′)

运行结果:

      x=
        a*t^3+(-a+1)*t

例1-34】 求一阶非线性微分方程的解。

MATLAB命令如下

      >>x=dsolve(′Dx=x*(1-x^2)′)

运行结果:

      x=
        [  1/(1+exp(-2*t)*C1)^(1/2)]
        [ -1/(1+exp(-2*t)*C1)^(1/2)]

即该微分方程的解为

例1-35】 求微分方程组在边界条件x1(0)=0,x1(π/2)=1,x2(0)=0,x2(π/2)=-1的解。

MATLAB命令如下

      >>f=′D2x1-x2=0,D2x2-x1=0′;
      >>[x1,x2] =dsolve(f,′x1(0)=0,x1(pi/2)=1,x2(0)=0,x2(pi/2)=-1′,′t′)

运行结果:

      x1=
            sin(t)
      x2=
          -sin(t)