第三章 MATLAB符号运算.doc

第3章 MATLAB符号计算符号计算则是可以对未赋值的符号对象(可以是常数、变量、表达式)进行运算和处理。MATLAB具有符号数学工具箱(Symbolic Math Toolbox)，将符号运算结合到MATLAB的数值运算环境。符号数学工具箱是建立在Maple软件基础上的。3.1 符号表达式的建立3.1.1 创建符号变量和表达式Symbolic Math Toolbox规定在进行符号计算时，首先要定义基本的符号对象然后才能进行符号运算。创建符号变量和符号表达式可以使用sym和syms命令。1. 使用sym命令创建符号变量和表达式语法：sym(变量,参数) %把变量定义为符号对象2.使用syms命令创建符号变量和符号表达式语法：syms(arg1, arg2, ,参数) %把字符变量定义为符号变量syms arg1 arg2 ,参数%把字符变量定义为符号变量的简洁形式说明：syms用来创建多个符号变量，这两种方式创建的符号对象是相同的。参数设置和前面的sym命令相同，省略时符号表达式直接由各符号变量组成。说明：参数用来设置限定符号变量的数学特性，可以选择为positive、real和unreal， positive 表示为“正、实”符号变量，real表示为“实”符号变量，unreal 表示为“非实”符号变量。如果不限定则参数可省略。【例3.1】创建符号变量，用参数设置其特性。 syms x y real %创建实数符号变量 z=x+i*y; %创建z为复数符号变量 real(z) %复数z的实部是实数x ans =x 【例3.2】创建符号表达式。 f1=sym(a*x2+b*x+c) f1 =a*x2+b*x+c 【例3.3】使用syms命令创建符号变量和符号表达式。 syms a b c x %创建多个符号变量f2=a*x2+b*x+c %创建符号表达式 f2 =a*x2+b*x+c 3.1.2符号表达式的代数运算符号运算与数值运算的区别主要有以下几点： 传统的数值型运算因为要受到计算机所保留的有效位数的限制，它的内部表示法总是采用计算机硬件提供的8位浮点表示法，因此每一次运算都会有一定的截断误差，重复的多次数值运算就可能会造成很大的累积误差。符号运算不需要进行数值运算，不会出现截断误差，因此符号运算是非常准确的。 符号运算可以得出完全的封闭解或任意精度的数值解。 符号运算的时间较长，而数值型运算速度快。 符号表达式的运算符和基本函数都与数值计算中的几乎完全相同。1. 符号运算中的运算符(1) 基本运算符 运算符“”，“”，“*”，“”，“/”，“”分别实现符号矩阵的加、减、乘、左除、右除、求幂运算。 运算符“.*”，“./”，“.”，“.”分别实现符号数组的乘、除、求幂，即数组间元素与元素的运算。 运算符“”，“.”分别实现符号矩阵的共轭转置、非共轭转置。(2) 关系运算符 在符号对象的比较中，没有“大于”、“大于等于”、“小于”、“小于等于”的概念，而只有是否“等于”的概念。 运算符“= =”、“=”分别对运算符两边的符号对象进行“相等”、“不等”的比较。当为“真”时，比较结果用1表示；当为“假”时，比较结果则用0表示。2. 函数运算(1) 三角函数和双曲函数三角函数包括sin、cos、tan；双曲函数包括sinh、cosh、tanh；三角反函数除了atan2函数仅能用于数值计算外，其余的asin、acos、atan函数在符号运算中与数值计算的使用方法相同。(2) 指数和对数函数指数函数sqrt、exp、expm的使用方法与数值计算的完全相同；对数函数在符号计算中只有自然对数log(表示ln)，而没有数值计算中的log2和log10。(3) 复数函数复数的共轭conj、求实部real、求虚部imag和求模abs函数与数值计算中的使用方法相同。但注意，在符号计算中，MATLAB没有提供求相角的命令。(4) 矩阵代数命令MATLAB提供的常用矩阵代数命令有diag，triu，tril，inv，det，rank， poly，expm，eig等，它们的用法几乎与数值计算中的情况完全一样。【例3.4】符号表达式f=2x2+3x+4与g=5x+6的代数运算。 f=sym(2*x2+3*x+4) f =2*x2+3*x+4 g=sym(5*x+6) g =5*x+6 f+g %符号表达式相加 ans =2*x2+8*x+10 f*g %符号表达式相乘 ans =(2*x2+3*x+4)*(5*x+6) 3.2 符号的基本运算命令1 合并同类项函数 collect格式 R = collect(S) %对于多项式S中的每一函数，collect(S)按缺省变量x的次数合并系数。R = collect(S,v) %对指定的变量v计算，操作同上。【例3.5】syms x y;R1 = collect(exp(x)+x)*(x+2)R2 = collect(x+y)*(x2+y2+1),y)R3 = collect(x+1)*(y+1),x+y)计算结果为：R1 = x2+(exp(x)+2)*x+2*exp(x)R2 = y3+x*y2+(x2+1)*y+x*(x2+1)R3 = (y+1)*x+y+1, x+y命令2 求反函数函数 finverse 格式 R = finverse(f,v) %对指定自变量v的函数f(v)求反函数说明：当v省略，则对默认的自由符号变量求反函数。【例3.6】求tex的反函数。f=sym(t*ex) %原函数 f =t*ex g=finverse(f) %对默认自由变量求反函数 g =log(x/t)/log(e) g=finverse(f,t) %对t求反函数 g =t/(ex) 程序分析：如果先定义t为符号变量，则参数t的单引号可去掉： syms t g=finverse(f,t) 命令3 复合函数计算函数 compose格式 compose(f,g) %返回复合函数fg(y)，其中f=f(x)，g=g(y)。其中符号x为函数f中由命令findsym(f) 确定的符号变量，符号y为函数g中由命令findsym(g) 确定的符号变量。compose(f,g,z) %返回复合函数fg(z)，其中f=f(x)，g=g(y)，符号x、y为函数f、g中由命令findsym确定的符号变量。compose(f,g,x,z) %返回复合函数fg(z)，而令变量x为函数f中的自变量f=f(x)。令x=g(z)，再将x=g(z)代入函数f中。compose(f,g,x,y,z) %返回复合函数fg(z)。而令变量x为函数f中的自变量f=f(x)，而令变量y为函数g中的自变量g=g(y)。令x=g(y)，再将x=g(y)代入函数f=f(x)中，得fg(y)，最后用指定的变量z代替变量y，得fg(z)。【例3.7】计算tex与ay2+by+c的复合函数。 f=sym(t*ex); %创建符号表达式 g=sym(a*y2+b*y+c); %创建符号表达式 h1=compose(f,g) %计算f(g(x) h1 =t*e(a*y2+b*y+c) h2=compose(g,f) %计算g(f(x) h2 =a*t2*(ex)2+b*t*ex+c h3=compose(f,g,z) %计算f(g(z) h3 =t*e(a*z2+b*z+c) 命令4 符号复数的共轭函数 conj格式 conj(X) %返回符号复数X的共轭复数【例3.8】X=real(X) + i*imag(X)，则conj(X)=real(X) - i*imag(X)命令5 符号复数的实数部分函数 real格式 real(Z) %返回符号复数z的实数部分命令6 符号复数的虚数部分函数 imag格式 imag(Z) %返回符号复数z的虚数部分命令7 符号表达式的展开函数 expand格式 R = expand(S) %对符号表达式S中每个因式的乘积进行展开计算。该命令通常用于计算多项式函数、三角函数、指数函数与对数函数等表达式的展开式。【例3.9】syms x y a b c tE1 = expand(x-2)*(x-4)*(y-t)E2 = expand(cos(x+y)E3 = expand(exp(a+b)3) E4 = expand(log(a*b/sqrt(c) E5 = expand(sin(2*t), cos(2*t)计算结果为：E1 = x2*y-x2*t-6*x*y+6*x*t+8*y-8*tE2 = cos(x)*cos(y)-sin(x)*sin(y)E3 = exp(a3)*exp(a2*b)3*exp(a*b2)3*exp(b3)E4 = log(a*b/c(1/2)E5 = 2*sin(t)*cos(t), 2*cos(t)2-1命令8 符号因式分解函数 factor格式 factor(X) %参量x可以是正整数、符号表达式阵列或符号整数阵列。若X为一正整数，则factor(X)返回X的质数分解式。若x为多项式或整数矩阵，则factor(X)分解矩阵的每一元素。若整数阵列中有一元素位数超过16位，用户必须用命令sym生成该元素。【例3.10】求f(x)=x3+6x2+11x-6 得因式分解。f=sym(x3-6*x2+11*x-6);factor(f) ans =(x-1)*(x-2)*(x-3) 命令9 符号表达式的分子与分母函数 numden格式 N,D = numden(A) 说明 将符号或数值矩阵A中的每一元素转换成整系数多项式的有理式形式，输出的参量N为分子阵，输出的参量D为分母。【例3.11】用numden函数来提取符号表达式和的分子、分母。 f1=sym(1/(s2+3*s+2) f1 =1/(s2+3*s+2) f2=sym(1/s2+3*s+2) f2 =1/s2+3*s+2 n1,d1=numden(f1) n1 =1d1 =s2+3*s+2 n2,d2=numden(f2) n2 =1+3*s3+2*s2d2 =s2 命令10 搜索符号表达式的最简形式函数 simple格式 r = simple(S) %该命令试图找出符号表达式S的代数上的简单形式，显示任意的能使表达式S长度变短的表达式，且返回其中最短的一个。若S为一矩阵，则结果为整个矩阵的最短形式，而非是每一个元素的最简形式。若没有输出参量r，则该命令将显示所有可能使用的算法与表达式，同时返回最短的一个。r,how = simple(S) %没有显示中间的化简结果，但返回能找到的最短的一个。输出参量r为一符号，how为一字符串，用于表示算法。【例3.12】syms xR1 = simple(cos(x)4+sin(x)4)R2 = simple(2*cos(x)2-sin(x)2)R3 = simple(cos(x)2-sin(x)2)R4 = simple(cos(x)+(-sin(x)2)(1/2)R5 = simple(cos(x)+i*sin(x)R6 = simple( (x+1)*x*(x-1)R7 = simple(x3+3*x2+3*x+1)R8,how = simple(cos(3*acos(x)计算的结果为：R1 = 1/4*cos(4*x)+3/4R2 = 3*cos(x)2-1R3 = cos(2*x)R4 = cos(x)+i*sin(x)R5 = exp(i*x)R6 = x 3-xR7 = (x+1)3R8 = 4*x3-3*xhow = expand命令11 符号表达式的化简函数 simplify格式 R = simplify(S) 说明 使用Maple软件中的化简规则，将化简符号矩阵S中每一元素。【例3.13】syms x a b cR1 = simplify(sin(x)4 + cos(x)4) R2 = simplify(exp(c*log(sqrt(a+b)S = (x2+5*x+6)/(x+2),sqrt(16);R3 = simplify(S)计算结果为：R1 = 2*cos(x)4+1-2*cos(x)2R2 = (a+b)(1/2*c)R3 = x+3, 4 3.3 符号极限、微积分和级数求和3.3.1符号极限假定符号表达式的极限存在，Symbolic Math Toolbox提供了直接求表达式极限的函数limit，函数limit的基本用法如表3.2所示。表3.2 limit函数的用法表达式函数格式说明limt(f)对x求趋近于0的极限 limt(f,x,a)对x求趋近于a的极限，当左右极限不相等时极限不存在。limt(f,x,a, left)对x求左趋近于a的极限 limt(f,x,a, right)对x求右趋近于a的极限【例3.14】分别求1/x在0处从两边趋近、从左边趋近和从右边趋近的三个极限值。 f=sym(1/x) f =1/x limit(f) %对x求趋近于0的极限 ans =NaN limit(f,x,0) %对x求趋近于0的极限 ans =NaN limit(f,x,0,left) %左趋近于0 ans =-inf limit(f,x,0,right) %右趋近于0 ans =inf 程序分析：当左右极限不相等，表达式的极限不存在为NaN。3.3.2符号求导函数diff是用来求符号表达式的导数。语法： diff(f,x)%求导 diff(f,x,n)%求n阶导数 diff(diff(f,x,m),y,n)或diff(diff(f,y,n),x,m) %多元函数的偏导数【例3.15】设， 求y syms x %定义X为符号变量 y=x10+10x+log(x) y = x10+10x+log(x) diff(y) %求ans = 10*x9+10x*log(10)+1/x【例3.16】设 ，求 syms x y;z=exp(2*x)*(x+y2+2*y);a=diff(z,x)b=diff(z,y)c=diff(z,x,2)d=diff(z,y,2)e=diff(a,y) 其中： 结果为：a =2*exp(2*x)*(x+y2+2*y)+exp(2*x) b =exp(2*x)*(2*y+2) c =4*exp(2*x)*(x+y2+2*y)+4*exp(2*x) d =2*exp(2*x) e =2*exp(2*x)*(2*y+2)3.3.3符号积分积分有定积分和不定积分，运用函数int可以求得符号表达式的积分。语法：int(f,t) %求符号变量t的不定积分int(f,t,a,b) %求符号变量t的积分int(f,t,m,n) %求符号变量t的积分说明：t为符号变量，当t省略则为默认自由变量；a和b为数值，a,b为积分区间；m和n为符号对象，m,n为积分区间；与符号微分相比，符号积分复杂得多。因为函数的积分有时可能不存在，即使存在，也可能限于很多条件，MATLAB无法顺利得出。当MATLAB不能找到积分时，它将给出警告提示并返回该函数的原表达式。【例3.17】求积分和 。 f=sym(cos(x); int(f) %求不定积分 ans =sin(x) int(f,0,pi/3) %求定积分 ans =1/2*3(1/2) int(f,a,b) %求定积分 ans =sin(b)-sin(a) int(int(f) %求多重积分 ans =-cos(x) 【练习1】已知，求。syms t xft=t2*cos(t)sx=int(ft,t,0,x) ft =t2*cos(t)sx =x2*sin(x)-2*sin(x)+2*x*cos(x) 【练习2】已知，求。syms t xft=exp(-sin(t)sx=int(ft,t,0,4) ft =exp(-sin(t)Warning: Explicit integral could not be found. In at 58sx =int(exp(-sin(t),t = 0 . 4) 3.3.4符号级数1. symsum函数语法：symsum(s,x,a,b)%计算表达式s的级数和说明：x为自变量，x省略则默认为对自由变量求和；s为符号表达式；a,b为参数x的取值范围。【例3.18】求级数和1+x+x2+xk+的和。 syms x k s1=symsum(1/k2,1,10) %计算级数的前10项和 s1 =1968329/1270080 s2=symsum(1/k2,1,inf) %计算级数和,inf:无穷大s2 =1/6*pi2 s3=symsum(xk,k,0,inf) %计算对k为自变量的级数和 s3 =-1/(x-1) 2. taylor函数语法：taylor (F,x,n) %求泰勒级数展开说明：x为自变量，F为符号表达式；对F进行泰勒级数展开至n项，参数n省略则默认展开前5项。【例3.18续】求ex的泰勒展开式为：。 syms x s1=taylor(exp(x),8) %展开前8项 s1 =1+x+1/2*x2+1/6*x3+1/24*x4+1/120*x5+1/720*x6+1/5040*x7 s2=taylor(exp(x) %默认展开前5项 s2 =1+x+1/2*x2+1/6*x3+1/24*x4+1/120*x5 上机练习1上机练习：例3.14，3.15，3.16，3.17，3.18，练习1，练习2，熟悉matlab的符号规则以及极限、导数、积分和级数的运算。3.5 符号积分变换*3.5.1傅里叶(Fourier)变换及其反变换fourier变换和反变换可以利用积分函数int来实现，也可以直接使用fourier或ifourier函数实现。1. fourier变换语法：Ffourier(f,t ,w) %求时域函数f(t)的fourier变换F说明：返回结果F是符号变量w的函数，当参数w省略，默认返回结果为w的函数；f为t的函数，当参数t省略，默认自由变量为x。2. fourier反变换语法：f=ifourier (F) %求频域函数F的fourier反变换f(t)f=ifourier (F,w,t) 说明：ifourier函数的用法与fourier函数相同。【例3.19】计算f(t)= 的fourier变换F以及F的fourier反变换。 syms t w F=fourier(1/t,t,w) %fourier变换 F =i*pi*(Heaviside(-w)-Heaviside(w) f=ifourier(F,t) %fourier反变换 f =1/t f=ifourier(F) %fourier反变换默认x为自变量 f =1/x 程序分析：其中Heaviside(t)是单位阶跃函数，函数名为数学家Heaviside的名字。【例3.19续】单位阶跃函数的fourier变换。 fourier(sym(Heaviside(t) ans =pi*Dirac(w)-i/w 程序分析：Dirac函数为单位脉冲函数。3.5.2拉普拉斯(Laplace)变换及其反变换1. Laplace变换语法：F=laplace(f,t,s) %求时域函数f的Laplace变换F说明：返回结果F为s的函数，当参数s省略，返回结果F默认为s的函数；f为t的函数，当参数t省略，默认自由变量为t。【例3.20】求sin(at)和阶跃函数的Laplace变换。 syms a t s F1=laplace(sin(a*t),t,s) %求sinat的Laplace变换 F1 =a/(s2+a2) F2=laplace(sym(Heaviside(t) %求阶跃函数的Laplace变换 F2 =1/s 2. Laplace反变换语法：filaplace(F,s,t)%求F的Laplace反变换f【例3.20续】求和1的Laplace反变换。 syms s a t f1=ilaplace(1/(s+a),s,t) %求1/s+a的Laplace反变换 f1 =exp(-a*t) f2=ilaplace(1,s,t) %求1的Laplace反变换是脉冲函数 f2 =Dirac(t) 3.5.3 Z变换及其反变换1. ztrans函数语法：Fztrans(f,n, z) %求时域序列f的Z变换F说明：返回结果F是以符号变量z为自变量；当参数n省略，默认自变量为n；当参数z省略，返回结果默认为z的函数。【例3.21】求阶跃函数、脉冲函数和e-at的Z变换。 syms a n z t Fz1=ztrans(sym(Heaviside(t),n,z) %求阶跃函数的Z变换 Fz1 =Heaviside(t)*z/(z-1) Fz2=ztrans(sym(Dirac(t),n,z) %求脉冲函数的Z变换 Fz2 =Dirac(t)*z/(z-1) Fz3=ztrans(exp(-a*t),n,z) %求e-at的Z变换 Fz3 =exp(-a*t)*z/(z-1) 2. iztrans函数语法：fiztrans(F,z,n) %求F的z反变换f【例3.21续】用Z反变换验算阶跃函数、脉冲函数和e-at的Z变换。 syms n z t f1=iztrans(Fz1,z,n) f1 =Heaviside(t) f2=iztrans(Fz2,z,n) f2 =Dirac(t) f3=iztrans(Fz3,z,n) f3 =exp(-a*t) 3.6符号方程的求解3.6.1代数方程当方程不存在解析解又无其他自由参数时，MATLAB可以用solve命令给出方程的数值解。语法：solve(eq,v) %求方程关于指定变量的解solve(eq1, eq2,v1,v2,) %求方程组关于指定变量的解 x,y,z= solve(eq1,eq2,eq3,) % 求方程组变量x,y,z的解说明：eq可以是含等号的符号表达式的方程，也可以是不含等号的符号表达式，但所指的仍是令eq=0的方程；当参数v省略时，默认为方程中的自由变量。【例3.22】求方程ax2+bx+c=0和sinx=0的解。 f1=sym(a*x2+b*x+c)%无等号 f1 =a*x2+b*x+c solve(f1) %求方程的解x ans =-(b + (b2 - 4*a*c)(1/2)/(2*a) -(b - (b2 - 4*a*c)(1/2)/(2*a) f2=sym(sin(x) f2 =sin(x) solve(f2,x) ans =0 程序分析：当sinx=0有多个解时，只能得出0附近的有限几个解。 【例3.23】求三元非