《MATLAB课件第6章MATLAB符号计算.ppt》由会员分享，可在线阅读，更多相关《MATLAB课件第6章MATLAB符号计算.ppt（50页珍藏版）》请在三一办公上搜索。
1、第6章 MATLAB符号计算,目录,在工程、应用数学和科学上经常要用到符号计算功能。MATLAB开发商Mathwork公司以maple的内核为符号计算的引擎，依赖MAOLE已有的库函数，开发了在MATLAB环境下实现符号计算的工具箱Symbolic Math Toolbox(符号数学工具箱)。,第6章 MATLAB符号计算目录 在工程、应用数学和科学上,6.1 符号计算基础6.2 符号导数及其应用6.3 符号积分6.4 级数6.5 代数方程的符号求解6.6 常微分方程的符号求解,6.1 符号计算基础,6.1 符号计算基础,目录,符号数学工具箱是操作和解决符号表达式的版本号数学工具的集合，它包含
2、复合、简化、微分、积分以及求代数方程和微分方程的工具。,6.1 符号计算基础目录符号数学工具箱是操作和解决符号表达式,6.1.1 符号对象1. 建立符号变量和符号常数(1)sym函数 sym函数用来建立单个符号量，例如，a=sym(a)建立符号变量a，此后，用户可以在表达式中使用变量a进行各种运算。,6.1.1 符号对象,例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,A
3、= a, b c, d, det(B) %计算数值矩阵B的行列式,B=w,x;y,z %建立数值矩阵B,B = 10 5 -8 11,ans = a*d-b*c, det(A) %计算符号矩阵A的行列式,ans = 150,例6.1考察符号变量和数值变量的差别。目录a=sym(,例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) % 计算
4、数值表达式值sqrt(k1) % 计算符号表达式值sqrt(r1) % 计算数值表达式值sqrt(k3+sqrt(k2) % 计算符号表达式值sqrt(r3+sqrt(r2) % 计算数值表达式值,目录,ans = 1/2*3(1/2),ans = 0.8660,ans = 2*2(1/2),ans =2.8284,ans = (3+2(1/2)(1/2),ans = 2.1010,例6.2比较符号常数与数值在代数运算时的差别。目录ans =,(2)syms函数syms函数的一般调用格式为：syms var1 var2 varn 函数定义符号变量var1,var2,varn等。用这种格式定义符
5、号变量时不要在变量名上加字符分界符()，变量间用空格而不要用逗号分隔。,目录,(2)syms函数目录,2. 建立符号表达式例6.3用两种方法建立符号表达式。在MATLAB窗口，输入命令：U=sym(3*x2+5*y+2*x*y+6) %定义符号表达式Usyms x y; %建立符号变量x、yV=3*x2+5*y+2*x*y+6 %定义符号表达式V 2*U-V+6 %求符号表达式的值,U = 3*x2+5*y+2*x*y+6,V = 3*x2+5*y+2*x*y+6,ans = 3*x2+5*y+2*x*y+12,2. 建立符号表达式U =V =ans =,例6.4计算3阶范得蒙矩阵行列式的值。
6、设A是一个由符号变量a,b,c确定的范得蒙矩阵。命令如下：syms a b c;U=a,b,c;A=1,1,1;U;U.2 %建立范得蒙符号矩阵det(A) %计算A的行列式值,目录,A = 1, 1, 1 a, b, c a2, b2, c2,ans = b*c2-c*b2-a*c2+a*b2+a2*c-a2*b,ans = -(-c+b)*(a-c)*(a-b),factor(ans),例6.4计算3阶范得蒙矩阵行列式的值。设A是一个由符号变量a,例6.5建立x,y的一般二元函数。在MATLAB命令窗口，输入命令：syms x y;f=sym(f(x,y);,目录,例6.5建立x,y的一般
7、二元函数。目录,6.1.2 基本的符号运算1. 符号表达式运算(1)符号表达式的四则运算例6.6符号表达式的四则运算示例。在 MATLAB命令窗口，输入命令：syms x y z; f=2*x+x2*x-5*x+x3 %符号表达式的结果为最简形式f=2*x/(5*x) %符号表达式的结果为最简形式f=(x+y)*(x-y) %符号表达式的结果不是x2-y2，而是(x+y)*(x-y),目录,f = -3*x+2*x3,f = 2/5f = (x+y)*(x-y),6.1.2 基本的符号运算目录f =f = 2/5,(2)因式分解与展开factor(S) 对S分解因式，S是符号表达式或符号矩阵。
8、expand(S) 对S进行展开，S是符号表达式或符号矩阵。collect(S) 对S合并同类项，S是符号表达式或符号矩阵。collect(S,v) 对S按变量v合并同类项，S是符号表达式或符号矩阵。,目录,(2)因式分解与展开目录,例6.7 对符号矩阵A的每个元素分解因式。命令如下：syms a b x y;A=2*a2*b3*x2-4*a*b4*x3+10*a*b6*x4,3*x*y-5*x2;4,a3-b3;factor(A) %对A的每个元素分解因式,目录,ans = 2*a*b3*x2*(5*b3*x2-2*b*x+a), -x*(-3*y+5*x) 4, (a-b)*(a2+b*a
9、+b2),例6.7 对符号矩阵A的每个元素分解因式。目录ans =,例6.8 计算表达式S的值。命令如下：syms x y;s=(-7*x2-8*y2)*(-x2+3*y2);expand(s) %对s展开collect(s,x) %对s按变量x合并同类项(无同类项)factor(ans) % 对ans分解因式,ans = 7*x4-13*x2*y2-24*y4,ans = 7*x4-13*x2*y2-24*y4,ans = (8*y2+7*x2)*(x2-3*y2),例6.8 计算表达式S的值。ans =ans =ans =,(3)表达式化简MATLAB提供的对符号表达式化简的函数有：sim
10、plify(S) 应用函数规则对S进行化简。simple(S) 调用MATLAB的其他函数对表达式进行综合化简，并显示化简过程。例6.9化简命令如下：syms x y;s=(x2+y2)2+(x2-y2)2;simple(s) %MATLAB自动调用多种函数对s进行化简，并显示每步结果,目录,(3)表达式化简目录,2. 符号矩阵运算transpose(S) 返回S矩阵的转置矩阵。determ(S) 返回S矩阵的行列式值。colspace(S) 返回S矩阵列空间的基。Q,D=eigensys(S) Q返回S矩阵的特征向量，D返回S矩阵的特征值。,目录,2. 符号矩阵运算目录,6.1.3 符号表达
11、式中变量的确定MATLAB中的符号可以表示符号变量和符号常数。findsym可以帮助用户查找一个符号表达式中的的符号变量。该函数的调用格式为：findsym(S,n)函数返回符号表达式S中的n个符号变量，若没有指定n，则返回S中的全部符号变量。在求函数的极限、导数和积分时，如果用户没有明确指定自变量，MATLAB将按缺省原则确定主变量并对其进行相应微积分运算。可用findsym(S,1)查找系统的缺省变量，事实上，MATLAB按离字符x最近原则确定缺省变量。,目录,6.1.3 符号表达式中变量的确定目录,6.2 符号导数及其应用,6.2.1函数的极限limit函数的调用格式为：limit(f,
12、x,a)limit函数的另一种功能是求单边极限，其调用格式为：limit(f,x,a,right) 或 limit(f,x,a,left),目录,6.2 符号导数及其应用6.2.1函数的极限目录,例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
13、-a)/sqrt(x*x-a*a);limit(f,x,a,right) %求极限(4),目录,ans = a(1/m)/a/mans = 2*cos(a) ans =0ans = 0,ans = -1/2*2(1/2)/a(1/2),例6.10求极限目录ans =ans =,6.2.2 符号函数求导及其应用MATLAB中的求导的函数为：diff(f,x,n)diff函数求函数f对变量x的n阶导数。参数x的用法同求极限函数limit，可以缺省，缺省值与limit相同，n的缺省值是1。,目录,6.2.2 符号函数求导及其应用目录,例6.11求函数的导数。命令如下：syms a b t x y z
14、;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)/y2;diff(f,x) %求(4)。z对x的偏导数diff(f,y) %求(4)。z对
15、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的偏导数,目录,例6.11求函数的导数。目录,例6.12在曲线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,目录,ans = 1/3*3(1/2) -1/3*3(1/2),例6.12在曲线y=x
16、3+3x-2上哪一点的切线与直线y=4,6.3 符号积分,6.3.1不定积分在MATLAB中，求不定积分的函数是int，其调用格式为：int(f,x)int函数求函数f对变量x的不定积分。参数x可以缺省，缺省原则与diff函数相同。,目录,6.3 符号积分6.3.1不定积分目录,例6.13求不定积分。命令如下：x=sym(x);f=(3-x2)3;int(f) %求不定积分(1)f=sqrt(x3+x4);int(f) %求不定积分(2)g=simple(ans) %调用simple函数对结果化简,目录,例6.13求不定积分。目录,ans = 27*x-1/7*x7+9/5*x5-9*x3 a
17、ns = -1/48*(x3+x4)(1/2)*(-16*(x2+x)(3/2)+12*(x2+x)(1/2)*x+6*(x2+x)(1/2)-3*log(1/2+x+(x2+x)(1/2)/x/(x+1)*x)(1/2) g = 1/3*(x+1)*x)(1/2)*x2+1/12*(x+1)*x)(1/2)*x-1/8*(x+1)*x)(1/2)+1/16*log(1/2+x+(x+1)*x)(1/2),ans =,6.3.2 符号函数的定积分定积分在实际工作中有广泛的应用。在MATLAB中，定积分的计算使用函数：int(f,x,a,b)例6.14求定积分。命令如下：x=sym(x);t=s
18、ym(t);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)double(I) %将上述符号结果转换为数值,目录,6.3.2 符号函数的定积分目录,ans = 1/2 ans = pians = 2*t*(sin(t)2-4)I = 97893129180187301565519001875382615/1192978373971185320372138406360121344
19、ans = 0.08205775671718,ans =,例6.15求椭球的体积。命令如下：syms a b c z;f=pi*a*b*(c2-z2)/c2;V=int(f,z,-c,c)V =4/3*pi*a*b*c,目录,例6.15求椭球的体积。目录,例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;再在M
20、ATLAB命令窗口，输入命令：m=quad(fx,0,10,1e-6),目录,m = 75,例6.16轴的长度为10米，若该轴的线性密度计算公式是f(x,例6.17求空间曲线c从点(0，0，0)到点(3，3，2)的长度。求曲线c的长度是曲线一型命令如下：syms t;x=3*t;y=3*t2;z=2*t3;f=diff(x,y,z,t) %求x,y,z对参数t的导数g=sqrt(f*f) %计算一型积分公式中的根式部分l=int(g,t,0,1) %计算曲线c的长度,目录,例6.17求空间曲线c从点(0，0，0)到点(3，3，2)的,6.3.3 积分变换(不讲）1. 傅立叶(Fourier)变
21、换在MATLAB中，进行傅立叶变换的函数是：fourier(fx,x,t) 求函数f(x)的傅立叶像函数F(t)。ifourier(Fw,t,x) 求傅立叶像函数F(t)的原函数f(x)。,目录,6.3.3 积分变换(不讲）目录,例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
22、(Fw,t,x) 求拉普拉斯像函数F(t)的原函数f(x)。,目录,例6.18求函数的傅立叶变换及其逆变换。目录,例6.19计算y=x2的拉普拉斯变换及其逆变换.命令如下：x=sym(x);y=x2;Ft=laplace(y,x,t) %对函数y进行拉普拉斯变换fx=ilaplace(Ft,t,x) %对函数Ft进行拉普拉斯逆变换,目录,例6.19计算y=x2的拉普拉斯变换及其逆变换.目录,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
23、的Z变换及其逆变换。命令如下：syms n zfn=exp(-n);Fz=ztrans(fn,n,z) %求fn的Z变换f=iztrans(Fz,z,n) %求Fz的逆Z变换,目录,3. Z变换目录,4. 积分变换的应用例6.21用拉普拉斯方法解微分方程。,目录,4. 积分变换的应用目录,6.4 级数6.4.1 级数的符号求和级数符号求和函数symsum，调用格式为：symsum(a,n,n0,nn) 例6.22求级数之和。命令如下：n=sym(n);s1=symsum(1/n2,n,1,inf) %求s1s2=symsum(-1)(n+1)/n,1,inf) %求s2。未指定求和变量，缺省为
24、ns3=symsum(n*xn,n,1,inf) %求s3。此处的求和变量n不能省略。s4=symsum(n2,1,100) %求s4。计算有限级数的和,s1 = 1/6*pi2 s2 = log(2) s3 = x/(x-1)2 s4 = 338350,6.4 级数s1 =,6.4.2 函数的泰勒级数MATLAB中提供了将函数展开为幂级数的函数taylor，其调用格式为：taylor(f,v,n,a)其中f是一个符号表达式，代表要展开的解析函数；n 表示展开的阶数直到n-1 阶，默认值展开到5阶为止；v是独立变量；a是要展开的点，缺省时表示a=0.输入参数n,v,a的顺序可以是任意的，系统会
25、根据输入参数值的特征自动进行识别。例6.23求函数在指定点的泰勒展开式。命令如下：x=sym(x);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),ans = 1+2*x+2*x2-2*x4 ans = 1/6*x2+x3+119/72*x4+239/72*x5,6.4.2 函数的泰勒级数ans =,例6.24将多项式表示成x+1的幂的多项式。命令如下：x=sym(x);p=1+3*x+5*x2-2*x3;f=taylor(
26、p,x,-1,4)例6.25应用泰勒公式近似计算 。命令如下：x=sym(x);f=(1-x)(1/12); %定义函数，4000(1/12)=2f(96/212)g=taylor(f,4) %求f的泰勒展开式g，有4000(1/12)2g(96/212)b=96/212;a=1-b/12-11/288*b2-253/10368*b3 %计算g(b)2*a %求4000(1/12)的结果4000(1/12) %用MATLAB的乘方运算直接计算,目录,例6.24将多项式表示成x+1的幂的多项式。目录,6.4.3 函数的傅立叶级数MATLAB 5.x版中，尚未提供求函数傅立叶级数的内部函数。下面我
27、们自己设计一个简化的求任意函数的傅立叶级数的函数文件。function mfourier=mfourier(f,n)syms x a b c;mfourier=int(f,-pi,pi)/2; %计算a0for 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);endreturn调用该函数时，需给出被展开的符号函数f和展开项数n，不可缺省。,目录,6.4.3 函数的傅立叶级数目录,例6.26在-，区间展开函数为傅立叶级数。命令
28、如下：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.26在-，区间展开函数为傅立叶级数。目录,6.5代数方程的符号求解,6.5.1线性方程组的符号求解MATLAB中提供了一个求解线性代数方程组的函数linsolve，其调用格式为
29、：linsolve(A,b),目录,6.5代数方程的符号求解6.5.1线性方程组的符号求解目录,例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)的解Ab %用另一种方法求(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函
30、数求(2)的解XX=Ab %用左除运算求(2)的解,目录,例6.27求线性方程组AX=b的解。目录,6.5.2 非线性方程组的符号求解求解非线性方程组的函数是solve，调用格式为：solve(eqn1,eqn2,eqnN,var1,var2,varN)例6.28 解方程。命令如下：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)x=solve(2*sin(3*x-pi/4)=1) %解方程(3)x=solve(x+x*exp(x)-10,x) %解方程(4)
31、。仅标出方程的左端,目录,6.5.2 非线性方程组的符号求解目录,例6.29】求方程组的解。命令如下：x y=solve(1/x3+1/y3=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:MATLABR11toolboxsymbolicsolve.m at line 136x = empty sym y = 对方程组(2)MATLAB给出了无解的结论，显然错误，请看完全与其同构的方程组(3)。输
32、入命令如下：u,v=solve(u3+v3-98,u+v-2,u,v) %解方程组(3)x v=solve(x2+y2-5,2*x2-3*x*y-2*y2) %解方程组(4),目录,例6.29】求方程组的解。目录,6.6常微分方程的符号求解,MATLAB的符号运算工具箱中提供了功能强大的求解常微分方程的函数dsolve。该函数的调用格式为：dsolve(eqn1,condition,var)该函数求解微分方程eqn1在初值条件condition下的特解。参数var描述方程中的自变量符号，省略时按缺省原则处理，若没有给出初值条件condition，则求方程的通解。dsolve在求微分方程组时的调
33、用格式为：dsolve(eqn1,eqn2,eqnN,condition1,conditionN,var1,varN) 函数求解微分方程组eqn1、eqnN在初值条件conditoion1、conditionN下的解，若不给出初值条件，则求方程组的通解，var1、varN给出求解变量。,目录,6.6常微分方程的符号求解MATLAB的符号运算工具箱中提供,例6.30 求微分方程的通解。命令如下：y=dsolve(Dy-(x2+y2)/x2/2,x) %解(1)。方程的右端为0时可以不写y=dsolve(Dy*x2+2*x*y-exp(x),x) %解(2)y=dsolve(Dy-x/y/sqrt
34、(1-x2),x) %解(3),目录,例6.30 求微分方程的通解。目录,例6.31 求微分方程的特解。命令如下：y=dsolve(Dy=2*x*y2,y(0)=1,x) %解(1)y=dsolve(Dy-x2/(1+y2),y(2)=1,x) %解(2),目录,例6.31 求微分方程的特解。目录,例6.32用微分方程的数值解法和符号解法解方程，并对结果进行比较。在MATLAB命令窗口，输入命令：y=dsolve(Dy+2*y/x-4*x,y(1)=2,x) %用符号方法得到方程的解析解为了求方程的数值解，需要按要求建立一个函数文件fxyy.m：function f=fxyy(x,y)f=(4
35、*x2-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.32用微分方程的数值解法和符号解法解方程，并对结果进,MATLAB课件第6章MATLAB符号计算,6.6.3 常微分方程组求解例6.33 求微分方程组的解。命令如下：x,y=dsol
36、ve(Dx=4*x-2*y,Dy=2*x-y,t) %解方程组(1)x,y=dsolve(D2x-y,D2y+x,t) %解方程组(2),目录,x =C1+C2*exp(3*t)y =1/2*C2*exp(3*t)+2*C1x =-C1*exp(1/2*2(1/2)*t)*sin(1/2*2(1/2)*t)-C2*exp(-1/2*2(1/2)*t)*sin(1/2*2(1/2)*t)+C3*exp(1/2*2(1/2)*t)*cos(1/2*2(1/2)*t)+C4*exp(-1/2*2(1/2)*t)*cos(1/2*2(1/2)*t)y =-C1*exp(1/2*2(1/2)*t)*cos(1/2*2(1/2)*t)+C2*exp(-1/2*2(1/2)*t)*cos(1/2*2(1/2)*t)-C3*exp(1/2*2(1/2)*t)*sin(1/2*2(1/2)*t)+C4*exp(-1/2*2(1/2)*t)*sin(1/2*2(1/2)*t),6.6.3 常微分方程组求解目录x =C1+C2*exp(3,
