MATLAB课件 第6章MATLAB符号计算.ppt

单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,第,6,章,MATLAB,符号计算,目录,在工程、应用数学和科学上经常要用到符号计算功能。,MATLAB,开发商,Mathwork,公司以,maple,的内核为符号计算的引擎，依赖,MAOLE,已有的库函数，开发了在,MATLAB,环境下实现符号计算的工具箱,Symbolic Math Toolbox(,符号数学工具箱,),。,6.1,符号计算基础,6.2,符号导数及其应用,6.3,符号积分,6.4,级数,6.5,代数方程的符号求解,6.6,常微分方程的符号求解,6.1,符号计算基础,目录,符号数学工具箱是操作和解决符号表达式的版本号数学工具的集合，它包含复合、简化、微分、积分以及求代数方程和微分方程的工具。,6.1.1,符号对象,1.,建立符号变量和符号常数,(1)sym,函数,sym,函数用来建立单个符号量，例如，,a=sym(a),建立符号变量,a,，此后，用户可以在表达式中使用变量,a,进行各种运算。,例,6.1,考察符号变量和数值变量的差别。,在,MATLAB,命令窗口，输入命令：,目录,a=,sym(a);b,=,sym(b);c,=,sym(c);d,=,sym(d,);,%,定义,4,个符号变量,w=10;x=5;y=-8;z=11;,%,定义,4,个数值变量,A=,a,b;c,d,%,建立符号矩阵,A,A=,a,b,c,d,det(B,)%,计算数值矩阵,B,的行列式,B=,w,x;y,z,%,建立数值矩阵,B,B=,10 5,-8 11,ans=,a*d-b*c,det(A,)%,计算符号矩阵,A,的行列式,ans=,150,例,6.2,比较符号常数与数值在代数运算时的差别。,在,MATLAB,命令窗口，输入命令：,pi1=sym(pi);k1=sym(8);k2=sym(2);k3=sym(3);%,定义符号变量,pi2=pi;r1=8;r2=2;r3=3;%,定义数值变量,sin(pi1/3)%,计算符号表达式值,sin(pi2/3)%,计算数值表达式值,sqrt(k1)%,计算符号表达式值,sqrt(r1)%,计算数值表达式值,sqrt(k3+sqrt(k2)%,计算符号表达式值,sqrt(r3+sqrt(r2)%,计算数值表达式值,目录,ans,=1/2*3(1/2),ans,=0.8660,ans,=2*2(1/2),ans,=2,.8284,ans,=(3+2(1/2)(1/2),ans,=2,.1010,(2)syms,函数,syms,函数的一般调用格式为：,syms,var1 var2,varn,函数定义符号变量,var1,var2,varn,等。用这种格式定义符号变量时不要在变量名上加字符分界符,(),，变量间用空格而,不要,用逗号分隔。,目录,2.,建立符号表达式,例,6.3,用两种方法建立符号表达式。,在,MATLAB,窗口，输入命令：,U=sym(3*x2+5*y+2*x*y+6)%,定义符号表达式,U,syms,x y;%,建立符号变量,x,、,y,V=3*x2+5*y+2*x*y+6%,定义符号表达式,V,2*U-V+6%,求符号表达式的值,U=,3*x2+5*y+2*x*y+6,V=,3*x2+5*y+2*x*y+6,ans,=,3*x2+5*y+2*x*y+12,例,6.4,计算,3,阶范得蒙矩阵行列式的值。设,A,是一个由符号变量,a,b,c,确定的范得蒙矩阵。,命令如下：,syms,a b c;,U=a,b,c;,A=1,1,1;U;U.2%,建立范得蒙符号矩阵,det(A,)%,计算,A,的行列式值,目录,A=,1,1,1,a,b,c,a2,b2,c2,ans,=,b*c2-c*b2-a*c2+a*b2+a2*c-a2*b,ans,=,-(-,c+b,)*(a-c)*(a-b),factor(ans,),例,6.5,建立,x,y,的一般二元函数。,在,MATLAB,命令窗口，输入命令：,syms,x y;,f=sym(f(x,y);,目录,6.1.2,基本的符号运算,1.,符号表达式运算,(1),符号表达式的四则运算,例,6.6,符号表达式的四则运算示例。,在,MATLAB,命令窗口，输入命令：,syms,x y z;,f=2*x+x2*x-5*x+x3%,符号表达式的结果为最简形式,f=2*x/(5*x)%,符号表达式的结果为最简形式,f=(x+y)*(x-y)%,符号表达式的结果不是,x2-y2,，而是,(x+y)*(x-y),目录,f=,-3*x+2*x3,f=2/5,f=,(,x+y,)*(,x-y,),(2),因式分解与展开,factor(S),对,S,分解因式，,S,是符号表达式或符号矩阵。,expand(S),对,S,进行展开，,S,是符号表达式或符号矩阵。,collect(S),对,S,合并同类项，,S,是符号表达式或符号矩阵。,collect(S,v),对,S,按变量,v,合并同类项，,S,是符号表达式或符号矩阵。,目录,例,6.7,对符号矩阵,A,的每个元素分解因式。,命令如下：,syms,a b x y;,A=2*a2*b3*x2-4*a*b4*x3+10*a*b6*x4,3*x*y-5*x2;4,a3-b3;,factor(A)%,对,A,的每个元素分解因式,目录,ans,=,2*a*b3*x2*(5*b3*x2-2*b*,x+a,),-x*(-3*y+5*x),4,(a-b)*(a2+b*a+b2),例,6.8,计算表达式,S,的值。,命令如下：,syms,x y;,s=(-7*x2-8*y2)*(-x2+3*y2);,expand(s)%,对,s,展开,collect(s,x,)%,对,s,按变量,x,合并同类项,(,无同类项,),factor(ans,)%,对,ans,分解因式,ans,=,7*x4-13*x2*y2-24*y4,ans,=,7*x4-13*x2*y2-24*y4,ans,=,(8*y2+7*x2)*(x2-3*y2),(3),表达式化简,MATLAB,提供的对符号表达式化简的函数有：,simplify(S),应用函数规则对,S,进行化简。,simple(S),调用,MATLAB,的其他函数对表达式进行综合化简，并显示化简过程。,例,6.9,化简,命令如下：,syms,x y;,s=(x2+y2)2+(x2-y2)2;,simple(s)%MATLAB,自动调用多种函数对,s,进行化简，并显示每步结果,目录,2.,符号矩阵运算,transpose(S),返回,S,矩阵的转置矩阵。,determ(S,),返回,S,矩阵的行列式值。,colspace(S,),返回,S,矩阵列空间的基。,Q,D=,eigensys(S,)Q,返回,S,矩阵的特征向量，,D,返回,S,矩阵的特征值。,目录,6.1.3,符号表达式中变量的确定,MATLAB,中的符号可以表示符号变量和符号常数。,findsym,可以帮助用户查找一个符号表达式中的的符号变量。该函数的调用格式为：,findsym(S,n,),函数返回符号表达式,S,中的,n,个符号变量，若没有指定,n,，则返回,S,中的全部符号变量。,在求函数的极限、导数和积分时，如果用户没有明确指定自变量，,MATLAB,将按缺省原则确定主变量并对其进行相应微积分运算。可用,findsym(S,1),查找系统的缺省变量，事实上，,MATLAB,按离字符,x,最近原则确定缺省变量。,目录,6.2,符号导数及其应用,6.2.1,函数的极限,limit,函数的调用格式为：,limit(f,x,a),limit,函数的另一种功能是求单边极限，其调用格式为：,limit(f,x,a,right),或,limit(f,x,a,left),目录,例,6.10,求极限,在,MATLAB,命令窗口，输入命令：,syms,a m x;,f=(x(1/m)-a(1/m)/(x-a);,limit(f,x,a)%,求极限,(1),f=(sin(a+x)-sin(a-x)/x;,limit(f)%,求极限,(2),limit(f,inf,)%,求,f,函数在,x,(,包括,+,和,-),处的极限,limit(f,x,inf,left,)%,求极限,(3),f=(,sqrt(x)-sqrt(a)-sqrt(x-a)/sqrt(x,*x-a*a);,limit(f,x,a,right)%,求极限,(4),目录,ans=,a(1/m)/a/m,ans=,2*cos(a),ans=,0,ans=,0,ans=,-1/2*2(1/2)/a(1/2),6.2.2,符号函数求导及其应用,MATLAB,中的求导的函数为：,diff(f,x,n),diff,函数求函数,f,对变量,x,的,n,阶导数。参数,x,的用法同求极限函数,limit,，可以缺省，缺省值与,limit,相同，,n,的缺省值是,1,。,目录,例,6.11,求函数的导数。,命令如下：,syms,a b t x y z;,f=sqrt(1+exp(x);,diff(f)%,求,(1),。未指定求导变量和阶数，按缺省规则处理,f=x*,cos(x,);,diff(f,x,2)%,求,(2),。求,f,对,x,的二阶导数,diff(f,x,3)%,求,(2),。求,f,对,x,的三阶导数,f1=a*cos(t);f2=b*sin(t);,diff(f2)/diff(f1)%,求,(3),。按参数方程求导公式求,y,对,x,的导数,(diff(f1)*diff(f2,2)-diff(f1,2)*diff(f2)/(diff(f1)3%,求,(3),。求,y,对,x,的二阶导数,f=x*exp(y)/y2;,diff(f,x)%,求,(4),。,z,对,x,的偏导数,diff(f,y)%,求,(4),。,z,对,y,的偏导数,f=x2+y2+z2-a2;,zx,=-diff(f,x)/diff(f,z)%,求,(5),。按隐函数求导公式求,z,对,x,的偏导数,zy,=-diff(f,y)/diff(f,z)%,求,(5),。按隐函数求导公式求,z,对,y,的偏导数,目录,例,6.12,在曲线,y=x3+3x-2,上哪一点的切线与直线,y=4x-1,平行。,命令如下：,x=sym(x);,y=x3+3*x-2;%,定义曲线函数,f=diff(y);%,对曲线求导数,g=f-4;,solve(g)%,求方程,f-4=0,的根，即求曲线何处的导数为,4,目录,ans=,1/3*3(1/2),-1/3*3(1/2),6.3,符号积分,6.3.1,不定积分,在,MATLAB,中，求不定积分的函数是,int,，其调用格式为：,int(f,x,),int,函数求函数,f,对变量,x,的不定积分。参数,x,可以缺省，缺省原则与,diff,函数相同。,目录,例,6.13,求不定积分。,命令如下：,x=sym(x);,f=(3-x2)3;,int(f,)%,求不定积分,(1),f=sqrt(x3+x4);,int(f,)%,求不定积分,(2),g=,simple(ans,)%,调用,simple,函数对结果化简,目录,ans=,27*x-1/7*x7+9/5*x5-9*x3,ans=,-1/48*(x3+x4)(1/2)*(-16*(x2+x)(3/2)+12*(x2+x)(1/2)*x+6*(x2+x)(1/2)-3*log(1/2+x+(x2+x)(1/2)/x/(x+1)*x)(1/2),g=,1/3*(x+1)*x)(1/2)*x2+1/12*(x+1)*x)(1/2)*x-1/8*(x+1)*x)(1/2)+1/16*log(1/2+x+(x+1)*x)(1/2),6.3.2,符号函数的定积分,定积分在实际工作中有广泛的应用。在,MATLAB,中，定积分的计算使用函数：,int(f,x,a,b,),例,6.14,求定积分。,命令如下：,x=sym(x);t=sym(t);,int(abs(1-x),1,2)%,求定积分,(1),f=1/(1+x2);,int(f,-inf,inf,)%,求定积分,(2),int(4*t*x,x,2,sin(t)%,求定积分,(3),f=x3/(x-1)100;,I=int(f,2,3)%,用符号积分的方法求定积分,(4),double(I)%,将上述符号结果转换为数值,目录,ans=,1/2,ans=,pi,ans=,2*t*(sin(t)2-4),I=,97893129180187301565519001875382615/1192978373971185320372138406360121344,ans=,0.08205775671718,例,6.15,求椭球的体积。,命令如下：,syms,a b c z;,f=pi*a*b*(c2-z2)/c2;,V=,int(f,z,-c,c,),V=,4/3*pi*a*b*c,目录,例,6.16,轴的长度为,10,米，若该轴的线性密度计算公式是,f(x)=6+0.3x,千克,/,米,(,其中,x,为距轴的端点距离,),，求轴的质量。,(1),符号函数积分。在,MATLAB,命令窗口，输入命令：,syms,x;,f=6+0.3*x;,m=int(f,0,10),(2),数值积分。,先建立一个函数文件,fx.m,：,function,fx,=,fx(x,),fx,=6+0.3*x;,再在,MATLAB,命令窗口，输入命令：,m=quad(fx,0,10,1e-6),目录,m=,75,例,6.17,求空间曲线,c,从点,(0,，,0,，,0),到点,(3,，,3,，,2),的长度。求曲线,c,的长度是曲线一型,命令如下：,syms,t;,x=3*t;y=3*t2;z=2*t3;,f=diff(x,y,z,t)%,求,x,y,z,对参数,t,的导数,g=,sqrt(f,*f)%,计算一型积分公式中的根式部分,l=int(g,t,0,1)%,计算曲线,c,的长度,目录,6.3.3,积分变换,(,不讲）,1.,傅立叶,(Fourier),变换,在,MATLAB,中，进行傅立叶变换的函数是：,fourier(fx,x,t,),求函数,f(x),的傅立叶像函数,F(t),。,ifourier(Fw,t,x,),求傅立叶像函数,F(t),的原函数,f(x),。,目录,例,6.18,求函数的傅立叶变换及其逆变换。,命令如下：,syms,x t;,y=abs(x);,Ft=,fourier(y,x,t,)%,求,y,的傅立叶变换,fx,=,ifourier(Ft,t,x,)%,求,Ft,的傅立叶逆变换,2.,拉普拉斯,(Laplace),变换,在,MATLAB,中，进行拉普拉斯变换的函数是：,laplace(fx,x,t,),求函数,f(x),的拉普拉斯像函数,F(t),。,ilaplace(Fw,t,x,),求拉普拉斯像函数,F(t),的原函数,f(x),。,目录,例,6.19,计算,y=x2,的拉普拉斯变换及其逆变换,.,命令如下：,x=sym(x);y=x2;,Ft=,laplace(y,x,t,)%,对函数,y,进行拉普拉斯变换,fx,=,ilaplace(Ft,t,x,)%,对函数,Ft,进行拉普拉斯逆变换,目录,3.Z,变换,对数列,f(n),进行,z,变换的,MATLAB,函数是：,ztrans(fn,n,z,),求,fn,的,Z,变换像函数,F(z),iztrans(Fz,z,n,),求,Fz,的,z,变换原函数,f(n),例,6.20,求数列,fn=e-n,的,Z,变换及其逆变换。,命令如下：,syms,n z,fn=exp(-n);,Fz,=,ztrans(fn,n,z,)%,求,fn,的,Z,变换,f=,iztrans(Fz,z,n,)%,求,Fz,的逆,Z,变换,目录,4.,积分变换的应用,例,6.21,用拉普拉斯方法解微分方程。,目录,6.4,级数,6.4.1,级数的符号求和,级数符号求和函数,symsum,，调用格式为：,symsum(a,n,n0,nn),例,6.22,求级数之和。,命令如下：,n=,sym(n,);,s1=symsum(1/n2,n,1,inf),%,求,s1,s2=symsum(-1)(n+1)/n,1,inf),%,求,s2,。未指定求和变量，缺省为,n,s3=,symsum(n,*xn,n,1,inf),%,求,s3,。此处的求和变量,n,不能省略。,s4=symsum(n2,1,100),%,求,s4,。计算有限级数的和,s1=,1/6*pi2,s2=,log(2),s3=,x/(x-1)2,s4=,338350,6.4.2,函数的泰勒级数,MATLAB,中提供了将函数展开为幂级数的函数,taylor,，其调用格式为：,taylor(f,v,n,a,),其中,f,是一个符号表达式，代表要展开的解析函数；,n,表示展开的阶数直到,n-1,阶，默认值展开到,5,阶为止；,v,是独立变量；,a,是要展开的点，缺省时表示,a=0.,输入参数,n,v,a,的顺序可以是任意的，系统会根据输入参数值的特征自动进行识别。,例,6.23,求函数在指定点的泰勒展开式。,命令如下：,x=sym(x);,f1=(1+x+x2)/(1-x+x2);,f2=sqrt(1-2*x+x3)-(1-3*x+x2)(1/3);,taylor(f1,x,5),%,求,(1),，展开到,x,的,4,次幂时应选择,n=5,taylor(f2,6)%,求,(2),ans,=,1+2*x+2*x2-2*x4,ans,=,1/6*x2+x3+119/72*x4+239/72*x5,例,6.24,将多项式表示成,x+1,的幂的多项式。,命令如下：,x=sym(x);,p=1+3*x+5*x2-2*x3;,f=taylor(p,x,-1,4),例,6.25,应用泰勒公式近似计算 。,命令如下：,x=sym(x);,f=(1-x)(1/12);%,定义函数，,4000(1/12)=2f(96/212),g=taylor(f,4)%,求,f,的泰勒展开式,g,，有,4000(1/12)2g(96/212),b=96/212;,a=1-b/12-11/288*b2-253/10368*b3%,计算,g(b),2*a%,求,4000(1/12),的结果,4000(1/12)%,用,MATLAB,的乘方运算直接计算,目录,6.4.3,函数的傅立叶级数,MATLAB 5.x,版中，尚未提供求函数傅立叶级数的内部函数。下面我们自己设计一个简化的求任意函数的傅立叶级数的函数文件。,function,mfourier,=,mfourier(f,n,),syms,x a b c;,mfourier,=int(f,-pi,pi)/2;%,计算,a0,for i=1:n,a(i)=,int(f,*,cos(i,*x),-pi,pi);,b(i)=,int(f,*sin(i*x),-pi,pi);,mfourier,=,mfourier+a(i,)*,cos(i,*x)+b(i)*sin(i*x);,end,return,调用该函数时，需给出被展开的符号函数,f,和展开项数,n,，不可缺省。,目录,例,6.26,在,-,，,区间展开函数为傅立叶级数。,命令如下：,x=sym(x);a=sym(a);,f=x;,mfourier(f,5)%,求,f(x)=x,的傅立叶级数的前,5,项,f=abs(x);,mfourier(f,5)%,求,f(x)=|x|,的傅立叶级数的前,5,项,syms,a;,f=,cos(a,*x);,mfourier(f,6)%,求,f(x)=,cos(ax,),的傅立叶级数的前,6,项,f=sin(a*x);,mfourier(f,4)%,求,f(x)=sin(ax),的傅立叶级数的前,4,项,目录,6.5,代数方程的符号求解,6.5.1,线性方程组的符号求解,MATLAB,中提供了一个求解线性代数方程组的函数,linsolve,，其调用格式为：,linsolve(A,b,),目录,例,6.27,求线性方程组,AX=b,的解。,解方程组,(1),的命令如下：,A=34,8,4;3,34,3;3,6,8;,b=4;6;2;,X=,linsolve(A,b,)%,调用,linsolve,函数求,(1),的解,Ab%,用另一种方法求,(1),的解,解方程组,(2),的命令如下：,syms,a11 a12 a13 a21 a22 a23 a31 a32 a33 b1 b2 b3;,A=a11,a12,a13;a21,a22,a23;a31,a32,a33;,b=b1;b2;b3;,X=,linsolve(A,b,)%,调用,linsolve,函数求,(2),的解,XX=Ab%,用左除运算求,(2),的解,目录,6.5.2,非线性方程组的符号求解,求解非线性方程组的函数是,solve,，调用格式为：,solve(eqn1,eqn2,eqnN,var1,var2,varN,),例,6.28,解方程。,命令如下：,x=solve(1/(x+2)+4*x/(x2-4)=1+2/(x-2),x)%,解方程,(1),f=sym(x-(x3-4*x-7)(1/3)=1);,x=solve(f)%,解方程,(2),x=solve(2*sin(3*x-pi/4)=1)%,解方程,(3),x=solve(x+x*exp(x)-10,x)%,解方程,(4),。仅标出方程的左端,目录,例,6.29】,求方程组的解。,命令如下：,x y=solve(1/x3+1/y3=28,1/x+1/y=4,x,y)%,解方程组,(1),x y=solve(x+y-98,x(1/3)+y(1/3)-2,x,y)%,解方程组,(2),Warning:Explicit solution could not be found.,In C:MATLABR11toolboxsymbolicsolve.m at line 136,x=,empty sym,y=,对方程组,(2)MATLAB,给出了无解的结论，显然错误，请看完全与其同构的方程组,(3),。输入命令如下：,u,v=solve(u3+v3-98,u+v-2,u,v)%,解方程组,(3),x v=solve(x2+y2-5,2*x2-3*x*y-2*y2)%,解方程组,(4),目录,6.6,常微分方程的符号求解,MATLAB,的符号运算工具箱中提供了功能强大的求解常微分方程的函数,dsolve,。该函数的调用格式为：,dsolve(eqn1,condition,var),该函数求解微分方程,eqn1,在初值条件,condition,下的特解。参数,var,描述方程中的自变量符号，省略时按缺省原则处理，若没有给出初值条件,condition,，则求方程的通解。,dsolve,在求微分方程组时的调用格式为：,dsolve(eqn1,eqn2,eqnN,condition1,conditionN,var1,varN,),函数求解微分方程组,eqn1,、,、,eqnN,在初值条件,conditoion1,、,、,conditionN,下的解，若不给出初值条件，则求方程组的通解，,var1,、,、,varN,给出求解变量。,目录,例,6.30,求微分方程的通解。,命令如下：,y=dsolve(Dy-(x2+y2)/x2/2,x)%,解,(1),。方程的右端为,0,时可以不写,y=,dsolve(Dy,*x2+2*x*y-exp(x),x)%,解,(2),y=dsolve(Dy-x/y/sqrt(1-x2),x)%,解,(3),目录,例,6.31,求微分方程的特解。,命令如下,：,y=,dsolve(Dy,=2*x*y2,y(0)=1,x)%,解,(1),y=dsolve(Dy-x2/(1+y2),y(2)=1,x)%,解,(2),目录,例,6.32,用微分方程的数值解法和符号解法解方程，并对结果进行比较。,在,MATLAB,命令窗口，输入命令：,y=dsolve(Dy+2*y/x-4*x,y(1)=2,x)%,用符号方法得到方程的解析解,为了求方程的数值解，需要按要求建立一个函数文件,fxyy.m,：,function f=,fxyy(x,y,),f=(4*x2-2*y)/x;%,只能是,y=f(x,y),的形式，当不是这种形式时，要变形。,return,输入命令：,t,w=ode45(fxyy,1,2,2);%,得到区间,1,，,2,中的数值解，以向量,t,、,w,存储。,为了对两种结果进行比较，在同一个坐标系中作出两种结果的图形。输入命令：,x=linspace(1,2,100);,y=x.2+1./x.2;%,为作图把符号解的结果离散化,plot(x,y,b.,t,w,r-);,目录,6.6.3,常微分方程组求解,例,6.33,求微分方程组的解。,命令如下：,x,y=,dsolve(Dx,=4*x-2*,y,Dy,=2*x-y,t)%,解方程组,(1),x,y=dsolve(D2x-y,D2y+x,t)%,解方程组,(2),目录,x=C1+C2*exp(3*t),y=1/2*C2*exp(3*t)+2*C1,x=-C1*exp(1/2*2(1/2)*t)*sin(1/2*2(1/2)*t)-C2*exp(-1/2*2(1/2)*t)*sin(1/2*2(1/2)*t)+C3*exp(1/2*2(1/2)*t)*cos(1/2*2(1/2)*t)+C4*exp(-1/2*2(1/2)*t)*cos(1/2*2(1/2)*t),y=-C1*exp(1/2*2(1/2)*t)*cos(1/2*2(1/2)*t)+C2*exp(-1/2*2(1/2)*t)*cos(1/2*2(1/2)*t)-C3*exp(1/2*2(1/2)*t)*sin(1/2*2(1/2)*t)+C4*exp(-1/2*2(1/2)*t)*sin(1/2*2(1/2)*t),