第5讲 Matlab符号计算

第5讲MATLAB的符号计算,符号计算基础5.1符号对象一种新的数据类型符号对象、符号表达式5.2基本的符号运算表达式运算(四则、化简、因式分解)、矩阵运算5.3符号表达式中变量的确定符号导数及其应用5.4函数的极限5.5符号函数求导及其应用（注：其help见SymbolicMathToolbox）,5.6不定积分5.7符号函数的定积分5.8积分变换(*)5.9级数的符号求和；5.10函数的泰勒级数；5.11函数的傅立叶级数(*)；5.12线性方程组的符号求解；5.13非线性方程组的符号求解；5.14常微分方程的符号求解；5.15常微分方程组求解,符号计算基础,Matlab中数学表达式的三种存在形式及其运用：、由数值对象参与运算的数学表达式（数值型）：组成：(数值)常数、数值型变量、函数、数学运算符；要求：执行表达式之前必须对其中所有数值型变量赋值；运算结果：一个数值或数值数组例如：x=1/3y=exp(x)*sin(2*x)-1/5,、由符号对象参与运算的数学表达式(符号型)：组成：(符号)常数、符号型变量、符号型常量、函数、数学运算符；要求：执行表达式之前必须对其中所有符号型变量、符号型常量进行类型定义；运算结果：一个符号型数学表达式例如：symsx;y=exp(x)*sin(2*x)-1/5,、以字符串形式表达的数学表达式(字符串型)：组成：由数值型或符号型数学表达式两端加就转换为字符串型数学表达式；要求：(1)若调用eval函数，可将其转化为数值型数学表达式并执行，但调用eval之前必须对其中所有数值型变量赋值；(2)若调用sym函数，可将其转化为符号型数学表达式；,运算结果：(1)若调用eval函数，可获得一个数值或数值数组;(2)若调用sym函数，获得一个符号型数学表达式;例如：d=1/3c=eval(d2-1/5*sin(x)e=sym(d2-1/5*sin(x)class(e)diff(e),5.1符号对象1.建立符号对象（类型定义）(1)syms函数：syms函数的一般调用格式为：symsvar1var2varn定义符号对象var1,var2,varn等。用这种格式定义符号量时符号量间用空格而不要用逗号分隔。还可以用如下方式定义：syms(var1,var2,varn)例如：symsxyzsyms(x,y,z),(2)sym函数：sym函数的一般调用格式为：f=sym(arg)将数值、数值表达式或字符串表达式arg转换为相应的符号对象后赋给f，则f成为一个符号量，它的内容可以是符号常数、符号变量甚至符号表达式。例如：f=sym(x2-5)diff(f),例1考察符号变量和数值变量的差别。在MATLAB命令窗口，输入命令：symsabcd;%定义4个符号变量w=10;x=5;y=-8;z=11;%定义4个数值变量A=a,b;c,d%建立符号矩阵AB=w,x;y,z%建立数值矩阵Bdet(A)%计算符号矩阵A的行列式det(B)%计算数值矩阵B的行列式class(A)%判断变量A的类型class(B),例2比较符号常数与数值在代数运算时的差别。在MATLAB命令窗口，输入命令：pi1=sym(pi);k1=sym(8);k2=sym(2);k3=sym(3);%定义4个符号变量，其内容由4个符号常数指定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)%计算数值表达式值,2.建立符号表达式例3用两种方法建立符号表达式。在MATLAB窗口，输入命令：(1)先定义符号量，再直接写出符号表达式：symsxy;%建立符号变量x、yV=3*x2+5*y+2*x*y+6%定义符号表达式V(2)用sym函数把字符串表达式转化为符号表达式：U=sym(x3+5*y)%定义符号表达式UU-V%求符号表达式的值,例4常数与符号常数的差异a1=1/3,pi/7,sqrt(5),pi+sqrt(5)%a1是数值常数a2=sym(1/3,pi/7,sqrt(5),pi+sqrt(5)%最接近的有理表示a3=sym(1/3,pi/7,sqrt(5),pi+sqrt(5)%准确的符号数值表示a23=a2-a3,例5计算3阶范得蒙矩阵行列式的值。设A是一个由符号变量a,b,c确定的范得蒙矩阵。命令如下：symsabc;U=a,b,c;A=1,1,1;U;U.2%建立范得蒙符号矩阵det(A)%计算A的行列式值,5.2基本的符号运算1.符号表达式运算(1)符号表达式的四则运算例6符号表达式的四则运算示例：symsxyz;f=2*x+x2*x-5*x+x3%符号表达式的结果为最简形式f=(x+y)*(x-y)%符号表达式的结果不是x2-y2，而是(x+y)*(x-y)simple(f),(2)表达式化简MATLAB提供的对符号表达式化简的函数有：simplify(S)应用函数规则对S进行化简。simple(S)调用MATLAB的其他函数对表达式进行综合化简，并显示化简过程。,例7已知向量t=1:5,试构造符号函数命令如下：t=1:5;symsxy;y=0;fork=1:5s=1;forl=1:ks=s*(x-t(l);endy=y+t(k)*s;endsimple(y)%MATLAB自动调用多种函数对s进行化简，并显示每步结果,(3)因式分解与展开令S是符号表达式或符号矩阵。factor(S)对S分解因式。expand(S)对S进行展开。collect(S)对S合并同类项。collect(S,v)对S按变量v合并同类项。,例8对符号矩阵A的每个元素分解因式。命令如下：symsabxy;A=2*a2*b3*x2-4*a*b4*x3+10*a*b6*x4,3*x*y-5*x2;4,a3-b3;factor(A)%对A的每个元素分解因式,例9计算表达式S的值。命令如下：symsxy;s=(-7*x2-8*y2+1)*(-x2+3*y2);expand(s)%对s展开collect(s,x)%对s按变量x合并同类项collect(s,y)%对s按变量y合并同类项factor(ans)%对ans分解因式,2.符号矩阵运算S返回S矩阵的转置矩阵。rank(S)返回S矩阵的秩det(S)返回S矩阵的行列式值。inv(S)返回S矩阵的逆Q,D=eig(S)Q返回S矩阵的特征向量，D返回S矩阵的特征值。例10a=sym(magic(4);b=sym(ones(4)c=a+b/2crank(c)det(c)Q,D=eig(c),5.3符号表达式中变量的确定MATLAB中的符号可以表示符号变量和符号常量。findsym可以帮助用户查找一个符号表达式中的符号变量。该函数的调用格式为：findsym(S)找出表达式S中的所有符号变量findsym(S,n)从表达式S中找出最靠近x的n个符号变量,符号导数及其应用,5.4函数的极限limit函数的调用格式为：limit(f,x,a)x可缺省，缺省时使用findsym(f)找到的自变量；a也可缺省，缺省时a=0.limit函数的另一种功能是求单边极限，其调用格式为：limit(f,x,a,right)或limit(f,x,a,left),例13求极限,输入命令：symsamx;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)%求极限(3)f=(sqrt(x)-sqrt(a)+sqrt(x-a)/sqrt(x*x-a*a);limit(f,x,a,right)%求极限(4),5.5符号函数求导及其应用MATLAB中的求导的函数为：diff(f,x,n)求函数f对变量x的n阶导数。参数x的用法同求极限函数limit，可以缺省，缺省值与limit相同，n的缺省值是1。,例14求函数的导数：,symsabtxyz;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的三阶导数x=a*cos(t);y=b*sin(t);diff(y)/diff(x)%求(3),按参数方程求导公式求y对x的一阶导数(diff(x)*diff(y,2)-diff(x,2)*diff(y)/(diff(x)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的偏导数,例15在曲线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,符号积分,5.6不定积分在MATLAB中，求不定积分的函数是int，其调用格式为：int(f,x)int函数求函数f对变量x的不定积分。参数x可以缺省，缺省原则与diff函数相同。,例15求不定积分(1)(2)symsx;f1=(3-x2)3;f1_=int(f1)%求不定积分(1)f2=sqrt(x3+x4);f2_1=int(f2)%求不定积分(2)f2_2=simple(f2_1)%调用simple函数对结果化简f2_3=diff(f2_2)%验证f2_4=simple(f2_3),5.7符号函数的定积分定积分在实际工作中有广泛的应用。在MATLAB中，定积分的计算使用函数：int(f,x,a,b)求函数f对变量x的在区间a,b上的定积分。参数x可以缺省，缺省原则与diff函数相同。例16求定积分,命令如下：symsxt;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)formatlonggdouble(I)%将上述符号结果转换为数值,例17求椭球的体积命令如下：symsabcx;f=pi*c*b*(1-x2/a2);V=int(f,x,-a,a),例18轴的长度为10米，若该轴的线性密度计算公式是f(x)=6+0.3x千克/米(其中x为距轴的端点距离)，求轴的质量。(1)符号函数积分。在MATLAB命令窗口，输入命令：symsx;f=6+0.3*x;m=int(f,0,10)(2)数值积分。先建立一个函数文件fx.m：functionfx=fx(x)fx=6+0.3*x;再在MATLAB命令窗口，输入命令：m=quad(fx,0,10,1e-6),例19(I型曲线积分)求空间曲线c:从点(0，0，0)到点(3，3，2)的长度。命令如下：symst;x=3*t;y=3*t2;z=2*t3;f=diff(x,y,z,t)%求x,y,z对参数t的导数g=sqrt(f*f)%计算I型曲线积分公式中的根式部分l=int(g,t,0,1)%计算曲线c的长度,5.8积分变换(*)1.傅立叶(Fourier)变换:在MATLAB中，进行傅立叶变换的函数是：fourier(fx,x,t)求函数f(x)的傅立叶像函数F(t)。ifourier(Fw,t,x)求傅立叶像函数F(t)的原函数f(x)。例求函数y=|x|的傅立叶变换及其逆变换。命令如下：symsxt;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)。例计算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)例求数列fn=e-n的Z变换及其逆变换。命令如下：symsnzfn=exp(-n);Fz=ztrans(fn,n,z)%求fn的Z变换f=iztrans(Fz,z,n)%求Fz的逆Z变换,级数5.9级数的符号求和级数符号求和函数symsum，调用格式为：symsum(a,n,n0,nn)a为级数通项表达式，n为求和变量，n0、nn为求和的起点和终点。例20求级数之和,命令如下：symsnx;s1=symsum(1/n2,n,1,inf)%求(1),常数项级数s2=symsum(-1)(n+1)/n,1,inf)%求(2),交错级数。未指定求和变量，缺省为ns3=symsum(n*xn,n,1,inf)%求(3),函数项级数。此处的求和变量n不能省略s4=symsum(n2,1,100)%求(4),计算有限级数的和,5.10函数的泰勒级数MATLAB中提供了将函数展开为幂级数的函数taylor，其调用格式为：taylor(f,v,n,a)f为待展开函数表达式，v为自变量，n为展开阶数(正整数)，a则指定对f在v=a处进行泰勒展开。,命令如下：symsx;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=5taylor(f2,6)%求(2)。展开到x的5次幂,例21求函数的泰勒展开式,例22将多项式1+3x+5x2-2x3表示成x+1的幂的多项式。命令如下：symsx;p=1+3*x+5*x2-2*x3;f=taylor(p,x,4,-1),例23应用泰勒公式近似计算,命令如下：symsx;f=(1-x)(1/12);%定义函数，4000(1/12)=2f(96/212)g=taylor(f,4)%求f的4阶泰勒展开式gv=subs(g,x,96/212)%计算g(96/212),v1=simple(v)%化简计算结果2*vpa(v1,8)%求4000(1/12)的近似值,小数点后保%留8位有效数字，4000(1/12)2*g(96/212)class(ans)4000(1/12)%用MATLAB的乘方运算直接计算class(ans),5.11函数的傅立叶级数(*)MATLAB5.x版中，尚未提供求函数傅立叶级数的内部函数。下面我们自己设计一个简化的求任意函数的傅立叶级数的函数文件。functionmfourier=mfourier(f,n)symsxabc;mfourier=int(f,-pi,pi)/2/pi;%计算a0fori=1:na(i)=int(f*cos(i*x),-pi,pi)/pi;b(i)=int(f*sin(i*x),-pi,pi)/pi;mfourier=mfourier+a(i)*cos(i*x)+b(i)*sin(i*x);endreturn调用该函数时，需给出被展开的符号函数f和展开项数n，不可缺省。,例24在-，区间展开函数为傅立叶级数。命令如下：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项symsa;f=cos(a*x);mfourier(f,6)%求f(x)=cos(ax)的傅立叶级数的前6项f=sin(a*x);mfourier(f,4)%求f(x)=sin(ax)的傅立叶级数的前4项,代数方程的符号求解,5.12线性方程组的符号求解对于形如Ax=b的线性方程组，符号求解如下：A=sym(数值型系数矩阵)b=sym(值向量)x=Ab,例25求线性方程组AX=b的解。,解方程组(1)的命令如下：A=sym(34,8,4;3,34,3;3,6,8);b=sym(4;6;2);X=Ab解方程组(2)的命令如下：symsa11a12a13a21a22a23a31a32a33b1b2b3;A=a11,a12,a13;a21,a22,a23;a31,a32,a33;b=b1;b2;b3;X=Ab,5.13一般代数方程组的符号求解求解非线性方程组的函数是solve，调用格式：solve(eqn1,eqn2,eqnN,var1,var2,varN),例26解非线性方程,命令如下：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)，方程为符号型，自变量缺省symsx;x=solve(2*sin(3*x-pi/4)-1)%解方程(3)，这一用法不允许“”号x=solve(x+x*exp(x)-10,x)%解方程(4)，无符号解时，给出数值解,例27求非线性方程组的解:,命令如下：xy=solve(1/x3+1/y3=28,1/x+1/y=4,x,y)%解方程组(1)s=solve(x2+y2-5,2*x2-3*x*y-2*y2)%解方程组(2)，s是结构数组,5.14常微分方程的符号求解,MATLAB的符号运算工具箱中提供了功能强大的求解常微分方程的函数dsolve。该函数的调用格式为：dsolve(eqn1,condition,var)该函数求解微分方程eqn1在初值条件condition下的特解。参数var描述方程中的自变量符号，省略时按缺省原则处理，若没有