1 / 44

第 6 章 MATLAB 符号计算

第 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 进行各种运算。. 例 6.1 考察符号变量和数值变量的差别。 在 MATLAB 命令窗口，输入命令：

irish
Télécharger la présentation

第 6 章 MATLAB 符号计算

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. 第6章 MATLAB符号计算 6.1 符号计算基础 6.2 符号导数及其应用 6.3 符号积分 6.4 级数 6.5 代数方程的符号求解 6.6 常微分方程的符号求解

  2. 6.1 符号计算基础 6.1.1 符号对象 1. 建立符号变量和符号常数 (1)sym函数 sym函数用来建立单个符号量，例如，a=sym('a')建立符号变量a，此后，用户可以在表达式中使用变量a进行各种运算。

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

  4. 例6.2比较符号常数与数值在代数运算时的差别。例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)) % 计算数值表达式值

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

  6. 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 %求符号表达式的值

  7. 例6.4计算3阶范得蒙矩阵行列式的值。设A是一个由符号变量a,b,c确定的范得蒙矩阵。例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的行列式值

  8. 例6.5建立x,y的一般二元函数。 在MATLAB命令窗口，输入命令： syms x y; f=sym('f(x,y)');

  9. 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)

  10. (2)因式分解与展开 factor(S) 对S分解因式，S是符号表达式或符号矩阵。 expand(S) 对S进行展开，S是符号表达式或符号矩阵。 collect(S) 对S合并同类项，S是符号表达式或符号矩阵。 collect(S,v) 对S按变量v合并同类项，S是符号表达式或符号矩阵。

  11. 例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的每个元素分解因式

  12. 例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分解因式

  13. (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进行化简，并显示每步结果

  14. 2. 符号矩阵运算 transpose(S) 返回S矩阵的转置矩阵。 determ(S) 返回S矩阵的行列式值。 colspace(S) 返回S矩阵列空间的基。 [Q,D]=eigensys(S) Q返回S矩阵的特征向量，D返回S矩阵的特征值。

  15. 6.1.3 符号表达式中变量的确定 MATLAB中的符号可以表示符号变量和符号常数。findsym可以帮助用户查找一个符号表达式中的的符号变量。该函数的调用格式为： findsym(S,n) 函数返回符号表达式S中的n个符号变量，若没有指定n，则返回S中的全部符号变量。 在求函数的极限、导数和积分时，如果用户没有明确指定自变量，MATLAB将按缺省原则确定主变量并对其进行相应微积分运算。可用findsym(S,1)查找系统的缺省变量，事实上，MATLAB按离字符'x'最近原则确定缺省变量。

  16. 6.2 符号导数及其应用 6.2.1函数的极限 limit函数的调用格式为： limit(f,x,a) limit函数的另一种功能是求单边极限，其调用格式为： limit(f,x,a,'right') 或 limit(f,x,a,'left')

  17. 例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)

  18. 6.2.2 符号函数求导及其应用 MATLAB中的求导的函数为： diff(f,x,n) diff函数求函数f对变量x的n阶导数。参数x的用法同求极限函数limit，可以缺省，缺省值与limit相同，n的缺省值是1。

  19. 例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的偏导数

  20. 例6.12在曲线y=x3+3x-2上哪一点的切线与直线y=4x-1平行。例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

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

  22. 例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函数对结果化简

  23. 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) %将上述符号结果转换为数值

  24. 例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

  25. 例6.16轴的长度为10米，若该轴的线性密度计算公式是f(x)=6+0.3x千克/米(其中x为距轴的端点距离)，求轴的质量。例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)

  26. 例6.17求空间曲线c从点(0，0，0)到点(3，3，2)的长度。求曲线c的长度是曲线一型例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的长度

  27. 6.3.3 积分变换 1. 傅立叶(Fourier)变换 在MATLAB中，进行傅立叶变换的函数是： fourier(fx,x,t) 求函数f(x)的傅立叶像函数F(t)。 ifourier(Fw,t,x) 求傅立叶像函数F(t)的原函数f(x)。

  28. 例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)。

  29. 例6.19计算y=x2的拉普拉斯变换及其逆变换. 命令如下： x=sym('x');y=x^2; Ft=laplace(y,x,t) %对函数y进行拉普拉斯变换 fx=ilaplace(Ft,t,x) %对函数Ft进行拉普拉斯逆变换

  30. 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变换

  31. 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。计算有限级数的和

  32. 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)。

  33. 例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应用泰勒公式近似计算 。 命令如下： x=sym('x'); 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的乘方运算直接计算

  34. 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，不可缺省。

  35. 例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项

  36. 6.5代数方程的符号求解 6.5.1线性方程组的符号求解 MATLAB中提供了一个求解线性代数方程组的函数linsolve，其调用格式为： linsolve(A,b)

  37. 例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)的解

  38. 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)。仅标出方程的左端

  39. 例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)

  40. 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给出求解变量。

  41. 例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)

  42. 例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)

  43. 例6.32用微分方程的数值解法和符号解法解方程，并对结果进行比较。例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-');

  44. 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)

More Related