第8章符号计算64.ppt

第八章,符号计算,荆华,MATLAB使用详解2012学年选修课,第八章 符号计算,与一般的数值计算不同，符号计算是对字符串符号进行分析和运算，为了便于理解，可以将符号计算看作“由计算机实现的数学公式推导”。进行符号计算时，MATLAB负责将计算请求提交给其内置的MAPLE组件并返回MAPLE的计算结果MATLAB的符号计算历经多次的改进和完善，其功能已经非常强大，尤其是在大规模的简单公式推导、逻辑推导等应用中有重要应用。,本章介绍的内容如下：符合对象的创建 符号函数符号对象的复合和化简 符号对象的运算符号矩阵符号微积分符号方程求解,第八讲 符号计算,8.1 符号类型与数值类型,符号对象（Symbolic Objet）是MATLAB符号工具箱中定义的一种新的数据类型。与一般的简单数值类型（如logical、int、float、double等）不同，符号对象是一种保存了符号字符串表示的复杂数据结构。首先，讨论符号类型与普通数值类型的差别。,例：对double型数值2求自然对数在Matlab命令窗口输入以下代码：log(2)ans=0.6931 log(sym(2)%2被定义为符号常量ans=log(2)符号计算在计算过程中没有采用近似，也没有字长限制引起的截断效应等。即符号计算能够得到问题的解析解，而不是数值近似解。,8.1 符号类型与数值类型,8.2 符号对象及其创建,MATLAB符号工具箱定义的数据类型符号对象（Symbolic Object），是一种存储了符号字符串表示的复杂数据结构。MATLAB利用符号对象表示符号变量、符号表达式、符号函数、符号矩阵。,8.2.1 创建符号变量,sym函数可用于创建符号常量和变量，其调用格式为：sym_variable=sym(sym_expression)该指令创建符号变量（或常量）sym_variable，其中sym_expression为符号变量的字符串表示。如果sym_expression为常数（例如pi、1、2.01等），则函数返回符号常量；如果sym_expression是符合命名规则的字符串（即该字符串只能包含数字、英文字母和下画线，且不能以数字开头），则函数返回符号变量。,另外，对符号变量的创建，Matlab提供了syms指令，syms可以一次创建多个符号变量。syms调用格式为:syms var1 var2 var3注意：变量之间只能以空格分开，不能是“，”或“；”。例：创建符号变量 syms a b c 不能是syms a,b,c,8.2.1 创建符号变量,8.2.2 创建符号表达式,符号表达式是符号变量或常量的组合，也可以将符号变量和符号常量认为是符号表达式。符号表达式的创建，依然使用函数sym，其调用格式为：sym_variable=sym(sym_expression)这里的输入字符串参数sym_expression是由多个子字符串通过数值运算符组合而成的。,例：创建符号表达式ax2+bx+c。有两种创建符号表达式的方法，第一种是直接输入法：f=sym(a*x2+x*x+c)f=a*x2+x*x+c第二种是先创建符号变量a、b、c、x，然后用数值运算符组合这些符号变量得到符号表达式。syms a b c x;f=a*x2+b*x+c f=a*x2+b*x+c,应用举例：,8.2.3 符号对象与数值的转换,MATLAB符号计算工具箱提供如下三类运算操作：Numberic，即MATLAB浮点运算；Rational，即Maple符号运算；VPA，即Maple精度可变运算。浮点运算：运算速度最快、内存使用最少，但其结果的精度也较低。符号运算：运算速度慢、内存使用最多、但结果是精确、无近似的。VPA是浮点运算与符号运算的折中，用户可以通过函数digits控制计算结果的精度，默认情况下digits=32。,求一元二次方程x2+x-1=0的数值解和解析解。在Matlab命令窗口中输入以下代码：sym x；%定义符号变量x s=solve(x2+x-1)%求解析解s=1/2*5(1/2)-1/2-1/2-1/2*5(1/2)%精确解，符号计算在计算过程中没 有采用近似 vpa(s,4)%将解析解转换为数值解，小数位个数是4ans=.6180-1.618,应用举例：,应用举例*：,例：分别用浮点运算、符号运算和VPA计算1/2+1/3得值,format long 1/2+1/3ans=0.83333333333333 sym(1/2)+sym(1/3)ans=5/6 digits(32);vpa(1/2+1/3)ans=.83333333333333333333333333333333,8.3 创建符号函数,可以通过下面两种方法创建符号函数：符号表达式；M文件。利用符号表达式创建符号函数的方法已经介绍，这种方法适合于创建一些比较简单的函数，对于比较复杂的符号函数的创建，应该利用M文件。注意：符号运算不支持函数句柄作为输入参数。,8.3.1 复合函数,已知函数f(x)、g(x)，复合函数f.g定义为(f.g)(z)=f(g(z)。MATLAB符号计算工具箱提供函数compose用于符号函数复合，其调用格式为：compose(f,g)compose(f,g,z)compose(f,g,x,z),compose(f,g,x,y,z)其中f、g为符号表达式，compose(f,g,x,y,z)首先将g中的y替换为z得到g(z)，然后将f中的x替换为g(z)。如果x、y未指定，则xfindsym(f,1)、yfindsym(g,1);如果z未指定，则zy。,8.3.1 复合函数,应用举例：,已知f=(1+x)2、g=tan(y)、hx*t、pe-y/u1)创建符号变量及符号函数；2）分别求compose(f、g)、compose(f、g、t);3）分别求compose(h、g、x、z)、compose(h、g、t、z)在Matlab中输入如下代码：,syms x y z t u;%创建符号变量及符号函数 f=(1+x)2;g=tan(y);h=x*t;p=exp(-y/u);compose(f,g)%将f中的x用g（y）代替ans=(1+tan(y)2 compose(f,g,t)%将f中的x用g（t）代替ans=(1+tan(t)2,compose(h,g,x,z)%将h中的x替换为g（z）ans=tan(z)*t compose(h,g,t,z)%将h中的t替换为g（z）ans=x*tan(z),函数findsym,findsym用于查询符号表达式中的变量，其使用方法如下：findsym(S)返回S中所有的符号变量。findsym(S,n)依Ascll码表顺序返回S中n个距x最近的符号变量。例：已知S=ax+b,利用findsym查询S的符号变量 syms a b x;s=a*x2+b;findsym(s)ans=a,b,x,8.3.2 反函数,MATLAB符号工具箱的函数finverse用于计算函数的反函数，其调用格式为：g=finverse(f)g=finverse(f,v)其中f为符号表达式，该函数求f关于v的反函数，如果未指定v，那么v=findsym(f)。,应用举例：,例：求函数f(x)=e-x、g(x)=x+e-x的反函数在Matlab命令窗口中输入以下代码：syms x;f=exp(-x);g=x+exp(-x);finverse(f)%f为符号表达式ans=-log(x)finverse(g)ans=lambertw(-1/exp(x)+x,8.4 符号表达式的复合和化简,符号表达式的显示、展开、分解、化简等是符号运算的重要内容，MATLAB符号工具箱对这些内容提供了强大的支持。本节将详细介绍MATLAB符号工具箱中符号表达式的显示、展开、分解、替换、化简等操作。,8.4.1 显示符号表达式,MATLAB显示符号表达式时，默认的方式是MATLAB式的显示，例如在MATLAB命令窗口输入以下代码。f=a*x2+b*x+c运行结果如下。f=a*x2+b*x+c,8.4.2 展开、分解及化简符号表达式,对数学表达式的展开、分解、化简是学习数学和应用数学中经常要做的工作，下面列出了一些相关的函数及其使用方法。函数collect用于同类项合并，其调用格式为：collect(f,s)collect(f)collect(f,s)将符号表达式f看作符号变量s的多项式（如果未指定s，则默认s为findsym(f,1)），对具有相同阶数的项进行合并。,函数expand对符号表达式进行展开，其调用格式为R=expand（S）函数factor对符号表达式进行因式分解，其调用格式为：Rfactor(X),其中X可以为符号表达式或正整数符号常量。,8.4.2 展开、分解及化简符号表达式,应用举例：,例1：对符号表达式f(x)=x2-x-3x+3合并同类项在Matlab命令窗口中输入：syms x;f=x2-x-3*x+3 collect(f)ans=x2-4*x+3,应用举例：,例：展开符号表达式：f(x)=(x-1)(x-3)syms x;f=(x-1)*(x-3)expand(f)ans=x2-4*x+3,应用举例,对符号表达式3+x24x作因式分解 syms x;factor(3+x2-4*x)ans=(x-1)*(x-3),应用举例*：,例：合并符号表达式 1)x2y+yx-x2-2x 2)4xe(-2x)+3e(-2x);,syms x y collect(x2*y+y*x-x2-2*x)ans=(y-1)*x2+(y-2)*x f=4*x*exp(-2*x)+3*exp(-2*x);collect(f)ans=4*x*exp(-2*x)+3*exp(-2*x),8.4.3 符号表达式替换,符号表达式的替换有两种，一是利用符号变量或符号表达式替换目标符号表达式中的符号变量（包括符号常量），利用函数subs实现；二是用符号变量替代目标符号表达式中的某个子符号表达式，利用函数subexpr实现。函数subs的一般调用格式为：R=subs(S)R=subs(S,new)R=subs(S,old,new),应用举例：,例：已知f(x)=ex,a=logb,A=1 2,求f(a)、f(A)syms x b a;f=exp(x);a=log(b)a=log(b)subs(f,a)ans=b A=1,2;3,1,subs(f,A)ans=2.7183 7.3891 20.0855 2.7183,应用举例*：,例：求符号表达式x2+5*x*sqrt(x)当x=3的值,syms x f=x2+5*x*sqrt(x)subs(f,x,3)ans=34.98076211353316,8.5 符号矩阵,符号矩阵是数值矩阵的扩展，符号矩阵的每个元素为符号表达式。与数值矩阵一样，符号矩阵也有代数运算、逻辑运算等。对于符号方阵，可以求行列式、特征值；符号方阵可逆时，还可以求符号方阵的逆。本节涉及的许多符号矩阵运算与数值矩阵相同。,8.5.1 创建符号矩阵,MATLAB符号工具箱中没有提供专门用于创建符号矩阵的函数，用户可以通过以下三种方式创建符号矩阵：直接输入符号矩阵的各个元素符号表达式；利用sym函数将数值矩阵转化为符号矩阵，这种方法只能创建符号常量矩阵；简单符号矩阵的组合创建复杂符号矩阵。,例：利用sym函数将数值矩阵转化为符号矩阵 V=vander(1,2);V=sym(V)V=1,1 2,1例：直接输入法创建符号矩阵 syms x;s=sin(x),cos(x);-cos(x),sin(x)s=sin(x),cos(x)-cos(x),sin(x),8.5.1 创建符号矩阵,9.5.2 符号矩阵的代数运算,同数值矩阵一样，符号矩阵也能进行加、减、乘、除四则运算，并且具有与数值矩阵几乎相同的运算规则，只是矩阵元素由数值变成了符号表达式。,应用举例：,已知A=a,b;c,d,x=u,v,求A*x、Ax、A/x syms a b c d u v;A=a,b;c,d;x=u;v;A*x ans=a*u+b*v c*u+d*v Ax ans=(d*u-b*v)/(-c*b+a*d)-(c*u-a*v)/(-c*b+a*d)A/x?Error using=sym.mrdivide Second dimensions must agree.矩阵右除要求x的列数与A的行数相等。,9.5.3 符号矩阵的逻辑运算,MATLAB符号工具箱对符号的逻辑运算支持很有限，仅支持=和=。不能对两个符号表达式作大于或小于关系运算，否则Matlab会给出错误信息。例：syms a b;(a+b)2a2+2*a*b+b2运行结果如下。?Function gt is not defined for values of class sym.Error in=gt at 16 builtin(gt,varargin:);,9.5.4 符号矩阵的行列式,符号矩阵的行列式与数值矩阵的行列式具有相同的定义，只是符号矩阵的行列式为符号函数。同样，行列式仅对符号方阵才有意义。MATLAB利用同名函数det计算符号矩阵的行列式，该函数在符号工具箱中被重载。,应用举例：,例：矩阵 被称为单位旋转矩阵，求该矩阵的行列式。syms x;det(sin(x),cos(x);-cos(x),sin(x)ans=sin(x)2+cos(x)2,9.5.5 符号矩阵的逆,对阶符号方阵，满足的符号方阵称为的逆，记为。其中为单位符号方阵，对角线上元素为符号常量，其余为0。MATLAB符号工具箱提供同名函数inv用于计算符号方阵的逆，这里的inv在符号工具箱中被重载。下面是符号矩阵求逆应用的例子。例：已知符号矩阵Q=sin(x)cos(x);-cos(x)sin(x)B=a;b 求X。,应用举例：,syms x a b;Q=sin(x)cos(x);-cos(x)sin(x);B=a;b;inv(Q)ans=sin(x)/(sin(x)2+cos(x)2),-cos(x)/(sin(x)2+cos(x)2)cos(x)/(sin(x)2+cos(x)2),sin(x)/(sin(x)2+cos(x)2)X=ans*B X=sin(x)/(sin(x)2+cos(x)2)*a-cos(x)/(sin(x)2+cos(x)2)*b cos(x)/(sin(x)2+cos(x)2)*a+sin(x)/(sin(x)2+cos(x)2)*b simplify(X)ans=sin(x)*a-cos(x)*b cos(x)*a+sin(x)*b,9.5.6 符号矩阵的秩,符号矩阵秩的含义与数值矩阵相似，对符号矩阵A，取其任意 r 行 r 列得到 A 的 r 阶子式M，A的秩定义为使“至少有一个阶子式的行列式不为0”的最大 r，即所有的 r+1 阶子式行列式为0，且至少有一个 r 阶子式行列式不为0。MATLAB提供同名函数rank计算符号矩阵的秩，这里的rank在符号工具箱中被重载。,应用举例：,例：求符号矩阵A=a,b;c,d和B=a,b;c,d;e,f的秩。syms a b c d e f;A=a,b;c,d;B=a,b;c,d;e,f;rank(A)ans=2 rank(B)ans=2,9.5.7 符号矩阵的特征分解,特征值反映了符号矩阵的一些重要的性质，对符号矩阵特征值的分析在信号处理、自动控制中有重要的应用。符号方阵A的特征值是满足的|E-A|=0的，其中E为单位符号矩阵。由于|E-A|=0即为A的特征多项式，所以A的特征值即为特征多项式的根。利用poly(A)可以得到符号方阵的特征多项式，对特征多项式求根：solve(poly(A)即可得到符号方阵A的特征值。另外，可利用函数eig计算符号矩阵的特征值。其调用格式为：eig(A),返回结果为N*1阶向量，向量元素为A的特征值,应用举例：,已知A=a,b,0;c,d,0;0,0,e，求A的特征值。syms a b c d;A=a,b,0;c,d,0;0,0,e;solve(poly(A)ans=e 1/2*d+1/2*a+1/2*(d2-2*d*a+a2+4*b*c)(1/2)1/2*d+1/2*a-1/2*(d2-2*d*a+a2+4*b*c)(1/2)eig(A)ans=e 1/2*d+1/2*a+1/2*(d2-2*d*a+a2+4*b*c)(1/2)1/2*d+1/2*a-1/2*(d2-2*d*a+a2+4*b*c)(1/2),9.6 符号微积分,前面几讲已经介绍了数值微积分，数值微积分的适用范围广泛，但对问题的描述一般是近似的，而符号微积分可能得到问题的精确描述，因此在条件允许的情况下利用符号微积分解决理论和实际问题是一个很好的选择。本小节内容包括符号的极限、微分、积分以及级数求和等内容。,9.6.1 极限,极限是当变量无限接近特定值时函数的值，例如，一元函数的导数的定义为下面的极限：。MATLAB符号工具箱利用函数limit计算符号极限，,应用举例：,已知函数f(x)=(sin(x)-x)/x求下列极限limf(x)、Limf(x)、limf(x)syms x;f=(sin(x)-x)/x;limit(f,x,inf)ans=-1 limit(f,x,0,right)ans=0,x,x0+,x0-,limit(f,x,0,left）ans=0,9.6.2 微分,照片第256页,9.6.2 微分,应用举例：,已知f(x)=eax+b，F(x,y)=ax+by;by，求df/dx、df/da、df/db、d2f/dx2和Jocabian矩阵J=F/(x,y),diff(f,x)ans=a*exp(a*x)diff(f,a)ans=x*exp(a*x)diff(f,b)ans=1 J=jacobian(F,x,y)J=a,b 0,b J=jacobian(F)J=x,y,a,b 0,y,0,b,syms x y a b f=exp(a*x)+b f=exp(a*x)+b F=a*x+b*y;b*y F=a*x+b*y b*y,9.6.3 积分,函数积分是微分的反运算，对符号函数f，f的积分是使diff(F)=f的F，用F=int(f)表示。同diff相似，如果f为多变量符号表达式，用户应该指定对某个自变量求积分，即F=int(f,v)；如果v未指定，则v设定为由函数fingsym(f,1)决定的默认符号变量。,应用举例：,例：求下列积分1/xdx、cos(ax+b)dx、cosx2dx syms x a b n;int(1/x)ans=log(x)int(cos(a*x+b),x)ans=1/a*sin(a*x+b)int(cos(x2)ans=1/2*2(1/2)*pi(1/2)*FresnelC(2(1/2)/pi(1/2)*x),9.6.5 级数求和,应用举例：,例：symsum求级数：k2 syms k;symsum(k2,k,0,k-1)ans=1/3*k3-1/2*k2+1/6*k,k=0,n-1,9.7 符号方程求解,方程求解在数学理论研究、实际应用中都是一类非常重要的问题，也是符号运算关注的一个主要内容。MATLAB符号工具箱对符号方程求解提供了强大的支持。本节主要介绍一般代数方程、线性代数方程组、常微分方程的符号求解。,9.7.1 一般代数方程,应用举例：,例：求一元二次方程的ax2+bx+c=0符号解 solve(a*x2+b*x+c=0,x)ans=1/2/a*(-b+(b2-4*a*c)(1/2)1/2/a*(-b-(b2-4*a*c)(1/2),应用举例：,例：求代数方程组 的解 x,y=solve(x2+y2=a2,x*y=b,x,y)x=b/(1/2*(a2+2*b)(1/2)+1/2*(a2-2*b)(1/2)b/(1/2*(a2+2*b)(1/2)-1/2*(a2-2*b)(1/2)b/(-1/2*(a2+2*b)(1/2)+1/2*(a2-2*b)(1/2)b/(-1/2*(a2+2*b)(1/2)-1/2*(a2-2*b)(1/2)y=1/2*(a2+2*b)(1/2)+1/2*(a2-2*b)(1/2)1/2*(a2+2*b)(1/2)-1/2*(a2-2*b)(1/2)-1/2*(a2+2*b)(1/2)+1/2*(a2-2*b)(1/2)-1/2*(a2+2*b)(1/2)-1/2*(a2-2*b)(1/2),应用举例：,并不是所有的代数方程都存在解，当解不存在时，solve返回空的符号，并且Matlab会给出警告信息。例：对ex0求解 solve(exp(x)Warning:Explicit solution could not be found.In solve at 140 ans=empty sym,9.7.2 线性代数方程组,应用举例：,例：265页令A=a,-b,b;b,a,a，B=c;d,x=x1;x2;x3，上面的线性方程组简化为AxB，对该方程的求解要分别求通解和特解。syms a b c d;A=a,-b,b;b,a,a;B=c;d;rank(A,B)，rank(A)ans=ans=%所以rank(A,B)=rank(A)2 2 x_s=AB%求特解x_s=(b*d+c*a)/(a2+b2)(d*a-c*b)/(a2+b2)0,应用举例：,x_g=null(A)%利用null函数求通解x_g=-2*b*a/(a2+b2)-(a2-b2)/(a2+b2)1 syms r;x=r*x_g+x_s%解为通解+特解x=-2*r*b*a/(a2+b2)+(b*d+c*a)/(a2+b2)-r*(a2-b2)/(a2+b2)+(d*a-c*b)/(a2+b2)r A*x-B%验证解xans=a*(-2*r*b*a/(a2+b2)+(b*d+c*a)/(a2+b2)-b*(-r*(a2-b2)/(a2+b2)+(d*a-c*b)/(a2+b2)+b*r-c b*(-2*r*b*a/(a2+b2)+(b*d+c*a)/(a2+b2)+a*(-r*(a2-b2)/(a2+b2)+(d*a-c*b)/(a2+b2)+a*r-d simplify(ans)%化简ans=0 0,9.7.3 符号常微分方程,注意：微分方程eqn不能是符号函数,应用举例：,例：267在Matlab命令窗口中输入如下代码：dsolve(Dy+a*y=0,y(0)=1,x)ans=exp(-a*x)dsolve(Dy+a*y=0,y(0)=1)ans=%若不声明,默认自变量为t exp(-a*t)与solve函数不同的是函数dsolve不支持调用格式dsolve(Dy+a*y)。在Matlab命令窗口输入以下代码：dsolve(Dy+a*y)?Undefined function or variable Dy.,应用举例：,注意：上面得到的微分方程符号解虽然是t或其它指定的自变量v的函数，但是由于t和v未被声明符号变量，对其求关于t或v的微分时会出错。在Matlab命令窗口输入以下代码：diff(ans,t)?Input argument“t is undefined.syms t;%定义t为符号变量 dy=diff(ans,t)dy=-a*exp(-a*t),注意：微分方程eqn不能是符号函数,9.8 小结,符号计算的目的是实现数学推导的自动化，这对某些规模较大的问题非常有用，可以大大减少工作量；不同于数值计算，符号计算可能得到问题的精确解，可以应用于某些计算精度要求很高的场合。本节着重从应用的角度介绍MATLAB的符号计算，内容由浅入深，通过本节学习应掌握以下的内容。符号计算与数值计算的差别符号对象（包括符号变量、符号表达式、符号函数、符号矩阵）的创建符号表达式的复合、化简符号函数符号矩阵运算符号微积分符号方程求解,11.3.4 特殊符号函数,MATLAB符号工具箱对一些重要的数学函数提供了支持，,9.5.8 符号矩阵的SVD分解,关于符号矩阵SVD分解的相关概念，可以参考第7章数值矩阵SVD分解，这里不再赘述。需要强调的是，MATLAB符号工具箱仅支持常量符号矩阵（元素为符号常量的符号矩阵）的SVD分解，这其中的原因可能是对一般的符号矩阵来说，SVD分解过于复杂而变得不可用。MATLAB中利用同名函数svd对符号矩阵进行SVD分解，其一般的调用格式为：U，S，V=svd(A),11.6.4 Taylor展开,