Matlab学习指导第二章-符号运算解读课件

第二章符号计算11/27/20231符号计算功能matlab自产生起就在数值计算上功能卓著,深受各专业计算人员的欢迎.但由于在数学,物理等各种科研和工程应用中经常会遇到符号运算的问题.为此,公司于1993年购买了Maple软件的使用权,并在此基础上,开发了符号计算工具箱(SymbolicToolbox)MATLAB中实现符号计算功能的三种途径1)调用matlab自己开发的各种功能函数进行常用的符号运算.包括符号表达式与符号矩阵的基本操作,符号矩阵的运算,符号微积分,符号线性方程求解,符号微分方程求解等.2)为了特殊专业人员提供方便,matlab还保留mpa.m和maple.m两个函数与Maple接口．3)符号函数计算器功能．11/27/20232符号运算与数值运算的区别：符号运算中，解算数学表达式、方程时，不是在离散化的数值点上进行，而是凭借一系列恒等式和数学定理，通过推理和演绎，获得解析结果。这种计算建立在数值完全准确表达和推演严格解析的基础上，所得结果是完全准确的。符号运算----代数运算，公式推导数值运算---算术运算代值11/27/202332.1符号对象和符号表达式在matlab中,数值和数值变量用于数值的存储和各种数值计算.而符号常量,符号变量,符号函数,符号操作等则是用来形成符号表达式,严格按照代数,微积分等课程中的规则,公式进行运算,并尽可能给出解析表达式.2.1.1符号对象的生成和使用●

数值计算---变量先赋值,再使用.●

符号计算---先定义基本的符号对象(可以是常量,变量,表达式),然后用这些基本符号对象去构成新的表达式,再进行所需的符号运算．11/27/20234符号对象的创建和衍生1.生成符号对象的基本规则①任何基本符号对象（数字、参数、变量、表达式）都必须借助专门的符号函数指令sym或syms定义。②任何包含符号对象的表达式或方程，将继承符号对象的属性。即任何包含符号对象的表达式、方程也一定是符号对象。11/27/20235符号对象的定义：f=sym(arg)---把数字,字符串或表达式arg转为符号对象f=sym(argn,flagn)---把数值或数值表达式转换为flagn格式的符号对象f=sym(1/2),s=sym(1/3),y=sym(pi),a=sym('b*c')f=1/2,s=1/3,y=pi,a=b*cf=sym(1/2,'d'),s=sym(1/3,'d'),y=sym(pi,'d')f=sym(1/2,'r')=sym(1/2)f=.50000000000000000000000000000000s=.3333333333333333148296162562473911/27/20236syms(‘argv1’,‘argv2’,‘argvk’)---把字符argv1,argv2,argvk定义为基本符号对象symsargv1argv2argvk---上述格式的简洁形式,各符号对象间不得有逗号argv=sym(‘argv’,flagv)---按flagv指定的要求把字符串‘argv’定义为符号对象argvf=sym('pi'),s=sym('2*pi+sin(60*pi*180)+exp(2)')f=pi,s=2*pi+sin(60*pi*180)+exp(2)syms('a','b','c')a=a,b=b,c=csymsabca=a,b=b,c=c只能定义符号变量11/27/20237说明：•f=sym(argn,flagn)中的argn是数值或数值表达式时,flagn可选:‘d’---最接近的十进制浮点精确表示.‘r’---最接近的有理表示,缺省选项.有理是指用两个正整数p,q构成的p/q,p*pi/q,sqrt(p),2^q,10^q的形式之一．•argv=sym(‘argv’,flagv)中的‘argv’是字符时,flagv可取限定选项‘positive’---限定argv是“正，实”符号变量‘real’---限定argv是“实”符号变量．‘unreal’---argv是非实符号变量．11/27/202382符号数字的定义格式：sc=sym('num')

%sc为值为num的符号数字注意：i)单引号必须在英文状态下输入，构成字符串ii)num为一个具体的数字如：

sc=sym(‘2/3')sb=sym('pi+sqrt(5)')sc=2/3sb=pi+sqrt(5)11/27/202393符号参数定义格式：i)symspara

para=sym('para')

symsa;a=sym('a')ii)symsparaflagpara=sym('para','flag')

symsapositive;a=sym('a','positive')iii)symsabcsymsabcflagflag为参数属性：positive----参数取正实数real-----参数为实数unreal-----参数为限定的复数4符号变量表达式中的自变量无逗号11/27/202310推荐格式:◆定义符号变量:syms变量名symsabc◆定义符号常数:常数名=sym(‘常数值’)K=sym(‘2/3’);◆定义符号常数:表达式名=sym(‘表达式’);f=sym(‘a*sin(x)+b’);11/27/202311例2.1.1-1符号对象的生成a1=[1/3,pi/7,sqrt(5),sqrt(9),5.1,pi+sqrt(5)]<1>a2=sym([1/3,pi/7,sqrt(5),sqrt(9),5.1,pi+sqrt(5)])<2>a3=sym(‘[1/3,pi/7,sqrt(5),sqrt(9),5.1,pi+sqrt(5)]’)<3>a1=0.33330.44882.23613.00005.10005.3777a2=[1/3,pi/7,sqrt(5),3,51/10,q*2^(-50)]a3=[1/3,pi/7,sqrt(5),sqrt(9),5.1,pi+sqrt(5)]q=6054707603575008a23=a2-a3a23=[0,0,0,0,0,189209612611719/35184372088832-pi-5^(1/2)]11/27/202312例5.1.1-2把字符表达式转换为符号变量y=sym('2*sin(x)*cos(x)')y=simple(y)

y=sin(2*x)说明：1)<1>是数值常数，<2>是最接近的有理表示，<3>是绝对准确的符号数值表示2)<2><3>部分元素相同，是因为<2>中的那几个元素是有理表示的基本形式，所以也是绝对准确的．3)<3>指令产生的符号数值总是绝对准确的，因此建议：在产生符号常量时应优先使用这种输入方式，把数值放入单引号对中，数组元素间一定要采用逗号分隔在符号运算中,如果事先没有对表达式中的独立变量进行定义,那么系统将自动检查哪些字符是符号函数,哪些是符号变量,并且总是把在英文字母表中离x最近的字母认做独立符号变量11/27/202313例2.1.1-3用符号计算验证三角等式symsphi1phi2;y=simple(sin(phi1)*cos(phi2)-cos(phi1)*sin(phi2));y=sin(phi1-phi2)例2.1.1-4求矩阵的行列式,逆和特征值symsa11a12a21a22;A=[a11,a12;a21,a22];DA=det(A),IA=inv(A),EA=eig(A)A=[a11,a12][a21,a22]DA=a11*a22-a12*a2111/27/202314例2.1.1-4验证积分symsAttaow;yf=int(A*exp(-i*w*t),t,-tao/2,tao/2);yf=simple(yf)yf=2*A*sin(1/2*tao*w)/w练习:展开sin(4x)(expand)

11/27/2023152.1.2符号计算中的算符和基本函数由于新版matlab采用了重载技术,使得用来构成符号计算表达式的算符和基本函数,无论在形状,名称,还是使用方法上,都与数值计算中的算符和基本函数几乎完全相同,这给编程带来极大的方便.(2)关系运算符在符号对象的比较中,没有大于,大于等于,小于,小于等于的概念,而只有是否等于的概念.”==““~=“分别用来对算符两边的对象进行相等和不等的比较,返回为逻辑量(1)基本运算符算符”+”,”-”,”*”,”\”,“/”,“^”分别构成矩阵的加,减,乘,左除,右除,求幂运算.算符”.*”,“./”,“.\”,“.^”分别实现元素对元素的数组乘,除,求幂运算.算符”'”,“.'”分别实现矩阵的共轭转置,非共轭转置11/27/202316(4)指数,对数函数函数sqrt,exp,expm在两者中用法相同．符号计算中只有自然对数，而没有数值计算中的log2,log10(5)复数函数conj,imag,real,abs在两者中用法相同.但在符号计算中没有求相角的指令.(6)矩阵代数指令在符号计算中,matlab提供的常用矩阵代数指令有：diag,triu,tril,inv,det,rank,rref,null,colspace,expm,poly,eig,svd(3)三角函数,双曲线函数以及他们的反函数除atan2只能用于数值计算外,另外的在两种运算中使用方法相同.11/27/202317数值计算对象,符号计算对象,字符串是MATALB中最常用的数据对象.他们遵循各自不同的运算法则,但有时在外形上却十分相似.MATLAB提供了一些识别不同数据对象的指令,常用的有class,isa,whos例2.1.3-1数据对象及其识别指令的使用(1)生成三种不同类型的矩阵，给出不同的显示形式clear,a=1;b=2;c=3;d=4%产生四个数值变量Mn=[a,b,c,d]%利用已赋值变量构成数值矩阵Mc=‘[a,b,c,d]’%字符串中的a,b,c,d与前面输入的数值变量无关Ms=sym(Mc)%符号变量，与前面的各变量无关Mn=1234Mc=[a,b;c,d]Ms=[a,b][c,d]2.1.3对象类别的识别11/27/202318(2)三个矩阵的大小不同SizeMn=size(Mn),SizeMc=size(Mc),SizeMs=size(Ms)SizeMn=22SizeMc=19SizeMs=22(3)用class获得每种矩阵的类别CMn=class(Mn),CMc=class(Mc),CMs=class(Ms)CMn=doubleCMc=charCMs=sym(4)用isa判断每种矩阵的类别isa(Mn,‘double’),isa(Mc,'char'),isa(Ms,'sym')ans=1ans=1ans=111/27/202319(5)利用whos观察内存变量的类别和其他属性

NameSizeBytesClassMc1x918chararrayMn2x232doublearrayMs2x2312symobjecta=0:1:6;theta=sym('30*pi/180*a')alfa=sym('30*pi/180')*atheta=30*pi/180*aalfa=[0,1/6*pi,1/3*pi,1/2*pi,2/3*pi,5/6*pi,pi]a与数组a无关11/27/2023202.1.4符号表达式中自由变量的确定findsym(expr)---确认表达式expr中所有自由符号变量findsym(expr,n)---从expr中确认出靠x最近的n个独立自变量例2.1.4-1对独立自由符号变量的自动辨认(1)生成符号变量symsabxXY;k=sym('3');z=sym('c*sqrt(delta)+y*sin(theta)');expr=a*z*X+(b*x^2+k)*Y;(2)找出expr中的全部自由符号变量findsym(expr)ans=X,Y,a,b,c,delta,theta,x,y11/27/202321(3)从expr中确定1,2,3个自由变量findsym(expr,1),findsym(expr,2),findsym(expr,3)ans=xans=x,yans=x,y,thetafindsym确认的是表达式中的”自由”,”独立”的符号变量,由于k不是自由的,z不是独立的,所以该指令不将其作为自由变量。该指令把expr表达式中n个最靠近x的自由符号变量确认为“独立自由变量”,并且认为大写字母离x的距离总大于所有小写字母离x的距离．作用:在符号函数调用中,当没有明确指定所针对的变量时,系统自动确定所操作的自由变量11/27/2023222.2符号对象的操作与转换2.2.1符号表达式的操作collect合并同类项numden提取公因式expand展开指定项simplify恒等式简化factor因式分解pretty习惯方式显示horner转换为嵌套形式simple简化例2.2.1-1简化symx;f=(1/x^3+6/x^2+12/x+8)^(1/3);g1=simple(f)g2=sinple(g1)g1=(2*x+1)/xg2=2+1/x验证:f2=g2^3;expand(f2)11/27/202323例2.2.1-2因式分解9876542fx=sym('9876542');factor(fx)ans=(2)*(13)*(19)*(19993)

例2.2.1-3因式分解x12-1symsx;f=x^12-1;factor(f)例2.2.1-4pretty,horner的使用symsx;f=exp(-x);f1=taylor(f);pretty(f1)fh=horner(f1)ans=(x-1)*(x^2+x+1)*(x+1)*(1-x+x^2)*(x^2+1)*(x^4-x^2+1)11/27/202324f1=1-x+1/2*x^2-1/6*x^3+1/24*x^4-1/120*x^5

23451-x+1/2x-1/6x+1/24x-1/120x

fh=1+(-1+(1/2+(-1/6+(1/24-1/120*x)*x)*x)*x)*xexpand的使用symsxy;f=sin(x+y)u=expand(f)f=

sin(x+y)u=

sin(x)*cos(y)+cos(x)*sin(y)11/27/2023252.2.2置换操作作用:把复杂表达式中所含的多个相同子表达式用一个符号代替,使表达简洁[RS,ssub]=subexpr(S,ssub)---运用符号变量ssub置换子表达式,重写S为RS例2.2.2-1表达式置换symsabcdW;[V,D]=eig([a,b;c,d]);[RVD,W]=subexpr([V;D],W)说明:被置换的子表达式是机器自动寻找的,置换原则为,只有比较长的子表达式才被置换,比较短的子表达式,即使出现多次,也不被置换.1)子表达式置换操作指令格式：11/27/2023262)通用置换指令RES=subs(ES,old,new)---用new置换ES中的old后产生RESRES=subs(ES,new)---用new置换ES中的自由变量后产生RES例2.2.2-2subs置换规则示例(1)产生符号函数symsax;f=a*sin(x)+5;(2)符号变量置换f1=subs(f,'sin(x)',sym('y'))(3)符号常数置换f2=subs(f,{a,x},{2,sym(pi/3)})f1=a*y+5a=2x=pi/3,符号数字f2=3^(1/2)+5double(f2)代值11/27/202327(4)双精度数值转换(所有自由变量被双精度数值取代)自定义函数的赋值f3=subs(f,[a,x],[2,pi/3])f3=6.7321(5)数值数组置换一(取a=2,x=0:pi/6:pi)f4=subs(subs(f,a,2),x,0:pi/6:pi)f4=5.006.006.737.006.736.005.00(6)数值数组置换二(取a=0:6,x=0:pi/6:pi)f5=subs(f,{a,x},{0:6,0:pi/6:pi})f5=5.005.506.738.008.467.505.0011/27/202328f6=subs(f,[a,x],{0:6,0:pi/6:pi})f6=5.005.506.738.008.467.505.00比较：f7=subs(f,[a,x],[2,pi/6])f8=subs(f,[a,x],[2,sym(pi/6)])6.732050807568883^(1/2)+5●subs指令的属性取决于new的属性,当new中全部为数值数字时,所得结果为双精度数据,当new中有符号数字时,所得结果也为字符数字。●new可以是数组【说明】11/27/202329求n^0.5(n=1,2,….10)之值symsx;f=sqrt(x);res=subs(f,x,[1:1:10])?思考：如何求f=a*n^0.5+b(其中(a=2,3,4…11,n=1,2,3,…10,b=0.1,0.2,0.3,…1))的值11/27/202330例2.2.2-2猴子吃桃问题：猴子第一天摘下若干个桃子，当即吃了一半，还不瘾，又多吃了一个，第二天早上又将剩下的桃子吃掉一半，又多吃了一个。以后每天早上都吃了前一天剩下的一半零一个。到第10天早上想再吃时，见只剩下一个桃子了。求第一天共摘了多少。

symsr;f=2*(r+1);u=f;fork=1:1:8u=subs(f,r,u);endtotal=subs(u,r,1)u=512*r+1022total=1534a(1)=1534;%验算fork=2:1:10a(k)=a(k-1)/2-1;end11/27/2023312.2.3符号数值精度控制和任意精度的计算数值计算受计算机字长的限制,每次数值操作都带有截断误差,因此任何数值计算不管采用什么算法都将产生积累误差.在matlab中,每个算术操作结果的相对精度约为16位数字.但是符号计算的结果是绝对准确的,不包含任何计算误差.符号计算中与数值精度计算有关的指令有：double(x)---把符号常数转化为16位相对精度的浮点数值对象digits(n)---设置以后的数值计算以n位相对精度进行xs=vpa(x)---在digits指定精度下，给出x的数值型符号结果xsxs=vpa(x,n)---在n位相对精度下，给出x的数值型符号结果xs11/27/202332说明：●double指令运作所得结果一定是双精度数值对象●除了vpa(x,n)对特定符号对象指定具体精度外,所有vpa(x)的精度都受到其前面的digits指令控制,digits指令缺省精度为32位●x可以是符号对象或数值对象,但运行后所得结果xs一定是符号对象11/27/2023332.2.4符号对象与其他数据对象间的转换数值,符号,字符是matlab中三种不同的数据类型,matlab为每种数据类型提供了各自特定的生成指令和操作指令.为实现不同数据类型间的交互,matlab提供了一系列的转换指令．

数值型符号结果符号常数数值符号量(表达式)字符串(表达式)

ASCII码vpadoublevpasymdoublestr2num,str2double,sscanfint2str,num2str,mat2strcharsymcharabsdouble11/27/202334数值数字和符号数字之间的相互转换1)数值数字

符号数字sk=sym(num,'flag')flag:r----数值类数字的广义有理表达d---十进制浮点近似表达e---带误差的理性近似表达f---十六进制浮点近似表达PI=sym(pi,'r');PI=pi注意：在符号运算中，数值类数字会自动转换为符号数字symst;y=1/3*sin(t)+2/4y=1/3*sin(t)+1/2sk=sym('num')sk=sym(num)11/27/202335PI=sym(pi,'d');PI=3.14159265358979311599796346854422)符号数字数值数字double(num_sym)PI=sym('pi');double(PI)ans=3.141611/27/2023362.2.5符号表达式的复合函数和反函数1)复合函数运算●compose(f,g)---返回复合函数f(g(y)),这里f=f(x),g=g(y)●compose(f,g,z)---返回自变量为z的复合函数f(g(z)),这里2)反函数●g=finverse(f)---返回符号函数f的反函数g,满足g(f(x))=x●g=finverse(f,v)---返回自变量为v的符号函数f的反函数使得f=f(x),g=g(y)g(f(v))=v,当f包含不止一个符号变量时,使用这种格式.11/27/202337例2.2.5-1复合函数和反函数示例symsxy;f=1/x^3;g=tan(y);fg=compose(f,g)fin=finverse(fg)fg=1/tan(y)^3fin=atan(1/y^(1/3))11/27/2023382.3符号微积分1)符号极限limit(F,x,a)---计算F在xa时的极限limit(F,x,a,'left/right')---计算左/右极限例5.3-1极限示例symsx;f=3/(x^2+6);val1=limit(f,x,inf)val2=limit(f,x,0)val1=0val2=1/211/27/202339例已知f(x)=sin(|x|),求f'x(0),f'x(x)f(x)=sin(|x|)曲线图11/27/202340symsx;symsdxpositive;x0=sym('0');f_p=sin(x);df_p=limit((subs(f_p,x,x+dx)-f_p)/dx,dx,0)df_p0=limit((sin(x0+dx)-sin(x0))/dx,dx,0)数学分析：f(x)=sin(|x|)11/27/202341绘图要绘制原函数在x>0时的导函数曲线，应先将导函数离散化，采用通用置换指令进行离散化deltx=pi/200;xx=0:deltx:2*pi;y=subs(df_p,x,xx);plot(xx,y);11/27/2023422)符号序列求和对于数学上的求和可用matlab的求和指令解决.其指令为:s=symsum(f,v,a,b)---求通式f字指定变量v取遍[a,b]中所有整数时的和.说明:●f是矩阵时,对矩阵的元素逐个求和,但自变量定义在整个矩阵上●v缺省时f中的自变量由findsym自动辨认;b可以取有限整数,也可以取无穷大inf●a,b可同时缺省,此时默认求和的自变量区间为[0,v-1]11/27/202343例2.3-2求symsn;f=[n^2,1/n^2];s1=symsum(f(1),1,100)s2=symsum(f(2),1,inf)s1=338350s2=1/6*pi^211/27/202344例2.3-3求symskt;f1=[t,k^3];f2=[1/(2*k-1)^2,(-1)^k/k];s1=simple(symsum(f1));s2=simple(symsum(f2,1,inf))3)符号微分和Jacobian矩阵求导数，高阶导数，偏导数和多元向量函数的Jacobian矩阵是数学分析的重要内容，有机器实现求导的指令有：dfdvn=diff(f,v,n)--->求fjac=jacobian(f,v)--->求多元向量函数的Jacobian矩阵11/27/202345●f是矩阵时,求导对矩阵的元素逐个进行,但自变量定义在整个矩阵上.●v缺省时,自变量会自动由findsym确认,n缺省时,默认n=1●在数值计算中，指令diff是用来求差分的．●如f的Jacobian矩阵为说明：11/27/202346例2.3-4求symsatx;f=[a,t^3;t*cos(x),log(x)];df=diff(f),dfdt2=diff(f,t,2),dfdxdt=diff(diff(f,x),t)11/27/2023474)符号积分与数值积分相比,符号积分指令简单,适应性强,但可能占用机器时间很长.intf=int(f,v)--->f对v的不定积分intf=int(f,v,a,b)--->f对v的定积分说明：●当f是矩阵时，积分对矩阵的元素逐个进行●v缺省时，积分对findsym确认的变量进行●a,b分别是积分的上下限．例2.3-5求symsabx;f=[a*x,b*x^2;1/x,sin(x)];int(f)11/27/202348例2.3-6求积分symsxyz;f=x^2+y^2+z^2;s=int(int(int(f,z,sqrt(x*y),x^2*y),y,sqrt(x),x^2),x,1,2)第一重积分第二重积分11/27/202349ft=ilaplace(Fs,s,t)--->求频域函数Fs的Laplace反变换ftFZ=ztrans(fn,n,z)--->求时域系列fn的Z变换FZfn=iztrans(FZ,z,n)--->求频域系列FZ的Z变换fn5)符号积分变换傅立叶变换,拉普拉斯变换,Z变换在信号处理和系统动态特性研究中起着重要作用．Fw=fourier(ft,t,w)--->求时域函数ft的傅立叶变换Fwft=ifourier(Fw,w,t)--->求频域函数Fw的傅立叶反变换ftFs=laplace(ft,t,s)--->求时域函数ft的Laplace变换Fs11/27/202350symsAttaoxw;Fw1=int(A*exp(-i*w*t),t,-tao/2,tao/2);Fw1=simple(Fw1)ft=exp(x-t)*exp(-i*w*t);Fw2=int(ft,t,x,inf)%被积变量不能省略11/27/2023512.4符号代数方程求解1)线性方程组的符号解矩阵计算是求解线性方程组最简单有效的方法,在matlab中,不管数据对象是数值还是符号,实现矩阵运算的指令形式几乎完全相同.例2.4-1求:d+n/2+p/2=q,n+d+q-p=10,q+d-n/4=p,q+p-n-8d=1线性方程组的解A=sym([1,1/2,1/2,-1;1,1,-1,1;1,-1/4,-1,1;-8,-1,1,1]);b=sym([0;10;0;1])x=A\b11/27/2023522)一般代数方程组的解这里所讲的一般代数方程包括线性,非线性和超越方程等,求解指令为solve.S=solve('

eq1','

eq2',…,‘eqn’,‘v1’,‘v2’,…,’vn’)--->方程组指定变量的解(推荐格式)S=solve(exp1,exp2,…,expn,v1,v2,…,vn)--->方程组指定变量的解(可用格式)11/27/202353说明：●‘eq1’,‘eq2’,…,‘eqn’或是字符串表达的方程,或是字符串表达式.‘v1’,‘v2’,…,’vn’是字符串表达的求解变量名．●exp1,exp2,…,expn只能是符号表达式.v1,v2,…,vn是求解的符号变量．●如eq1,eq2,…,eqn是不含等号的表达式,则指令是对eq1=0,

eq2=0,…,eqn=0求解．●S是一个结构数组,如要显示求解结果,必须采用S.V1,S.V2,…,S.Vn的援引方式.11/27/202354例2.4-1求(x+2)x=2的解symsx;S=solve('(x+2)^x=2','x')解方程组:x+sin(y)=3/2x*sin(y)=1/2S=solve('x+sin(y)=3/2','x*sin(y)=1/2');S.x=[1][1/2]S.y=[1/6*pi][1/2*pi]11/27/2023552.5符号微分方程求解S=dsolve(‘eqn1’,‘eqn2’,…)--->符号微分方程组求解指令格式.说明：●输入宗量包括三部分内容:微分方程,初始条件,指定独立变量.其中微分方程是必不可少的.其余可以没有．●对于独立变量,若要指定独立变