《第5讲Matlab符号计算.ppt》由会员分享，可在线阅读，更多相关《第5讲Matlab符号计算.ppt（63页珍藏版）》请在三一文库上搜索。
1、第5讲 MATLAB的符号计算,符号计算基础 5.1 符号对象一种新的数据类型 符号对象、符号表达式 5.2 基本的符号运算 表达式运算(四则、化简 、因式分解)、矩阵运算 5.3 符号表达式中变量的确定 符号导数及其应用 5.4 函数的极限 5.5 符号函数求导及其应用 （注：其help见Symbolic Math Toolbox ）,5.6 不定积分 5.7 符号函数的定积分 5.8 积分变换(*) 5.9 级数的符号求和； 5.10 函数的泰勒级数； 5.11 函数的傅立叶级数(*)； 5.12 线性方程组的符号求解； 5.13 非线性方程组的符号求解； 5.14 常微分方程的符号求解；
2、 5.15 常微分方程组求解,符号计算基础,Matlab中数学表达式的三种存在形式及其运用： 、由数值对象参与运算的数学表达式（数值型）： 组成： (数值)常数、数值型变量、函数、数学运算符； 要求：执行表达式之前必须对其中所有数值型变量赋值； 运算结果：一个数值或数值数组 例如：x=1/3 y=exp(x)*sin(2*x)-1/5,、由符号对象参与运算的数学表达式(符号型)： 组成： (符号)常数、符号型变量、符号型常量、函数、数学运算符； 要求：执行表达式之前必须对其中所有符号型变量、符号型常量进行类型定义； 运算结果：一个符号型数学表达式 例如： syms x; y=exp(x)*si
3、n(2*x)-1/5,、以字符串形式表达的数学表达式(字符串型)： 组成： 由数值型或符号型数学表达式两端加 就转换为字符串型数学表达式； 要求： (1)若调用eval函数，可将其转化为数值型数学表达式并执行，但调用eval之前必须对其中所有数值型变量赋值；(2)若调用sym函数，可将其转化为符号型数学表达式；,运算结果： (1)若调用eval函数，可获得一个数值或数值数组; (2)若调用sym函数，获得一个符号型数学表达式; 例如：d=1/3 c=eval(d2-1/5*sin(x) e=sym(d2-1/5*sin(x) class(e) diff(e),5.1 符号对象 1. 建立符号对
4、象（类型定义） (1)syms函数： syms函数的一般调用格式为： syms var1 var2 varn 定义符号对象var1,var2,varn等。用这种格式定义符号量时符号量间用空格而不要用逗号分隔。 还可以用如下方式定义： syms(var1, var2, ,varn) 例如：syms x y z syms(x,y,z),(2)sym函数： sym函数的一般调用格式为： f=sym(arg) 将数值、数值表达式或字符串表达式arg转换为相应的符号对象后赋给f，则f成为一个符号量，它的内容可以是符号常数、符号变量甚至符号表达式。 例如： f=sym(x2-5) diff(f),例1 考
5、察符号变量和数值变量的差别。 在 MATLAB命令窗口，输入命令： syms a b c 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的行列式 class(A) %判断变量A的类型 class(B),例2 比较符号常数与数值在代数运算时的差别。 在 MATLAB命令窗口，输入命令： pi1=sym(pi);k1=sym(8);k2=sym(2);k3=sym(3); % 定义4个符号变量，其内容由4个符号常数
6、指定 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)先定义符号量，再直接写出符号表达式： syms x y; %建立符号变量x、y V=3*x2+5*y+2*x*y+6 %定义符号表达式V (2)用sym函数
7、把字符串表达式转化为符号表达式： U=sym(x3+5*y) %定义符号表达式U U-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确定的范得蒙矩阵。 命令如下： syms a b c; U=a,b,c; A=1,1,1;U;U.2 %建立
8、范得蒙符号矩阵 det(A) %计算A的行列式值,5.2 基本的符号运算 1. 符号表达式运算 (1)符号表达式的四则运算 例6 符号表达式的四则运算示例： syms x y z; 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, 试构造符号函
9、数 命令如下： t=1:5; syms x y; y=0; for k=1:5 s=1; for l=1:k s=s*(x-t(l); end y=y+t(k)*s; end simple(y) %MATLAB自动调用多种函数对s进行化简，并显示每步结果,(3)因式分解与展开 令S是符号表达式或符号矩阵。 factor(S) 对S分解因式。 expand(S) 对S进行展开。 collect(S) 对S合并同类项。 collect(S,v) 对S按变量v合并同类项。,例8 对符号矩阵A的每个元素分解因式。 命令如下： syms a b x y; A=2*a2*b3*x2-4*a*b4*x3+1
10、0*a*b6*x4,3*x*y-5*x2;4,a3-b3; factor(A) %对A的每个元素分解因式,例9 计算表达式S的值。 命令如下： syms x y; 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
11、返回S矩阵的特征值。 例10 a=sym(magic(4);b=sym(ones(4) c=a+b/2 c rank(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也可缺省，缺
12、省时a=0. limit函数的另一种功能是求单边极限，其调用格式为： limit(f,x,a,right) 或 limit(f,x,a,left),例13 求极限,输入命令： 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) %求极限(3) f=(sqrt(x)-sqrt(a)+sqrt(x-a)/sqrt(x*x-a*a); limit(f,x,a,right) %求极限(4),5.5 符号函数求导及其应用 MATLAB中
13、的求导的函数为： diff(f,x,n) 求函数f对变量x的n阶导数。参数x的用法同求极限函数limit，可以缺省，缺省值与limit相同，n的缺省值是1。,例14 求函数的导数：,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的三阶导数 x=a*cos(t); y=b*sin(t); diff(y)/diff(x) %求(3),按参数方程求导公式求y对x的一阶导数 (diff(x)*d
14、iff(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=di
15、ff(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) syms x; 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)
16、%验证 f2_4=simple(f2_3),5.7 符号函数的定积分 定积分在实际工作中有广泛的应用。在MATLAB中，定积分的计算使用函数： int(f,x,a,b) 求函数f对变量x的在区间a,b上的定积分。参数x可以缺省，缺省原则与diff函数相同。 例16求定积分,命令如下： syms x 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) format l
17、ong g double(I) %将上述符号结果转换为数值,例17求椭球的体积 命令如下： syms a b c x; f=pi*c*b*(1-x2/a2); V=int(f,x,-a,a),例18 轴的长度为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(f
18、x,0,10,1e-6),例19 (I型曲线积分)求空间曲线c: 从点(0，0，0)到点(3，3，2)的长度。 命令如下： 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) %计算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(
19、x)。 例 求函数y=|x|的傅立叶变换及其逆变换。 命令如下： 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)。 例 计算y=x2的拉普拉斯变换及其逆变换. 命令如下： x=sym(x);y=x2; Ft=laplace(y,x,t) %对函数y进行拉
20、普拉斯变换 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变换及其逆变换。 命令如下： syms n z fn=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为级数通项
21、表达式，n为求和变量，n0、nn为求和的起点和终点。 例20求级数之和,命令如下： syms n x; s1=symsum(1/n2,n,1,inf) %求(1),常数项级数 s2=symsum(-1)(n+1)/n,1,inf) %求(2),交错级数。未指定求和变量，缺省为n s3=symsum(n*xn,n,1,inf) %求(3),函数项级数。此处的求和变量n不能省略 s4=symsum(n2,1,100) %求(4),计算有限级数的和,5.10 函数的泰勒级数 MATLAB中提供了将函数展开为幂级数的函数taylor，其调用格式为： taylor(f,v,n,a) f为待展开函数表达式
22、，v为自变量，n为展开阶数(正整数)，a则指定对f在v=a处进行泰勒展开。,命令如下： syms 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=5 taylor(f2,6) %求(2)。展开到x的5次幂,例21 求函数的泰勒展开式,例22 将多项式1+3x+5x2-2x3表示成x+1的幂的多项式。 命令如下： syms x; p=1+3*x+5*x2-2*x3; f=taylor(p,x,4,-1),例23 应用泰勒公式近似计算,命令如下： syms
23、x; f=(1-x)(1/12); %定义函数，4000(1/12)=2f(96/212) g=taylor(f,4) %求f的4阶泰勒展开式g v=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 函数的傅立叶级数(*) MATLAB 5.x版中，尚未提供求函数傅立叶级数的内部函数。下面我们自己设
24、计一个简化的求任意函数的傅立叶级数的函数文件。 function mfourier=mfourier(f,n) syms x a b c; mfourier=int(f,-pi,pi)/2/pi; %计算a0 for i=1:n a(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); end return 调用该函数时，需给出被展开的符号函数f和展开项数n，不可缺省。,例24 在-，区间展开函数为傅立叶级数。 命令如下： x=sy
25、m(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项,代数方程的符号求解,5.12 线性方程组的符号求解 对于形如Ax=b的线性方程组，符号求解如下： A=sym(数值型系数矩阵) b=sym(值向量) x=Ab,例25 求线性方程组AX
26、=b的解。,解方程组(1)的命令如下： A=sym(34,8,4; 3,34,3; 3,6,8); b=sym(4;6;2); X=Ab 解方程组(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=Ab,5.13 一般代数方程组的符号求解 求解非线性方程组的函数是solve，调用格式： solve(eqn1,eqn2,eqnN,var1,var2,varN),例26 解非线性方程,命令如下： x=solve(1/(x+2
27、)+4*x/(x2-4)=1+2/(x-2),x) %解方程(1),方程与自变量都为字符串型表达式，结果为符号型 f=sym(x-(x3-4*x-7)(1/3)=1); x=solve(f) %解方程(2)，方程为符号型，自变量缺省 syms x; x=solve(2*sin(3*x-pi/4)-1) %解方程(3)，这一用法不允许“”号 x=solve(x+x*exp(x)-10,x) %解方程(4)，无符号解时，给出数值解,例27 求非线性方程组的解:,命令如下： x y=solve(1/x3+1/y3=28,1/x+1/y=4,x,y) %解方程组(1) s=solve(x2+y2-5,
28、2*x2-3*x*y-2*y2) %解方程组(2)，s是结构数组,5.14 常微分方程的符号求解,MATLAB的符号运算工具箱中提供了功能强大的求解常微分方程的函数dsolve。该函数的调用格式为： dsolve(eqn1,condition,var) 该函数求解微分方程eqn1在初值条件condition下的特解。参数var描述方程中的自变量符号，省略时按缺省原则处理，若没有给出初值条件condition，则求方程的通解。,命令如下： y=dsolve(Dy-(x2+y2)/x2/2,x) %解(1), 方程的右端为0时可以不写 y=dsolve(Dy*x2+2*x*y-exp(x),x)
29、%解(2) y=dsolve(Dy-x/y/sqrt(1-x2),x) %解(3),例28 求微分方程 的通解,命令如下： y=dsolve(Dy=2*x*y2,y(0)=1,x) %解(1) y=dsolve(Dy-x2/(1+y2),y(2)=1,x) %解(2),例29 求微分方程的特解,例30 用数值解法和符号解法解微分方程并比较结果 (1)用数值运算求方程的数值解： 首先，要建立一个函数文件fxyy.m： function f=fxyy(x,y) f=(4*x2-2*y)/x; %形式只能是y=f(x,y) ,不是这种形式要变形 输入数值函数命令： t,w=ode45(fxyy,1,
30、2,2) %得到区间1，2中的数值解，以向量t、w存储。,(2) 用符号运算求方程的解析解(即符号解)： y=dsolve(Dy+2*y/x-4*x, y(1)=2, x) (3) 作出两种结果的图形,对它们进行比较： x=linspace(1,2,100); y=eval(vectorize(char(y) %为作图把符号解的结果离散化;vectorize函数见P344,213 plot(x,y,b.,t,w,r-);,dsolve在求微分方程组时的调用格式为： dsolve(eqn1,eqn2,eqnN,condition1,conditionN,var1,varN) 函数求解微分方程组eqn1、eqnN在初值条件conditoion1, conditionN下的解，若不给出初值条件，则求方程组的通解，var1、varN给出求解变量。,5.15 常微分方程组求解,命令如下： x,y=dsolve(Dx=4*x-2*y,Dy=2*x-y,t) %解方程组(1) x,y=dsolve(D2x-y,D2y+x,t) %解方程组(2),例31 求微分方程组的解,