MATLAB符号计算.docx

1、MATLAB符号计算第5章 MATLAB符号计算Matlab本身并没有符号计算功能，1993年通过购买Maple的使用权后，开始具备符号运算的功能符号运算的类型很多，几乎涉及数学的所有分支符号运算使用一种特殊的数据类型，称为符号对象(Symbolic Object)，用字符串形式表达，但又不同于字符串(Char Array)。符号中的变量，函数和表达式都是符号对象。For personal use only in study and research; not for commercial use5.1 符号对象的定义5.1.1 建立符号变量MATLAB提供了两个建立符号对象的函数：sym和s

2、yms，两个函数的用法不同。For personal use only in study and research; not for commercial use1. sym函数sym函数用来建立单个符号对象，一般调用格式为：符号对象名=sym(符号字符串)该函数可以建立一个符号对象，符号字符串可以是常量、变量、函数或表达式。符号对象名=sym(num,d) 该函数将数值表达式num转化为符号表达式，数值用十进制表示；试比较分析下列命令: a=pi2, b=pi2, c=sym(pi2),d= sym(pi2,d), e=sym(b)a=pi2 %结果为一个数a = 9.8696b=pi2 %

3、结果为一个字符串b =pi2c=sym(pi2) %数值转化为符号对象，有理表示！c =*2(-49)c=sym(pi2,d) %数值转化为符号对象，十进制表示！c =9.e=sym(b) %或e=sym(pi2) 字符串转化为符号对象e =pi22. syms函数函数sym一次只能定义一个符号变量，使用不方便。MATLAB提供了另一个函数syms，一次可以定义多个符号变量（对象）。syms函数的一般调用格式为：syms 符号变量名1 符号变量名2 符号变量名n用这种格式定义符号变量时不要在变量名上加字符串分界符()，变量间用空格而不要用逗号分隔。syms x y zh=x3+2*y2+eh

4、=x3+2*y2+pi25.1.2 建立符号表达式含有符号对象的表达式称为符号表达式。建立符号表达式有以下3种方法：(1)利用单引号来生成符号表达式。(2)用sym函数建立符号表达式。(3)使用已经定义的符号变量组成符号表达式。A=c b;e-h h-x3A = *2(-49), pi2 -x3-2*y2, 2*y2+pi25.1.3 计算精度和数据类型转换 利用函数sym可以将数值表达式变换成它的符号表达式,函数numeric或eval可以将符号表达式变换成数值表达式。double(s) 将符号对象转化为双精度数值char(s) 将符号对象转化为字符串digits(n) 将数值计算精度设为n

5、位；x=vpa(s) 求s的数值结果x=vpa(s,n) 采用n位计算精度求s的数值结果210000ans = Infa=sym(2);b=a10000b =216268709376 %很长的整数，准确的，而不是近似的vpa(b)ans =.e3011 format long;pi2,format short %用长格式显示pi的值，再恢复为短格式ans = 9.108936试比较以下C=PI2的显示结果。c=sym(pi2);vpa(c,16) %16位ans =9.1089357vpa(c) %32位，默认ans =9.9998761vpa(c,100) %100位，高精度ans =9.0

6、2vpa(c,2) %2位，低精度ans =9.6double(c) %双精度，数值 ans = 9.8696syms x; f=sym(x3-1)f =x3-1fun=inline(f)fun = Inline function: fun(x) = x.3-15.2 符号矩阵和符号函数5.2.1 符号矩阵MATLAB大部分矩阵和数组运算符及指令都可以用于符号如：inv，./，.，eig等clear;A=sym(a,b;c,d);B=inv(A)B = d/(a*d-b*c), -b/(a*d-b*c) -c/(a*d-b*c), a/(a*d-b*c)A.B,ABans = d/(a*d-b

7、*c)/a, -1/(a*d-b*c) -1/(a*d-b*c), a/(a*d-b*c)/dans = (d2+b*c)/(d2*a2-2*a*d*b*c+b2*c2), -b*(a+d)/(d2*a2-2*a*d*b*c+b2*c2) -c*(a+d)/(d2*a2-2*a*d*b*c+b2*c2), (b*c+a2)/(d2*a2-2*a*d*b*c+b2*c2)A*Bans = a*d/(a*d-b*c)-b*c/(a*d-b*c), 0 0, a*d/(a*d-b*c)-b*c/(a*d-b*c)eig(A)ans = 1/2*d+1/2*a+1/2*(d2-2*a*d+a2+4*b

8、*c)(1/2) 1/2*d+1/2*a-1/2*(d2-2*a*d+a2+4*b*c)(1/2)2.符号函数计算大部分MATLAB的数学函数和逻辑关系运算也可以用于符号对象，另外还有：factor(expr) 对expr作因式分解expand(expr) 将expr展开collect(expr,v) 将expr按变量v合并同类项simplify(s)：应用函数规则对s进行化简。simple(s)：调用MATLAB的其他函数对表达式进行综合化简，并显示化简过程。g=finverse(f,v) 求函数f(v)的反函数g(v)fg=compose(f, g) 求函数f(v)和g(v)的复合函数f(

9、g(v)n,d=numden(expr) 分式通分，n返回分子，d返回分母如果符号表达式是一个有理分式或可以展开为有理分式，可利用numden函数来提取符号表达式中的分子或分母。例1. 已知 f (x,y) = (x-y)3, g(x,y) = (x+y)3，s=x6+1;考虑相关的运算，如求 h= f*g 并展开。clear;syms x y z ;f=(x-y)3;g=(x+y)3;h=f*g %两函数相乘h =(x-y)3*(x+y)3hs=expand(h) %函数展开hs =x6-3*x4*y2+3*x2*y4-y6s=x6+1;sf=factor(s) %因式分解sf =(x2+1

10、)*(x4-x2+1)sz=subs(sf,x,z2+z+1) %用z2+z+1替换sf中的xsz =(z2+z+1)2+1)*(z2+z+1)4-(z2+z+1)2+1)n,d=numden(x/y+y/x) %分式通分n =x2+y2d =x*y下面是复合函数与反函数的例子：clear;syms x y z t;f=1/(1+x2);g=sin(y);A1=compose(f,g) %求复合函数A1 =1/(sin(y)2+1)A2=compose(f,g,t) %求复合函数，变量用t表示A2 =1/(sin(t)2+1)v=finverse(f)Warning: finverse(1/(

11、x2+1) is not unique.In D:toolboxsymbolicsymfinverse.m at line 43v =1/x*(-x*(x-1)(1/2) %结果不唯一。5.3 符号表达式中变量的确定MATLAB中的符号可以表示符号变量和符号常量。findsym可以帮助用户查找一个符号表达式中的的符号变量。该函数的调用格式为：findsym(s,n)函数返回符号表达式s中的n个符号变量，若没有指定n，则返回s中的全部符号变量。符号变量确定原则(1)除了i和j之外，字母位置最接近x的字母；若距离相等，则取ASCII码大的；(2)若没有除了i 与j以外的字母，则视x为默认的符号变量

12、；(微分方程视t为默认的符号变量)syms a b x yfindsym(3*a*b+y2+1,1)ans =yfindsym(3*a*b+y2+1)ans =a, b, y5.4 符号微积分5.4.1 符号极限limit函数的调用格式为：(1)limit(f,x,a)：求符号函数f(x)的极限值。即计算当变量x趋近于常数a时，f(x)函数的极限值。(2)limit(f,a)：求符号函数f(x)的极限值。由于没有指定符号函数f(x)的自变量，则使用该格式时，符号函数f(x)的变量为函数findsym(f)确定的默认自变量，即变量x趋近于a。(3)limit(f)：求符号函数f(x)的极限值。符

13、号函数f(x)的变量为函数findsym(f)确定的默认变量；没有指定变量的目标值时，系统默认变量趋近于0，即a=0的情况。(4)limit(f,x,a,right)：求符号函数f的极限值。right表示变量x从右边趋近于a。(5)limit(f,x,a,left)：求符号函数f的极限值。left表示变量x从左边趋近于a。例5-1 求下列极限（改用课本例子）。(1) (2) （3） (4) syms h n xL=limit(log(x+h)-log(x)/h,h,0) %单引号可省略掉L =1/xM=limit(1-x/n)n,n,inf) M=exp(-x) limit(1-cos(x)/

14、(x*sin(x),x,0)ans =1/2syms x;f=(sqrt(x)-sqrt(2)+sqrt(x-2)/sqrt(x*x-4);limit(f,x,2,right)ans =1/2例5-2 （A）（B）（C）二元，其中C不存在！1. (A)求二元函数的极限： fun1=(x3*y+x*y3+x*y)/(x+y)3; limit(limit(fun1,x,2),y,1) ans =4/9（B）1. 求极限: syms m n limit(symsum(m/n2,m,1,n),n,inf) ans = 1/2（C）syms x ylimit(limit(x*y/(x2+y2),x,0)

15、,y,0) ans =0注：此二元函数的极限不存在，但计算结果为零，是不正确的。要注意二元函数的极限与二次极限的区别。故使用MATLAB计算极限，要在极限存在的前提下进行。5.4.2 符号导数diff函数用于对符号表达式求导数。该函数的一般调用格式为：diff(s)：没有指定变量和导数阶数，则系统按findsym函数指示的默认变量对符号表达式s求一阶导数。diff(s,v)：以v为自变量，对符号表达式s求一阶导数。diff(s,n)：按findsym函数指示的默认变量对符号表达式s求n阶导数，n为正整数。diff(s,v,n)：以v为自变量，对符号表达式s求n阶导数。例5-3 (1)设z=e2

16、xcos(3y)，求（化简过程中只写最后一个最简结果）。2. 设z=e2xcos(3y)，求,并求(2)y=ex(sinx+cosx)z=exp(2*x)*cos(3*y);f1=diff(diff(z,x),y) f1 =-6*exp(2*x)*sin(3*y)%或f1=diff(diff(z,y),x) 结果同上!ff2=subs(f1,x,1); ff2=subs(ff2,y,pi)ff2 = -1.6288e-014（2） clear;syms x y y=exp(x)*(sin(x)+cos(x); dy=diff(y,x)dy =exp(x)*(sin(x)+cos(x)+exp(

17、x)*(cos(x)-sin(x) simple(dy)ans =2*exp(x)*cos(x) simplify(dy)ans =2*exp(x)*cos(x)例5-3 求下列函数的导数。(同前面合并！)(2)， clear;syms x y s=diff(x2*exp(-y),x,2)s =2*exp(-y) t=diff(x2*exp(-y),x); t=diff(t,y); t=subs(t,x,1) t=subs(t,y,2)t = -0.2707例5-4. (1)求由方程y5+2y-x-3x7=0 所确定的隐函数的导数。(2).求函数z=ln(x2+y2)的全微分, 并化简（化简过

18、程中只写最后一个最简结果）。(3)求摆线，在t=所对应的点的切线斜率(4)求y=ln(sinx)的微分 syms x y; fxy=y5+2*y-x-3*x7;fx=diff(fxy,x);fy=diff(fxy,y); dv=-fx/fy dv = (1+21*x6)/(5*y4+2)syms x y dx dyz=log(x2+y2);dz=diff(z,x)*dx+diff(z,y)*dy;dz=simple(dz)dz =2*(x*dx+y*dy)/(x2+y2) syms a t; x=a*(t-sin(t);y=a*(1-cos(t); f1=diff(y,t);f2=diff(x

19、,t); f=f1/f2f =sin(t)/(1-cos(t) ff=inline(f); k=ff(pi)k = 6.1232e-017 syms x y dx dy y=log(sin(x); f1=diff(y,x); dy=f1*dxdy =cos(x)/sin(x)*dx5.4.3 符号积分符号积分由函数int来实现。该函数的一般调用格式为：int(s)：没有指定积分变量和积分阶数时，系统按findsym函数指示的默认变量对被积函数或符号表达式s求不定积分。int(s,v)：以v为自变量，对被积函数或符号表达式s求不定积分。int(s,v,a,b)：求定积分运算。a,b分别表示定积分

20、的下限和上限。该函数求被积函数在区间a,b上的定积分。a和b可以是两个具体的数，也可以是一个符号表达式，还可以是无穷(inf)。当函数f关于变量x在闭区间a,b上可积时，函数返回一个定积分结果。当a,b中有一个是inf时，函数返回一个广义积分。当a,b中有一个符号表达式时，函数返回一个符号函数。当系统求不出解析解，会自动求原点附近的一个近似解。例5-5 求下列积分。syms x f=(x2+1)/(x2-2*x+2)2;g=cos(x)/(sin(x)+cos(x);h=exp(-x2);I=int(f)J=int(g,0,pi/2)K=int(h,0,inf) 结果为： I =1/4*(2*

21、x-6)/(x2-2*x+2)+3/2*atan(x-1) J =1/4*pi K =1/2*pi(1/2)4.求不定积分： 并化简（化简过程中只写最后一个最简结果）。(4)f2=int(x*exp(a*x) f2 =1/a2*(a*x*exp(a*x)-exp(a*x) simple(f2)ans =exp(a*x)*(a*x-1)/a2例5-5 计算积分（无解析解）（3） (4) (5)重积分syms x;t3=int(3*sin(x2)/x,1,4)t3 =3/2*sinint(16)-3/2*sinint(1) %特殊函数sinint(a)表示积分int(sin(a*x)/x,0,1)

22、t3=vpa(t3,5) %用vpa求数值解t3 =1.0279t4=int(exp(-xsin(x),0,1) %求不出解析解Warning: Explicit integral could not be found.In D:toolboxsymbolicsymint.m at line 58t4 =int(exp(-xsin(x),x = 0 . 1) t4=vpa(t4,5) %也可用vpa求数值解t4 =.45491 syms x y iy=int(2*sqrt(1-x2),y,-sqrt(1-x2),sqrt(1-x2); t5=int(iy,x,-1,1)t5 =16/35.5

23、级数5.5.1 级数符号求和求无穷级数的和需要符号表达式求和函数symsum，其调用格式为：symsum(s,v,n,m)其中s表示一个级数的通项，是一个符号表达式。v是求和变量，v省略时使用系统的默认变量。n和m是求和的开始项和末项。例5-6求级数 的和S, 以及前十项的部分和S1 syms n x S=symsum(1/n2, 1, inf)S =1/6*pi2 S1=symsum(1/n2, 1, 10)S1 =1968329/1270080重要说明：当求函数项级数的和S2时，可用命令： syms n x S2=symsum(x/n2, n, 1, inf)S2 =1/6*x*pi2 S

24、2=symsum(x/n2, n, 1, 10)S2 =1968329/1270080*x两点说明：（1）注意观察S2与S1的细微区别！（2）当通项公式的Matlab表达式较长时，表达式要加上单引号后面的练习中会遇到此问题例5-7 求下列级数之和。3. 求幂级数的和函数(注意最后结果要化简!) (3) syms n x; % (1分)s1=symsum(1/(2*n+1)/(2*x+1)(2*n+1),n,0,inf); % (4分)s1=simple(s1) % (3分)s1 =1/2*log(x+1)/x)5.5.2 函数的泰勒级数MATLAB提供了taylor函数将函数展开为幂级数，其调

25、用格式为：taylor(f,v,n,a)该函数将函数f按变量v展开为泰勒级数，展开到第n项(即变量v的n-1次幂)为止，n的缺省值为6。v的缺省值与diff函数相同。参数a指定将函数f在自变量v=a处展开，a的缺省值是0。例5-8 求函数在指定点的泰勒级数展开式。4. 将函数在x=0点展开成x的幂级数（取前9项，实际结果为4项）(4) f4=log(1+x)/(1-x); % (2分)taylor(f4,x,0,9) % (6分)ans =2*x+2/3*x3+2/5*x5+2/7*x74. 将函数, 在x=-1点展开成x的幂级数(采用默认阶数) (4) taylor(1/x2,-1) % (

26、8分) ans = 3+2*x+3*(x+1)2+4*(x+1)3+5*(x+1)4+6*(x+1)55.6符号代数方程求解在MATLAB中，求解用符号表达式表示的代数方程可由函数solve实现，其调用格式为：solve(s)：求解符号表达式s的代数方程，求解变量为默认变量。solve(s,v)：求解符号表达式s的代数方程，求解变量为v。solve(s1,s2,sn,v1,v2,vn)：求解符号表达式s1,s2,sn组成的代数方程组，求解变量分别v1,v2,vn。注：线性方程可用linsolve函数当系统求不出解析解，会自动求原点附近的一个近似解。例5-9 解方程(1) 二次方程ax2+bx+

27、c=0;(2) 无解析解 x2-3x+ex=2;(3) 方程组，这里y, z是未知量。 syms a b c x;solve(a*x2+b*x+c,x)ans = 1/2/a*(-b+(b2-4*a*c)(1/2) 1/2/a*(-b-(b2-4*a*c)(1/2) solve(x2-3*x+exp(x)-2,x) %无解析解，得一近似解ans =1.8160 syms y z u v w;s=solve(u*y2+v*z+w,y+z+w,y,z)s = y: 2x1 sym z: 2x1 sym s.y,s.z %察看结构的具体内容，结果略，也可用下列方法更直接 y,z=solve(u*y2

28、+v*z+w,y+z+w,y,z)y = -1/2/u*(-2*u*w-v+(4*u*w*v+v2-4*u*w)(1/2)-w -1/2/u*(-2*u*w-v-(4*u*w*v+v2-4*u*w)(1/2)-wz = 1/2/u*(-2*u*w-v+(4*u*w*v+v2-4*u*w)(1/2) 1/2/u*(-2*u*w-v-(4*u*w*v+v2-4*u*w)(1/2)例5-9 解下列方程。*4. 问k为何值时，下面的方程组有非零解？请写出相应的MATLAB命令x1 -3x3=0 x1+2x2+kx3=0 2x1+kx2 - x3=0(4)syms k %(1分)A=1 0 -3;1 2 k;2 k -1; %(1分) D=det(A) %(2分) D =10-k2-3*k solve(D)ans = -5 2当k=-5或k=2时,det(A)=0,从而有非零解!。5.7 符号常微分方程求解在MATLAB中，用大写字母D表示导数。例如，Dy表示y，D2y表示y，Dy