第6章 MATLAB符号计算 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进行各种运算。

a=sym('a');b=sym('b');c=sym('c');d=sym('d'); %定义4个符号变量 例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 B=[w,x;y,z] %建立数值矩阵B det(A) %计算符号矩阵A的行列式 det(B) %计算数值矩阵B的行列式

例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)) % 计算数值表达式值

(2)syms函数 syms函数的一般调用格式为： syms var1 var2 … varn 函数定义符号变量var1,var2,…,varn等。用这种格式定义符号变量时不要在变量名上加字符分界符(')，变量间用空格而不要用逗号分隔。

2. 建立符号表达式 例6.3用两种方法建立符号表达式。 在MATLAB窗口，输入命令： U=sym('3*x^2+5*y+2*x*y+6') %定义符号表达式U syms x y; %建立符号变量x、y V=3*x^2+5*y+2*x*y+6 %定义符号表达式V 2*U-V+6 %求符号表达式的值

例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的行列式值

例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+x^2*x-5*x+x^3 %符号表达式的结果为最简形式 f=2*x/(5*x) %符号表达式的结果为最简形式 f=(x+y)*(x-y) %符号表达式的结果不是x^2-y^2，而是(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*a^2*b^3*x^2-4*a*b^4*x^3+10*a*b^6*x^4,3*x*y-5*x^2;4,a^3-b^3]; factor(A) %对A的每个元素分解因式

例6.8 计算表达式S的值。 命令如下： syms x y; s=(-7*x^2-8*y^2)*(-x^2+3*y^2); expand(s) %对s展开 collect(s,x) %对s按变量x合并同类项(无同类项) factor(ans) % 对ans分解因式

MATLAB提供的对符号表达式化简的函数有： simplify(S) 应用函数规则对S进行化简。 (3)表达式化简 MATLAB提供的对符号表达式化简的函数有： simplify(S) 应用函数规则对S进行化简。 simple(S) 调用MATLAB的其他函数对表达式进行综合化简，并显示化简过程。 例6.9化简 命令如下： syms x y; s=(x^2+y^2)^2+(x^2-y^2)^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)

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)/y^2; diff(f,x) %求(4)。z对x的偏导数 diff(f,y) %求(4)。z对y的偏导数 f=x^2+y^2+z^2-a^2; 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=x^3+3*x-2; %定义曲线函数 f=diff(y); %对曲线求导数 g=f-4; solve(g) %求方程f-4=0的根，即求曲线何处的导数为4

6.3 符号积分 6.3.1不定积分 在MATLAB中，求不定积分的函数是int，其调用格式为： int(f,x) int函数求函数f对变量x的不定积分。参数x可以缺省，缺省原则与diff函数相同。

例6.13求不定积分。 命令如下： x=sym('x'); f=(3-x^2)^3; int(f) %求不定积分(1) f=sqrt(x^3+x^4); int(f) %求不定积分(2) g=simple(ans) %调用simple函数对结果化简

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+x^2); int(f,-inf,inf) %求定积分(2) int(4*t*x,x,2,sin(t)) %求定积分(3) f=x^3/(x-1)^100; I=int(f,2,3) %用符号积分的方法求定积分(4) double(I) %将上述符号结果转换为数值

例6.15求椭球的体积。 命令如下： syms a b c z; f=pi*a*b*(c^2-z^2)/c^2; 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)

例6.17求空间曲线c从点(0，0，0)到点(3，3，2)的长度。求曲线c的长度是曲线一型 命令如下： syms t; x=3*t;y=3*t^2;z=2*t^3; 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=x^2; 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/n^2,n,1,inf) %求s1 s2=symsum((-1)^(n+1)/n,1,inf) %求s2。未指定求和变量，缺省为n s3=symsum(n*x^n,n,1,inf) %求s3。此处的求和变量n不能省略。 s4=symsum(n^2,1,100) %求s4。计算有限级数的和

MATLAB中提供了将函数展开为幂级数的函数taylor，其调用格式为： 6.4.2 函数的泰勒级数 MATLAB中提供了将函数展开为幂级数的函数taylor，其调用格式为： taylor(f,v,n,a) 例6.23求函数在指定点的泰勒展开式。 命令如下： x=sym('x'); f1=(1+x+x^2)/(1-x+x^2); f2=sqrt(1-2*x+x^3)-(1-3*x+x^2)^(1/3); taylor(f1,x,5) %求(1)。展开到x的4次幂时应选择n=5 taylor(f2,6) %求(2)。

例6.24将多项式表示成x+1的幂的多项式。 命令如下： x=sym('x'); p=1+3*x+5*x^2-2*x^3; f=taylor(p,x,-1,4) 例6.25应用泰勒公式近似计算 。 f=(1-x)^(1/12); %定义函数，4000^(1/12)=2f(96/2^12) g=taylor(f,4) %求f的泰勒展开式g，有4000^(1/12)≈2g(96/2^12) b=96/2^12; a=1-b/12-11/288*b^2-253/10368*b^3 %计算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)的解 A\b %用另一种方法求(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=A\b %用左除运算求(2)的解

6.5.2 非线性方程组的符号求解 求解非线性方程组的函数是solve，调用格式为： solve('eqn1','eqn2',…,'eqnN','var1,var2,…,varN') 例6.28 解方程。 命令如下： x=solve('1/(x+2)+4*x/(x^2-4)=1+2/(x-2)','x') %解方程(1) f=sym('x-(x^3-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/x^3+1/y^3=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:\MATLABR11\toolbox\symbolic\solve.m at line 136 x = [ empty sym ] y = [] 对方程组(2)MATLAB给出了无解的结论，显然错误，请看完全与其同构的方程组(3)。输入命令如下： [u,v]=solve('u^3+v^3-98','u+v-2','u,v') %解方程组(3) [x v]=solve('x^2+y^2-5','2*x^2-3*x*y-2*y^2') %解方程组(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-(x^2+y^2)/x^2/2','x') %解(1)。方程的右端为0时可以不写 y=dsolve('Dy*x^2+2*x*y-exp(x)','x') %解(2) y=dsolve('Dy-x/y/sqrt(1-x^2)','x') %解(3)

例6.31 求微分方程的特解。 命令如下： y=dsolve('Dy=2*x*y^2','y(0)=1','x') %解(1) y=dsolve('Dy-x^2/(1+y^2)','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*x^2-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)