>> simple(ans) % one more time convert(exp)： exp(i*x) convert(tan)： (1-tan(1/2*x)^2)/(1+tan(1/2*x)^2)+2*i*tan(1/2*x)/(1+tan(1/2*x)^2) ans = exp(i*x) 22.7 可变精度算术运算 因为数值的精度受每次操作所保留的数位的限制，所以数值的任何运算都会引入舍入误差，重复的多次数值运算会造成累积误差。而对符号表达式的运算是非常准确的，因为它们不需要进行数值运算，所以无舍入误差。对符号运算结果用函数eval或numeric，仅在结果转换时会引入舍入误差。 MATLAB对数的处理完全依靠计算机的浮点算术运算，显然在内存中进行运算，又快又好，只是浮点运算受到所支持字长的限制，每次操作会引入舍入误差，所以不能产生精确的结果。MATLAB中各个算术运算的相对精度大约是16位。相反，Maple的符号处理能力可以实现任何数位的运算。当缺省的数位增加时，每次计算就需要附加时间和计算机内存。 Maple缺省为16位的精度。函数digits返回全局Digits参数的当前值。Maple缺省准确度可以由digits(n)来改变，其中n是所期望的准确度数位。用这种方法增加准确度的副作用是，每个Maple函数随后进行的计算都以新的准确度为准，增加了计算时间。结果的显示不会改变，只有所用的Maple函数的缺省准确度受到影响。 另外有一个函数，它可以任何精度实行单个计算，而使全局的Digits参数不变。即：可变精度的算术或函数vpa，它以缺省的精度或任何指定的精度对单个符号表达式进行计算，并以同样的精度来显示结果。 >> format long % let ' s see all the usual digits >> pi % how about  to numeric accuracy ans= 3.14159265358979 >> digits % display the default ' Digits ' value Digits=16 >> vpa(' pi ') % how about  to ' Digits ' value ans= 3.141592653589793 >> digits(18) % change the default to 18 digits >> vpa(' pi ') % how about  to ' Digits ' accuracy ans= 3.14159265358979324 >> vpa( ' pi ' ，20) % how about  to 20 digits ans= 3.1415926535897932385 >> vpa( ' pi ' ，50) % how about  to 50 digits ans= 3.1415926535897932384626433832795028841971693993751 >> vpa( ' 2^(1/3) ' ，200) % the cube root of 2 to 200 digits ans= 1.2599210498948731647672106072782283505702514647015079800819751121552996765 1395948372939656243625509415431025603561566525939900240406137372284591103042 693552469606426166250009774745265654803068671854055 将函数vpa作用于符号矩阵，对它的每一个元素进行计算也同样达到所指定的位数。 >> A=sym( ' [1/4，log(sqrt(2))；exp(1)，3/7] ' ) A= [ 1/4 ， log(sqrt(2))] [exp(1) ， 3/7)] >> vpa(A，20) % evaluate to 20 digits ans= [ .2500000000000000000， .34657359027997265471] [2.7182818284590452354， .42857142857142857143] 22.8 方程求解 用MATLAB所具有的符号工具可以求解符号方程。有一些工具已经在前面介绍过，更多将在本节予以检验。 求解单个代数方程 我们在前面已经看到，MATLAB具有求解符号表达式的工具。如果表达式不是一个方程式(不含等号)，则在求解之前函数solve将表达式置成等于0。 >> solve( ' a*x^2+b*x+c ' ) % solve for the roots of the quadratic eqution ans= [1/2/a*(-b+(b^2-4*a*c)^1/2)] [1/2/a*(-b-(b^2-4*a*c)^1/2)] 结果是符号向量，其元素是方程的2个解。如果想对非缺省x变量求解，solve必须指定变量。 >> solve( ' a*x^2+b*x+c ' ， ' b ' ) % solve for b ans= -(a*x^2+c)/x 带有等号的符号方程也可以求解。 >> f=solve( ' cos(x)=sin(x) ' ) % solve for x f= 1/4*pi >> t=solve( ' tan(2*x)=sin(x) ' ) t= [ 0] [acos(1/2+1/2*3^(1/2))] [acos(1/2=1/2*3^(1/2))] 并得到数值解。 >> numeric(f) ans= 0.7854 >> numeric(t) ans= 0 0 + 0.8314i 1.9455 注意在求解周期函数方程时，有无穷多的解。在这种情况下，solve对解的搜索范围限制在接近于零的有限范围，并返回非唯一的解的子集。 如果不能求得符号解，就计算可变精度解。 >> x=solve( ' exp(x)=tan(x) ' ) x= 1.306326940423079 代数方程组求解 可以同时求解若干代数方程，语句solve(s1，s2，.....，sn)对缺省变量求解n个方程，语句solve(s1，s2，...，sn，' v1，v2，...，vn ')对n个' v1，v2，...vn '的未知数求解n个方程。 如何处理中小学典型的代数问题？ 黛安娜(Diane)想去看电影，她从小猪存钱罐倒出硬币并清点，她发现：  10美分的硬币数加上5美分的硬币总数的一半等于25美分的硬币数。  1美分的硬币数比5美分、10美分以及25美分的硬币总数多10。  25美分和10美分的硬币总数等于1美分的硬币数加上1/4的5美分的硬币数  25美分的硬币数和1美分的硬币数比5美分的硬币数加上8倍的10美分的硬币数多1。 如果电影票价为3.00美元，爆米花为1.00美元，糖棒为50美分，她有足够的钱去买这三样东西？ 首先，根据以上给出的信息列出一组线性方程，假如p，n，d和q分别表示1美分，5美分，10美分，和25美分的硬币数  然后，建立MATLAB符号方程并对变量求解。 >> eq1= ' d+(n+p)/2=q ' ； >> eq2= ' p=n+d+q-10 ' ； >> eq3= ' q+d=p+n/4 ' ； >> eq4= ' q+p=n+8*d-1 ' ； >>[pennies，nickles，dimes，quarters]=solve(equ1，equ2，equ3，equ4，' p，n，d，q ' ) pennies= 16 nickles= 8 dimes= 3 quarters= 15 所以，黛安娜有16枚1美分的硬币，8枚5美分的硬币，3枚10美分的硬币，15枚25美分的硬币，这就意味着 >> money=.01*16+.05*8+.10*3+.25*15 money= 4.6100 她就有足够的钱去买电影票，爆米花和糖棒并剩余11美分。 单个微分方程 常微分方程有时很难求解，MATLAB提供了功能强大的工具，可以帮助求解微分方程。函数dsovle计算常微分方程的符号解。因为我们要求解微分方程，就需要用一种方法将微分包含在表达式中。所以，dsovle句法与大多数其它函数有一些不同，用字母D来表示求微分，D2，D3等等表示重复求微分，并以此来设定方程。任何D后所跟的字母为因变量。方程=0用符号表达式D2y=0来表示。独立变量可以指定或由symvar规则选定为缺省。例如，一阶方程dy/dx=1+y2的通解为： >> dsolve( ' Dy=1+y^2 ' ) % find the general solution ans= -tan(-x+C1) 其中，C1是积分常数。求解初值y(0)=1的同一个方程就可产生： >> dsolve(' Dy=1+y^2 '，' y(0)=1 ') % add an initial condition y= tan(x+1/4*pi) 独立变量可用如下形式指定： >> dsolve(' Dy=1+y^2 '，' y(0)=1 '，' v ') % find solution to dy/dv ans= tan(v+1/4*pi) 让我们举一个二阶微分方程的例子，该方程有两个初始条件： =cos(2x)-y (0)=0 y(0)=1 >> y=dsolve(' D2y=cos(2*x)-y '，' Dy(0)=0 '，' y(0)=1 ') y= -2/3*cos(x)^2+1/3+4/3*cos(x) >> y=simple(y) % y looks like it can be simplified y= -1/3*cos(2*x)+4/3*cos(x) 通常，要求解的微分方程含有一阶以上的项，并以下述的形式表示： -2-3y=0 通解为： >> y=solve( 'D2y-2Dy-3*y=0 ') y= C1*exp(-x)+C2*exp(3*x) 加上初始条件：y(0)=0和y(1)=1可得到： >> y=solve( ' D2y-2Dy-3*y=0 ' ， ' y(0)=0，y(1)=1 ' ) y= 1/(exp(-1)-exp(3))*exp(-x)-1/(exp(-1)-exp(3))*exp(3*x) >> y=simple(y) % this looks like a candidate for simplification y= -(exp(-x)-exp(3*x))/(exp(3)-exp(-1)) >> pretty(y) % pretty it up exp(-x)-exp(3 x) - --------------------- exp(3) -exp(-1) 现在来绘制感兴趣的区域内的结果。 >> ezplot(y，[-6 2])  图22.3 符号函数y=-(exp(-x)-exp(3*x))/(exp(3)-exp(-1)) (-6≤x≤2) 微分方程组 函数dsolve也可同时处理若干个微分方程式，下面有两个线性一阶方程。 =3f+4g =-4f+3g 通解为： >> [f，g]=dsolve( ' Df=3*f+4*g ' ， ' Dg=-4*f+3*g ' ) f= C1*exp(3*x)*sin(4*x)+C2*exp(3*x)*cos(4*x) g= -C2*exp(3*x)*sin(4*x)+C1*exp(3*x)*cos(4*x) 加上初始条件：f(0)=0和g(0)=1，我们可以得到： >> [f，g]=dsolve( ' Df=3*f+4*g ' ， ' Dg=-4*f+3*g ' ， ' f(0)=0，g(0)=1 ' ) f= exp(3*x)*sin(4*x) g= exp(3*x)*cos(4*x) 22.9 线性代数和矩阵 在本节中，我们将介绍符号矩阵和MATLAB提供的工具，它用线性代数求解问题。 符号矩阵 符号矩阵和向量是数组，其元素为符号表达式，可用函数sym来产生。 >> A=sym( ' [a， b， c； b， c， a； c， a， b] ' ) A= [ a， b， c] [ b， c， a] [ c， a， b] >> G=sym( ' [cos(t)， sin(t)；-sin(t)，cos(t)] ' ) G= [ cos(t)， sin(t)] [-sin(t)， cos(t)] 函数sym也可以扩展成定义各元素的公式。注意，只有在这种情况下，i，j分别表示行列的位置；且不影响i，j的缺省值(它代表)。下面的例子建立了3×3的矩阵，其元素依赖于行和列的位置。 >> S=sym(3，3，' (i+j)+(i-j+s) ') % create a matrix using a formula S= [ 2/s， 3/(-1+s)， 4/(-2+s)] [1/(1+s)， 4/s， s/(-1+s)] [4/(2+s)， 5/(1+s)， 6/s] >> S=sym(3， 3， ' m '， ' n '， ' (m-n)/(m-n-t) ') % use m and n in another formula S= [ 0， -1/(-1-t)， -2/(-2-t)] [1/(1-t)， 0， -1/(-1-t)] [2/(2-t)， 1/(1-t)， 0] 函数sym也可以把数值矩阵转换成符号形式 >> M=[1.1， 1.2， 1.3； 2.1， 2.2， 2.3； 3.1， 3.2， 3.3] % a numeric matrix M= 1.1000 1.2000 1.3000 2.1000 2.2000 2.3000 3.1000 3.2000 3.3000 >> S=sym(M) % convert to symbolic form S= [11/10， 6/5， 13/10] [21/10， 11/5， 23/10] [31/10， 16/5， 33/10] 如果数值矩阵的元素可以指定为小的整数之比，则函数sym将采用有理分式表示。如果元素是无理数，则在符号形式中sym将用符号浮点数表示元素。 >> E=[exp(1) sqrt(2)] E= 2.7183 1.4142 >> sym(E) ans= [3060513257434036*2^(-50)， 3184525836262886*2^(-51)] 用函数symsize可以得到符号矩阵的大小(行，列数)。函数返回数值或向量，而不是符号表达式。symsize的四种形式说明如下： >> S=sym(' [a，b，c；d，e，f] ') % create a symbolic matrix S= [a，b，c] [d，e，f] >> d=symsize(S) % return the size of S as the 2-element vector d d= 2 3 >> [m，n]=symsize(S) % return the number of rows in m， and column in n m= 2 n= 3 >> m=symsize(S，1) % return the number of rows m= 2 >> n=symsize(S，2) % return the number of columns n= 3 数值数组用N(m，n)形式来访问单个元素，但符号数组元素必须用函数如sym(S，m，n)来获取。若能用相同的句法当然好，但MATLAB符号表达式的表示不方便。在内部，符号数组表示成一个字符串数组；而S(m，n)返回单个字母。所以，符号数组的各个元素，必须由符号函数，如sym，来指定,而不是直接指定。 >> G=sym(' [ab，cd；ef，gh] ') % create a 2-by-2 symbolic matrix G= [ab， cd] [ef ， gh] >> G(1，2) % this is the second character of the first row of G ans= a >> r=sym(G，1，2) % this is the second expressinon in the first row of G r= cd 记住，在上例中的符号矩阵G，实际是以2×7字符数组存储在计算机中。第一行为' [ab， cd] '，所以第二个元素是' a '。 最后，sym可以用于改变符号数组的一个元素。 >> sym(G，2，2，' pq ') % change the (2，2) element in G from ' gh ' to ' pq ' ans= [ab， cd] [ef， pq] 代数运算 用函数symadd，symsub，symmul和symdiv，对符号矩阵可以执行许多通用的代数运算，用sympow可计算乘幂，用transpose计算符号矩阵的转置。 >> G=sym(' [cos(t)，sin(t0；-sin(t0，cos(t)] ') % create a symbolic matrix G= [ cos(t)， sin(t)] [ -sin(t)， cos(t)] >> symadd(G， ' t ' ) % add ' t ' to each element ans= [ cos(t)+t， sin(t)+t] [ -sin(t)+t， cos(t)+t] >> symmul(G，G) % multiply G by G； Synpow(G，2) does the same thing ans= [cos(t)^2-sint(t)^2， 2*cos(t)*sin(t)] [ -2*cos(t)* sin(t)， cos(t)^2-sin(t)^2] >> simple(G) % try to simplify ans= [cos(2*t)， sin(2*t)] [-sin(2*t)， cos(2*t)] 下面通过证明G的转置是它的逆来证明G是正交阵。 >> I=symmul(G，transpose(G)) % multiply G by its transpose I= [cos(t)^2+sin(t)^2， 0] [ 0， cos(t)^2+sin(t)^2] >> simplify(I) % there appears to be a trig identity here ans= [1， 0] [0， 1] 正如所期望的那样，这是单位阵。 线性代数运算 用函数inverse和determ，可计算符号矩阵的逆阵以及行列式。 >> H=sym(hilb(3)) % the symbolic form of the numeric 3-by-3 Hilbert matrix H= [ 1， 1/2， 1/3] [1/2， 1/3， 1/4] [1/3， 1/4， 1/5] >> determ(H) % dind the determinant of H ans= 1/2160 >> J=inverse(H) % find the inverse of H J= [ 9， -36， 30] [-36， 192， -180] [ 30， -180， 180] >> determ(J) % find the determinant of the inverse ans= 2160 用函数linsolve求解齐次线性方程；这是等价于基本的MATLAB的逆斜杠\算子的符号，linsolve(A，B)对X方阵求解矩阵方程A*X=B。回到以前的硬币问题： 黛安娜想去看电影，她从小猪存钱罐倒出硬币并清点，她发现：  10美分的硬币数加上5美分的硬币总数的一半等于25美分的硬币数。  1美分的硬币数比5美分、10美分以及25美分的硬币总数多10。  25美分和10美分的硬币总数等于1美分的硬币数加上1/4的5美分的硬币数  25美分的硬币数和1美分的硬币数比5美分的硬币数加上8倍的10美分的硬币数多1。。 象上次所做的那样，列出线性方程组，令p，n，d和q分别为1美分，5美分，10美分，和25美分的硬币数  重新以p，n，d，q的顺序排列表达式 p/2+n/2+d-q=0 p-n-d-q=-10 -p-n/4+d+q=0 p-n-8d+q=1 接下来，列出方程系数的符号数组。 >> A=sym(' [1/2，1/2，1，-1；1，-1，-1，-1；-1，-1/4，1，1；1，-1，-8，1] ') A= [1/2， 1/2， 1， -1] [ 1， -1， -1， -1] [ -1， -1/4， 1， 1] [ 1， -1， -8， 1] >> B=sym(' [0；-10；0；-1] ') % Create the symbolic vector B B= [ 0] [-10] [ 0] [ -1] >> X=linsolve(A，B) % solve the symbolic system A ' X=B for x x= [16] [ 8] [ 3] [15] 结果是相同的，黛安娜有16枚1美分的硬币，8枚5美分的硬币，3枚10美分的硬币，15枚25美分的硬币。 其他特性 symop将其参量串接起来，并计算所得到的表达式。 >> f=' cos(x) ' % create an expression f= cos(x) >> symop( ' atan( ' ，f， ' + ' ，a， ' ) ' ， ' ^2 ' ) ans= atan(cos(x)+a)^2 在用函数symop时，若将数组和标量混合应慎重。例如： >> M=sym( ' [a b； c d] ' ) M= [a， b] [c， d] >> symop(M， ' + ' ， ' t ' ) ans= [a+t， b] [ c， d+t] 是把t加到M的对角线上。 函数charpoly求解矩阵的特征多项式。 >> G=sym(' [1，1/2；1/3，1/4] ') % create a symbolic matrix G= [ 1， 1/2] [1/3， 1/4] >> charpoly(G) % find the characteristic polynomial of G ans= x ' 2-5/4*x+1/12 用函数eigensys可以求得符号矩阵的特征根和特征向量 >> F=sym(' [1/2，1/4；1/4，1/2] ') % create a symbolic matrix F= [1/2，1/4] [1/4，1/2] >> eigensys(F) % find the eigenvalues of F ans= [3/4] [1/4] >> [V，E]=eigensys(F) % find eigenvalues E and eigenvectors v V= [-1， 1] [ 1， 1] E= [1/4] [3/4] 矩阵的约当(Jordan)标准型是特征值的对角矩阵；转换矩阵的列是特征向量。对于给定的矩阵A， jordan(A)求出非奇异的矩阵V，使得inv(V)*A*V成为约当标准型，函数jordan有两种形式 >> jordan(F) % find the Jordan form of F ，above ans= [1/4， 0] [ 0， 3/4] >> [V，J]=jordan(F) % find the Jordan form and eigenvectors V= [ 1/2， 1/2] [-1/2， 1/2] J= [1/4， 0] [ 0， 3/4] 标准形和特征向V的列是某些F可能的特征向量。 因为F是非奇异的，F的零空间的基是空矩阵，而列空间的基是单位阵。 >> F=sym( ' [1/2，1/4；1/4，1/ 2] ' ) % recreate F F= [1/2， 1/4] [1/4， 1/2] >> nullspace(F) % the nullspace of F is the empty matrix ans= [] >> colspace(F) % find the column space of F ans= [1， 0] [0， 1] 用函数singvals可求解矩阵奇异值。 >> A=sym(magic(3)) % Generate a 3-by-3 matrix A= [8， 1， 6] [3， 5， 7] [4， 9， 2] >> singvals(A) % find the singular expressions ans= [15.00000000000000] [6.928203230275511] [3.464101615137752] 函数jacobian(w，v)可相对于v求解w的雅可比(Jacobia)值。其结果的第(i，j)项是df(i)/dv(j)。注意，当f是标量时，f的雅可比值是f的梯度。 >> jacobian( ' u*exp(v) ' ，sym( ' u，v ' )) ans= [exp(v，u*exp(v)] 22.10 小结 下列各表(表22.2-表22.8)综合了符号数学工具箱的特性： 表22.2 符号表达式的运算   numeric 符号到数值的转换  pretty 显示悦目的符号输出  subs 替代子表达式  sym 建立符号矩阵或表达式  symadd 符号加法  symdiv 符号除法  symmul 符号乘法  symop 符号运算  sympow 符号表达式的幂运算  symrat 有理近似  symsub 符号减法  symvar 求符号变量   表22.3 符号表达式的简化   collect 合并同类项  expand 展开  factor 因式  simple 求解最简形式  simplify 简化  symsum 和级数   表22.4 符号多项式   charpoly 特征多项式  horner 嵌套多项式表示  numden 分子或分母的提取  poly2sym 多项式向量到符号的转换  sym2poly 符号到多项式向量的转换   表22.5 符号微积分   diff 微分  int 积分  jordan 约当标准形  taylor 泰勒级数展开   表22.6 符号可变精度算术   digits 设置可变精度    vpa 可变精度计算     表22.7 求解符号方程   compose 函数的复合  dsolve 微分方程的求解  finverse 函数逆  linsolve 齐次线性方程组的求解  solve 代数方程的求解   表22.8 符号线性代数   charploy 特征多项式  determ 矩阵行列式的值  eigensys 特征值和特征向量  inverse 矩阵逆  jordan 约当标准形  linsolve 齐次线性方程组的解  transpose 矩阵的转置