ch_符号计算

第 2 章符号计算符号计算：解算数学表达式、方程不是在离散化的数值点上进行，而是凭借一系列恒等式，数学定理，通过推理和演绎，获得解析结果.特点：一，相对于MATLAB的数值计算“引擎”和“函数库”而言，符号计算的“引擎”和“函数库”是独立的.二，在相当一些场合，符号计算解算问题的指令和过程，显得比数值计算更自然、更简明.三，大多数理工科的本科学生在学过高等数学和其他专业基础课以后，比较习惯符号计算的解题理念和模式.Symbolic Math Toolbox™The computational engine underlying the toolboxes is the kernel of Maple software, a system developed primarily at the University of Waterloo, Canada and, more recently, at the Eidgen?ssiche Technische Hochschule, Zürich, Switzerland. Maple software is marketed and supported by Waterloo Maple, Inc.运算引擎MuPADMuPAD作为MATLAB7.8的符号计算工具箱，是一具有人工智能的数学软件.方程式可以处理复数计算，完美的绘图功能，图型输入，输出,可以输入多个2-D函数或极坐标函数或3-D 函数,选择所要绘图参数，就可以完成图形，以及图形的动画制作也是非常方便.数值计算结果并不是MATLAB命令行窗口所得的类似代码形式，而是规范数学格式.并拥有一内建的程序语言，帮助文档以及文本操作，文本操作在一定程度上可以取代word.MathWorks自从2008年10开始，在Matlab的新版本(Matlab2008a，即7.6之后)中使用MuPAD 内核替换原来的Maple符号计算内核！2.1符号对象和符号表达式MATLAB依靠基本符号对象（包括数字、参数、变量）、运算符及一些预定义函数来构造和衍生符号表达式和符号方程.2.1.1符号对象的创建和衍生一生成符号对象的基本规则●任何基本符号对象都必须借助专门的符号函数指令sym或syms定义.●任何包含符号对象的表达式或方程，将继承符号对象的属性.二符号数字符号（类）数字的定义：sym('Num') 创建一个符号数字Numsc=sym('Num') 创建一个符号常数sc，该常数值准确等于Num说明：Num代表一个具体的数字Num必须处于（英文状态下的）单引号内，构成字符串（关于字符串参见附录A.1）. 【例2.1-1】符号（类）数字与数值（类）数字之间的差异.a=pi+sqrt(5) % 创建方式sa=sym('pi+sqrt(5)')Ca=class(a) % 类别判断Csa=class(sa)vpa(sa-a)a =5.3777sa =pi+sqrt(5)Ca =doubleCsa =symans =.138223758410852e-16三符号参数表达式e-ax sinbx中的a,b称为参数.定义格式：syms Para 定义符号参数ParaPara=sym('Para')syms Para Flag 定义具有Flag指定属性的符号参数ParaPara=sym('Para', 'Flag')syms Para1 Para2 ParaN 定义Para1 Para2 ParaN为符号参数syms Para1 Para2 ParaN Flag 定义Para1 Para2 ParaN为具有Flag指定属性的符号参数●符号参数名不要用处于“字母表中小写字母x及其两侧的英文字母”开头.●Flag表示参数属性，可具体取以下词条：positive 表示那些符号参数取正实数；real 表示那些符号参数限定为实时；unreal 表示那些符号参数为不限定的复数.syms x a bint(1/(x),a,b)Warning: Warning, unable to determine if 0 is between a andb; try to use assumptions or set _EnvAllSolutions to trueWarning: Explicit integral could not be found.> In sym.int at 58ans =int(1/x,x = a .. b)Var=sym('x');Upp=sym('a','real');Low=sym('b','real');Intergral=int(1/(x),a,b)Warning: Warning, unable to determine if 0 is between a andb; try to use assumptions or set option _EnvAllSolutions totrueWarning: Explicit integral could not be found.> In sym.int at 58Intergral =int(1/x,x = a .. b)Var=sym('x');Upp=sym('a','positive');Low=sym('b','positive');Intergral=int(1/(x),a,b)Intergral =-log(a)+log(b)四符号变量e-ax sinbx中的x称为变量，符号变量的定义同符号参数.确定自由符号变量的规则：●在专门指定变量名的符号运算中，解题一定围绕指定变量名进行.●自动识别符号变量时，字母的优先次序为x,y,w,z,v等.自动识别表达式中自由、独立的符号变量的指令：findsym(EXPR) 确认表达式EXPR 中所有自由符号变量findsym(EXPR, N) 确认表达式EXPR 中距离x 最近的N 个自由符号变量【例2.1-2】用符号计算研究方程02=++w vz uz 的解.（1）不指定变量情况syms u v w z % 定义符号参数/变量 Eq=u*z^2+v*z+w;result_1=solve(Eq)%findsym(Eq,1)result_1 = -u*z^2-v*z ans = w（2）指定变量情况result_2=solve(Eq,z) result_2 =-1/2*(v-(v^2-4*u*w)^(1/2))/u -1/2*(v+(v^2-4*u*w)^(1/2))/u【例2.1-3】对独立自由符号变量的自动辨认. （1）syms a b x X Y % 定义符号参数/变量k=sym('3'); % 符号常数z=sym('c*sqrt(delta)+y*sin(theta1)'); % 直接定义符号表达式 EXPR=a*z*X+(b*x^2+k)*Y;% 构成衍生符号表达式（2）findsym(EXPR) ans =X, Y, a, b, c, delta, theta1, x, y（3）findsym(EXPR,1) ans = x（4）findsym(EXPR,2),findsym(EXPR,9) ans = x,y ans =x,y,theta1,delta,c,b,a,X,Y【例2.1-4】findsym 确定自由变量是对整个矩阵进行的.syms a b t u v x yA=[a+b*x,sin(t)+u;x*exp(-t),log(y)+v] findsym(A,5) A =[ a+b*x, sin(t)+u] [ x*exp(-t), log(y)+v] ans =x,y,v,u,t2.1.2 符号计算中的算符● 与数值计算中的算符在形状、名称和使用方法上几乎完全相同. ● 仅注意：在符号对象的关系运算符中，只有算符“==”，“~=”比较结果为“真”时，用1表示; 否则用0表示.2.1.3符号计算中的函数指令2.1.4符号对象的识别为了函数指令与数据对象的适配，MATLAB提供了用于识别数据对象属性的指令：class(var) 给出变量var的数据类别（如double，sym等）isa(var, 'Obj') 若变量var是Obj代表的类型，给出1，表示“真”whos 给出所有MATLAB内存变量的属性【例2.1-5】数据对象及其识别指令的使用.（1）cleara=1;b=2;c=3;d=4; % 产生4个数值变量Mn=[a,b;c,d] % 利用已赋值变量构成数值矩阵Mc='[a,b;c,d]' % 字符串中的a,b,c,d与前面输入的数值变量无关Ms=sym(Mc) % Ms是一个符号矩阵，它与前面各变量无关Mn =1 23 4Mc =[a,b;c,d]Ms =[ a, b][ c, d]（2）SizeMn=size(Mn)SizeMc=size(Mc)SizeMs=size(Ms)SizeMn =2 2SizeMc =1 9SizeMs =2 2（3）CMn=class(Mn)CMc=class(Mc)CMs=class(Ms)CMn =doubleCMc =charCMs =sym（4）isa(Mn,'double')isa(Mc,'char')isa(Ms,'sym')ans =1ans =1ans =1（5）whos Mn Mc MsName Size Bytes Class AttributesMc 1x9 18 charMn 2x2 32 doubleMs 2x2 312 sym2.2符号数字及表达式的操作2.2.1数值数字与符号数字之间的转换一数值数字向符号数字的转换在符号运算中，“数值类数字”会自动地转换为符号数字.亦可借助sym函数：sym(Num,'r')或sym(Num)数值类数字Num的广义有理表达sym(Num,'d') 数值类数字Num的“十进制浮点”近似表达sym(Num,'e') 数值类数字Num的带eps误差的理性近似表达sym(Num,'f') 数值类数字Num的“十六进制浮点”近似表达a=pia1=sym(a,'r')a2=sym(a,'d')a3=sym(a,'e')a4=sym(a,'f')a =3.1416a1 =pia2 =3.1415926535897931159979634685442a3 =pi-198*eps/359a4 ='1.921fb54442d18'*2^(1)a=pi+0.1;a1=sym(a,'r')a2=sym(a,'d')a3=sym(a,'e')a4=sym(a,'f')a1 =7299417733396965*2^(-51)a2 =3.2415926535897932048158054385567a3 =7299417733396965*2^(-51)a4 ='1.9eec82110f9e5'*2^(1)1、Num只能是数值类数字2、sym('Num')和sym(Num)区别问题'Num'数字字符串，理论值，Num表示数字，近似（双精度值）3、在符号运算中，“数值数字”自动转换为符号运算二符号数字向双精度数字转换double(Num_sym) 把符号数字Num_sym转换为双精度数字2.2.2符号数字的任意精度计算digits 显示当前环境下符号数字“十进制浮点”表示的有效数字位数digits(n) 设定符号数字“十进制浮点”表示的有效数字位数（默认32位）xs=vpa(x) 据表达式x得到digits指定精度下的符号数字xsxs=vpa(x,n) 据表达式x得到n位有效数字的符号数字xs【例2.2-1】digits, vpa指令的使用.digitsp0=sym('(1+sqrt(5))/2')pr=sym((1+sqrt(5))/2)p d=sym((1+sqrt(5))/2,'d')e32r=vpa(abs(p0-pr) )e16=vpa(abs(p0-pd),16) %16位变精度运算计算误差e32d=vpa(abs(p0-pd))e16=vpa(abs(p0-pd),64)digits2.2.3符号表达式的基本操作collect（合并同类项）factor（进行因式分解）numden（提取公因式）等最常用：simple(EXPR) 把EXPR转换成最简形式【例2.2-2】简化32381261+++=xxxf.syms xf=(1/x^3+6/x^2+12/x+8)^(1/3); g1=simple(f)g2=simple(g1)有时需要多次使用simplesimple(f) %直接看化简过程simple(ans)pretty(f)2.2.4 表达式中的置换操作 一子表达式置换操作[RS,ssub]=subexpr(S,ssub) 运用符号变量ssub 置换子表达式，并重写S 为RS 【例2.2-3】对符号矩阵⎥⎦⎤⎢⎣⎡d c b a 进行特征向量分解. clear allsyms a b c d W[V,D]=eig([a b;c d]) % V ：特征向量阵 D ：特征值阵[RVD,W]=subexpr([V;D],W) %对矩阵元素中的公共子表达式进行置换表达只置换较长的表达式二 通用置换指令RES=subs(ES,old,new) 用new 置换ES 中的old 后产生RESRES=subs(ES,new) 用new 置换ES 中的自由变量后产生RES 【例2.2-4】用简单算例演示subs 的置换规则. （1）产生符号函数syms a x;f=a*sin(x)+5f =a*sin(x)+5（2）符号表达式置换f1=subs(f,'sin(x)',sym('y'))%<2>class(f1)（3）符号常数置换f2=subs(f,{a,x},{2,sym('pi/3')}) % <3> a 被双精度数字置换，x 被符号数字置换 class(f2)（4）双精度数值置换f3=subs(f,{a,x},{2,pi/3})%<4>class(f3)（5）数值数组置换之一f4=subs(subs(f,a,2),x,0:pi/6:pi)%<5>class(f4)（6）数值数组置换之二f5=subs(f,{a,x},{0:6,0:pi/6:pi}) %<6> class(f5)f6=subs(f,{a,x},{0:7,0:pi/6:pi})%<6>class(f5)2.3符号微积分2.3.1极限和导数的符号计算大学本科高等数学中的大多数微积分问题，都能用符号计算解决，手工笔算演绎的烦劳都可以由计算机完成.limit(f,v,a) 求极限 )(lim v f av →【例2.3-1】试求kxx x ⎪⎭⎫⎝⎛-∞→11lim .syms x kLim_f=limit ((1-1/x)^(k*x),x ,inf)diff(f,v,n) 求nn dv fd (n 缺省时，默认n=1)【例2.3-2】求⎥⎦⎤⎢⎣⎡=x x t t a f ln cos 3求dx df ， 22dt f d ，dtdx fd 2.syms a t x;f=[a,t^3;t*cos(x), log(x)]; df=diff (f) %f 对x 的导数dfdt2=diff(f,t,2) %f 对x 的二阶导数 dfdxdt=diff(diff(f,x),t) %二阶混合导数【例2.3-3】求⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=)sin()cos(),(2121212x x x e x x x x f 的Jacobian （雅可比）矩阵⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎣⎡∂∂∂∂∂∂∂∂∂∂∂∂231322122111x f x f x f x f x f x f . syms x1 x2;f=[x1*exp(x2);x2;cos(x1)*sin(x2)];v=[x1 x2];fjac=jacobian(f,v)【例2.3-4】x x f sin )(=，求)0(x f '，)(x f x '. （1）clearsyms xsyms d positivef_p=sin(x); % df_p=limit((subs(f_p,x,x+d)-f_p)/d,d,0) % <4>df_p0=limit((subs(f_p,x,d)-subs(f_p,x,0))/d,d,0) % <5>（2）f_n=sin(-x); df_n=limit((f_n-subs(f_n,x,x-d))/d,d,0) % <7> df_n0=limit((subs(f_n,x,0)-subs(f_n,x,-d))/d,d,0) %<8>（3）f=sin(abs(x)); dfdx=diff(f,x)% <10>dfdx =cos(abs(x))*abs(1,x) %abs(1,x)表示x 一阶导数（4）xn=-3/2*pi:pi/50:0;xp=0:pi/50:3/2*pi;xnp=[xn,xp(2:end)]; hold onplot(xnp,subs(f,x,xnp),'k','LineWidth',2.5) % <13> plot(xn,subs(df_n,x,xn),'--r','LineWidth',2.5) plot(xp,subs(df_p,x,xp),':r','LineWidth',2.5)legend(char(f),char(df_n),char(df_p),'Location','NorthEast')% <16> grid onxlabel('x') hold off2.3.2 序列/级数的符号求和symsum(f,v,a,b) 求f 在变量v 取遍[a, b]中所有整数时的和.a,b 缺省时默认求和区间[0, v-1].【例2.3-8】求∑-=103],[t t k t ，∑∞=⎥⎦⎤⎢⎣⎡--12)1(,)12(1k k k k . syms k t;f1=[t k^3];f2=[1/(2*k-1)^2,(-1)^k/k];s1=simple(symsum (f1)) % f1的自变量被确认为t s2=simple(symsum (f2,1,inf)) % f1的自变量被确认为k syms x y;f1=[y x^3];f2=[1/(2*y-1)^2,(-1)^y/y];s1=simple(symsum(f1)) % f1的自变量被确认为x2.3.3 符号积分int(f,v) 求f 对变量v 的不定积分 int(f,v,a,b) 求f 对变量v 的定积分 【例2.3-9】求dx x x x +⎰11.clear syms xf=sqrt((1+x)/x)/x s=int(f,x)s=simple(simple(s)) f =((x + 1)/x)^(1/2)/x s =- 2*(1/x + 1)^(1/2) - 2*atan((1/x + 1)^(1/2)*i)*i s =- 2*(1/x + 1)^(1/2) - 2*atan((1/x + 1)^(1/2)*i)*i【例2.3-10】求dx x x bx ax ⋅⎥⎥⎦⎤⎢⎢⎣⎡⎰sin 12. syms a b x;f=[a*x,b*x^2;1/x,sin(x)];disp('The integral of f is');pretty(int(f))【例2.3-11】求积分⎰⎰⎰++2122222)(x x yx xydzdydx z y x .syms x y zF2=int (int (int (x^2+y^2+z^2,z,sqrt(x*y),x^2*y),y,sqrt(x),x^2),x,1,2) VF2=vpa (F2) %积分结果用32位数字表示【例 2.3-12】求阿基米德（Archimedes ）螺线)0(,>⋅=a a r θ在0=θ到ϕ间的曲线长度函数，并求出πϕ2,1==a 时的曲线长度. （1）syms a r theta1 phi1 positive x=r*cos(theta1);x=subs(x,r,a*theta1); y=r*sin(theta1);y=subs(y,r,a*theta1);dLdth=sqrt(diff(x,theta1)^2+diff(y,theta1)^2); L=simple(int(dLdth,theta1,0,phi1))（2）L_2pi=subs(L,[a,phi1],sym('[1,2*pi]')) L_2pi_vpa=vpa(L_2pi)L_2pi =asinh(2*pi)/2 + pi*(4*pi^2 + 1)^(1/2) L_2pi_vpa =21.256294148209098800702512272566（3）L1=subs(L,a,sym('1'));ezplot(L1*cos(phi1),L1*sin(phi1),[0,2*pi]) grid on hold on x1=subs(x,a,sym('1')); y1=subs(y,a,sym('1'));h1=ezplot(x1,y1,[0,2*pi]);set(h1,'Color','r','LineWidth',5) title(' ')legend('螺线长度-幅角曲线','阿基米德螺线') hold off2.4 利用MAPLE 的深层符号计算资源【例2.8-1】求⎰xdt t 0 ln 1，)1 ,0(∈x .注意：被积函数在0=t 无定义，在-=1t 处为负无穷. （1）syms t, syms x positive ft=1/log(t);sx=int(ft,t,0,x) （2）x=0.5:0.1:0.9sx_n=-mfun('Ei',1,-log(x))2.4.1 发挥MAPLE 的计算潜力【例2.8-2】求递推方程f n f n f n ()()()=----3122的通解. （1）gs1=maple('rsolve(f(n)=-3*f(n-1)-2*f(n-2),f(k));')%<1>gs1 =(2*f(0)+f(1))*(-1)^k+(-f(0)-f(1))*(-2)^k（2）gs2=maple('rsolve','f(n)=-3*f(n-1)-2*f(n-2)','f(k)')%<2>gs2 =(2*f(0)+f(1))*(-1)^k+(-f(0)-f(1))*(-2)^k2.4.2 符号计算结果的数值化绘图【例2.9-4】在)1 ,0(区间，图示t ln 1和⎰x dt t0 ln 1（在例2.8-1中已计算过积分）. （1）syms t x positiveft=1/log(t);sx=int(ft,t,0,x) sx =PIECEWISE([NaN, 1 < x],[-Ei(1,-log(x)), otherwise])（2）x=0.1:0.05:0.9; sx_n=-mfun('Ei',1,-log(x));（3）图示被积函数和积分函数clfezplot('1/log(t)',[0.1,0.9]) grid on hold onplot(x,sx_n,'LineWidth',3) Char1='1/ln(t)';Char2='{\int_0^x} 1/ln(t) dt'; title([Char1,' and ',Char2])legend(Char1,Char2,'Location','SouthWest') 交互近似积分 syms xf=sin(x)+log(x); rsums(f,1,10)1． 说出以下三条指令产生的结果各属于哪种数据类型，是“双精度”对象，还是“符号”对象？3/7+0.1, sym(3/7+0.1), vpa(sym(3/7+0.1))〖答案〗2． 在不加专门指定的情况下，以下符号表达式中的哪一个变量被认为是独立自由变量.sym('sin(w*t)') , sym('a*exp(-X)' ) , sym('z*exp(j*theta)') 〖答案〗3． 求符号矩阵⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=333231232221131211a a a a a a a a a A 的行列式值和逆，所得结果应采用“子表达式置换”简洁化.〖答案〗 syms a11 a12 a13 a21 a22 a23 a31 a32 a33 ssubdet(A),B=A^-1,[RS,ssub]=subexpr(B,ssub)4． 对函数⎩⎨⎧<≥=000)(k k a k f k，当a 为正实数时，求∑∞=-0)(k k z k f .（实际上，这就是根据定义求Z 变换问题.）〖答案〗syms a positivesyms z kf=a^k;symsum(f*z^(-k),k,0,inf)5． 对于0>x ，求12011122+∞=∑⎪⎭⎫ ⎝⎛+-+k k x x k .（提示：理论结果为x ln ）〖答案〗 syms x positivesyms z kf=2/(2*k+1)*(((x-1)/(x+1))^(2*k+1));symsum(f,k,0,inf)ans =piecewise([abs(x - 1) < x + 1, 2*atanh((x - 1)/(x + 1))])6． （1）通过符号计算求t t y sin )(=的导数dtdy .（2）然后根据此结果，求-=0t dt dy 和2π=t dt dy. 〖答案〗syms ty=abs(sin(t));dydt=diff(y)d0=limit(dydt,t,0,'left')dpi_2=limit(dydt,pi/2) dydt =sign(sin(t))*cos(t)d0 =-1dpi_2 =d =abs(1,sin(t))*cos(t)d0_ =-1dpi_2 =7． 求出dx x e x sin 7.15⎰--ππ的具有64位有效数字的积分值.〖答案〗syms xy=int(exp(-abs(x))*abs(sin(x)),x,-5*pi,1.7*pi)vpa(y,64)matlab 2010by =1/exp(pi) + exp(2*pi) + exp(3*pi) + exp(4*pi) + exp(5*pi)/2 - 1/(8*exp((17*pi)/10)) + exp(pi) - 5^(1/2)/(8*exp((17*pi)/10)) + (2^(1/2)*(5 - 5^(1/2))^(1/2))/(8*exp((17*pi)/10)) + 1ans =3617514.6356470887071000183934655005542427350578351234317736807048． 计算二重积分⎰⎰+211222)(x dydx y x .〖答案〗 9． 在]2,0[π区间，画出dt tt x y x ⎰=0sin )(曲线，并计算)5.4(y . 〖答案〗。