山大MATLAB编程指导课件第2章 符号计算

1、第2章 符号计算教学目标教学重点教学内容教学目标一是讲述MATLAB符号计算基本知识，包括符号对象的创建、符号数字、符号表达式的操作； 二是介绍符号微积分的计算；三是介绍符号矩阵分析和代数方程（组）的符号解法；四是介绍符号计算结果的可视化。教学重点熟悉符号对象的创建、符号数字、符号表达式的操作。熟悉符号微积分的基本计算函数指令。熟悉代数方程（组）的符号解法。熟悉符号计算结果可视化的的基本指令。了解符号计算帮助系统及其帮助指令。 教学内容2.1符号对象和符号表达式2.2 符号数字及表达式的操作 2.3 符号微积分2.4 微分方程的符号解法2.5 符号变换和符号卷积 2.6 符号矩阵分析和代数方程

2、解 2.7 代数状态方程求符号传递函数 2.8 符号计算结果的可视化 2.9 符号计算资源深入利用Matlab的符号计算功能matlab自产生起就在数值计算上功能卓著,深受各专业计算人员的欢迎.但由于在数学,物理等各种科研和工程应用中经常会遇到符号运算的问题. 为此, 公司于1993年购买了 Maple 软件的使用权,并在此基础上,开发了符号计算工具箱 (Symbolic Math Toolbox)matlab 从2008b 开始与符号计算语言MuPAD 相结合，到2009b止仍然支持Maple引擎（需单独安装Maple软件） 。在此版本之间输入指令symengine弹出选择MuPAD和 Ma

3、ple引擎的窗口。从2010a开始不再支持Maple引擎。符号运算与数值运算的区别：符号运算中，解算数学表达式、方程时，不是在离散化的数值点上进行，而是凭借一系列恒等式和数学定理，通过推理和演绎，获得解析结果。这种计算建立在数值完全准确表达和推演严格解析的基础上，所得结果是完全准确的。符号运算代数运算，公式推导数值运算算术运算代值2.1 符号对象和符号表达式在matlab中,数值和数值变量用于数值的存储和各种数值计算.而符号常量,符号变量,符号函数,符号操作等则是用来形成符号表达式,严格按照代数,微积分等课程中的规则,公式进行运算,并尽可能给出解析表达式.2.1.1 符号对象的创建和衍生 数值

4、计算变量先赋值,再使用. 符号计算先定义基本的符号对象(可以是常量,变 量,表达式),然后用这些基本符号对象去构成新的表达式,再进行所需的符号运算2.1.1 符号对象的创建和衍生1. 生成符号对象的基本规则 任何基本符号对象（数字、参数、变量、表达式）都必须借助专门的符号函数指令sym或syms定义。 任何包含符号对象的表达式或方程，将继承符号对象的属性。即任何包含符号对象的表达式、方程也一定是符号对象。2 符号数字的定义格式：sc=sym(num) % sc为值为num的符号数字注意： i) 单引号必须在英文状态下输入，构成字符串 ii) num为一个具体的数字如： sc=sym(2/3)

5、sb=sym(pi+sqrt(5)sc=2/3sb =pi + 5(1/2)2 符号数字的定义【例2.1-1】符号（类）数字与数值（类）数字之间的差异。a=pi+sqrt(5)sa=sym(pi+sqrt(5)Ca=class(a)Csa=class(sa)vpa(sa-a) a = 5.3777sa = pi + 5(1/2)Ca = doubleCsa = symans =0.00000000000000001382237584108520004859354256418 本例表现符号数字总是被准确记录和运算，而数值数字并不总能保证被准确存储，运算时会引进截断误差。3. 基本符号变量：定义格

6、式：i) syms para para=sym(para) syms a; a=sym(a)ii) syms para flag para=sym(para, flag) syms a positive; a=sym(a, positive)flag为参数属性：positive参数取正实数real参数为实数unreal参数为限定的复数iii) syms a b c syms a b c flag无逗号符号参数（表达式中的参数）, “待解符号变量”或“自由符号变量” （表达式中的自变量x，默认为x）4. 自由符号变量symvar(expression) 列出表达式中所有基本符号变量symvar(

7、expression,n) 列出表达式中认定n个自由符号变量expression是符号表达式或符号表达式矩阵，x是首选自由符号变量，认定优先次序为x, y, w, z, v等解题结果是“用符号参数构成的表达式表述自由符号变量”。解题时自由符号变量可“人为指定”，也可“默认地自动认定”：与小写字母 x 的ASII码距离最小的变量。syms u v w z a5f=sym(3);Eq=sin(f)*u*z2+v*z+f*w-a5;【例2.1-2】 用符号计算研究方程 的解。symvar(Eq) %按字母表顺序列出基本符号变量, 无 fans = a5, u, v, w, z symvar(Eq,1

8、00) %按离x的距离列出所有自由符号变量ans = w, z, v, u, a5result_1=solve(Eq) result_1 =a5/3 - (v*z)/3 - (u*sin(3)*z2)/3 result_2=solve(Eq,z)result_2 = -(v - (v2 + 4*a5*u*sin(3) - 12*u*w*sin(3)(1/2)/(2*u*sin(3) -(v + (v2 + 4*a5*u*sin(3) - 12*u*w*sin(3)(1/2)/(2*u*sin(3)syms a b x X Yk=sym(3);z=sym(c*sqrt(d)+y*sin(t);E

9、XPR=a*z*X+(b*x2+k)*Y;【例2.1-3】元符号表达式、衍生符号表达式定义，基本符号变量、自由符号变量的机器认定。E3=sym(a*sqrt(theta) ? Error using = sym.sym convertExpression at 2515 E4=sym(a*sqrt(theta1)E5=sym(a*sqrt(theta*t) %在R2009b版本中还正确? Error using = sym.sym convertExpression at 2515symvar(EXPR)ans = X, Y, a, b, c, d, t, x, y 无 k zsymvar(E

10、XPR,10) ans = x, y, t, d, c, b, a, X, Y 4. 自由符号变量syms a b t u v x yA=a+b*x, sin(t)+u; x*exp(-t), log(y)+vsymvar(A,1) A = a + b*x, u + sin(t) x/exp(t), v + log(y)ans = x 【例2.1-4】symvar确定自由变量是对整个矩阵进行的。 2.1.2 符号计算中的算符由于新版matlab采用了重载(Overload)技术, 使得用来构成符号计算表达式的算符和基本函数,无论在形式,名称,还是使用方法上,都与数值计算中的算符和基本函数几乎完

11、全相同,这给编程带来极大的方便.(1) 基本运算符算符 ”+”, ”-”, ”*”, ”, “/”,“” 分别构成矩阵的加, 减, 乘,左除,右除,求幂运算. 算符 ”.*”, “./”, “.”, “.” 分别实现元素对元素的数组乘,除,求幂运算.算符” ”, “ . ” 分别实现矩阵的共轭转置,非共轭转置2.1.2 符号计算中的算符(2) 关系运算符在符号对象的比较中,没有大于,大于等于,小于,小于等于的概念,而只有是否等于的概念。 ”=“ “=“分别用来对算符两边的对象进行相等和不等的比较,返回为逻辑量。事实为真时，比较结果为1，事实为假时，结果为0.2.1.3 符号计算中的函数指令符号

12、计算中的函数分成三个层次：与数值类函数和指令对应的同名符号类函数和指令。约50个经典特殊函数（误差函数erf、贝塞尔函数besselj、椭圆积分ElliptiK等），借助mfun调用，用mfunlist可列出。数量较大的MuPAD库函数，借助evalin和feval调用。(2) 指数,对数函数 sqrt, exp, expm在两者中用法相同符号计算中只有自然对数log，而没有数值计算中的log2, log10(3) 复数函数 conj, imag, real, abs在两者中用法相同.但在符号计算中没有求相角的指令angle.(1) 三角函数,双曲线函数以及他们的反函数除atan2只能用于数值

13、计算外,另外的在两种运算中使用方法相同.与数值类函数和指令对应的同名符号类函数和指令(4) 矩阵代数指令在符号计算中, matlab提供的常用矩阵代数指令有：diag, tril, inv, det, rank, eig, svd（ Singular value decomposition奇异值分解）等与数值类函数和指令对应的同名符号类函数和指令(5) 方程求解指令solve，与数值类不同。(6) 微积分如diff，int，与数值类不完全相同。(7) 积分变换和反变换函数如laplace， ilaplace，数值类只有Fourier变换。(8) 绘图函数如ezplot，ezsurf，数值类绘图

14、指令更丰富。数值计算对象,符号计算对象,字符串是MATALB中最常用的数据对象.他们遵循各自不同的运算法则,但有时在外形上却十分相似.MATLAB提供了一些识别不同数据对象的指令,常用的有class, isa, whos 例2.1.3-1 数据对象及其识别指令的使用(1) 生成三种不同类型的矩阵，给出不同的显示形式clear, a=1;b=2;c=3;d=4 % 产生四个数值变量Mn=a,b; c,d % 利用已赋值变量构成数值矩阵Mc=a,b; c,d % 字符串中的a,b,c,d与前面输入的数值变量无关Ms=sym(Mc) % 符号变量，与前面的各变量无关2.1.4 符号对象的识别Mn =

15、 1 2 3 4Mc =a,b;c,dMs = a, b c, dSizeMn=size(Mn), SizeMc=size(Mc), SizeMs=size(Ms)SizeMn = 2 2SizeMc = 1 9SizeMs = 2 2(3) 用class获得每种矩阵的类别CMn=class(Mn), CMc=class(Mc), CMs=class(Ms)CMn = doubleCMc = charCMs = sym(4) 用isa判断每种矩阵的类别isa(Mn,double), isa(Mc,char), isa(Ms,sym)ans = 1ans = 1ans = 1Mn = 1 2 3

16、 4 Mc = a,b;c,dMs = a, b c, d(2) 三个矩阵的大小不同(5) 利用whos观察内存变量的类别和其他属性 Name Size Bytes Class Mc 1x9 18 char array Mn 2x2 32 double array Ms 2x2 312 sym objecta=0:1:6; theta=sym(30*pi/180*a)alfa=sym(30*pi/180)*atheta =30*pi/180*aalfa = 0, 1/6*pi, 1/3*pi, 1/2*pi, 2/3*pi, 5/6*pi, pia与数组a无关whos P26 Name Siz

17、e Bytes Class ans 1x1 8 double t 1x201 1608 double y 1x201 1608 double 2.1.5 符号运算机理和变量假设1.符号运算的工作机理Matlab借助sym或syms定义符号变量时，启动MuPAD引擎并启动一个专供MuPAD使用的内存工作空间执行符号运算；matlab内存空间只保存该符号变量和计算结果。该定义变量的限定性假设（assumption）被保存在MuPAD的内存空间。若变量不带限定性假设，则MuPAD默认为复数。2.1.5 符号运算机理和变量假设2. 对符号变量的限定性假设i) syms x para=sym(x)ii)

18、 syms x flag para=sym(x, flag) syms a positive; a=sym(a, positive)iii) syms a b c syms a b c flagflag为参数属性：positive参数取正实数real参数为实数unreal参数为限定的复数2.1.5 符号运算机理和变量假设3. 清除变量和撤销假设符号变量和其假设存放在不同的内存空间，因此删除符号变量和撤销关于变量的假设是两件需要分别处理的事情。 clear all 可同时删除clear x 删除matlab内存中的x变量syms x clear 撤销MuPAD内存中对变量x的假设evalin(s

19、ymengine,getprop(x) ) 获取x的限定性假设 evalin(symengine,anames (Properties) 列出MuPAD内存中带限定性假设的符号变量reset(symengine) 重启MuPAD引擎，清空MuPAD内存clear all 删除matlab及 MuPAD内存中的所有变量【例2.1-6】syms 对符号变量限定性假设的影响syms x clearf=x3+4.75*x+2.5;rf=solve(f,x)rf = -1/2 1/4 - (79(1/2)*i)/4 (79(1/2)*i)/4 + 1/4 evalin(symengine,getprop

20、(x) ans =C_ syms x realrfr=solve(f,x) rfr =-1/2 evalin(symengine,getprop(x) ans =R_ clear xsyms xg=x2+x+5;rg=solve(g,x)Warning: Explicit solution could not be found. In solve at 98rg = empty sym evalin(symengine,anames(Properties) ans =xsyms x clearrg=solve(g,x) rg = - (19(1/2)*i)/2 - 1/2 (19(1/2)*i

21、)/2 - 1/2 reset(symengine)clear all2.1.6 符号帮助体系 Matlab符号指令帮助系统与2008年前不同，因为用MuPAD取代maple作为符号计算的内核。help SymNamehelpwin SymNamedoc SymNamedoc mfunlistdoc(symengine)2.1.6 符号帮助体系 Matlab符号指令帮助系统与2008年前不同，因为用MuPAD取代maple作为符号计算的内核。help SymNamehelpwin SymNamedocsearch laplacedoc mfunlistdoc(symengine)2.1.6 符

22、号帮助体系 Matlab符号指令帮助系统与2008年前不同，因为用MuPAD取代maple作为符号计算的内核。help SymNamehelpwin SymNamedocsearch laplacedoc mfunlistdoc(symengine)开启MuPAD帮助浏览器。2.2 符号数字及表达式的操作2.2.1双精度数字与符号数字之间的转换sym(num, r) % “有理分数”表达的符号数字 : p/q, n(p/q)sym(num) % sym(num, r) 的简写形式1. 双精度数字向符号数字的转换 sym(num, f) %以数值 表示“浮点数”，N, e为整数 sym(1/10

23、,f) is 3602879701896397/36028797018963968 sym(num, e) % “有理分数”表达+机器实际浮点表达误差e sym(3*pi/4,e) is 3*pi/4 - 103*eps/249sym(num, d) %十进制小数近似表示，有效数字位数受 digits指令控制。默认为digits(32)情况.2.2.1双精度数字与符号数字之间的转换sym(num) % sym(num, r) 的简写形式sym(num) % num是字符串数字 1. 双精度数字向符号数字的转换 sym(num) 中num 用作符号计算函数的输入时，体现其理论真值，对sym(nu

24、m)中的num用作符号计算函数的输入时，体现其字面数子的双精度近似。2.2.1双精度数字与符号数字之间的转换double(num_sym) % 把符号数字转换为双精度数字2. 符号数字向双精度数字转换 注意：double(num) 把字符串数字转换为字符的ASCII码double(0.1)ans = 48 46 49f=sym(10/3)Df=double(f)class(Df) Df = 3.3333ans = doubledigits 显示当前环境下符号数字“十进制浮点” 表示的有效数字位数；digits(n) 设定 “十进制浮点”表示的有效数字位数； xs=vpa(x) 据表达式x得到d

25、igits指定精度下的符号数字xs xs=vpa(x,n) 据表达式x得到n位有效数字的符号数字xs2.2.2 符号数字的任意精度表达形式 数值计算与符号计算的最重要区别：数值计算存在截断误差，且在计算中不断传播，产生累计误差；符号计算完全准确，没有累计误差，但降低计算速度，增加内存。为兼顾精度和速度，采用“变精度”算法。reset(symengine) sa=sym(1/3+sqrt(2) sa =2(1/2) + 1/3 digits Digits = 32 format longa=1/3+sqrt(2)sa_Plus_a=vpa(sa+a,20)sa_Minus_a=vpa(sa-a,

26、20)a = 1.747546895706428sa_Plus_a =3.4950937914128567869sa_Minus_a =-0.000000000000000022658064826339973669 例2.2-1 digits, vpa, symengine指令演示sa32=vpa(sa)digits(48)sa5=vpa(sa,5)sa48=vpa(sa)sa32 =1.747546895706428382135022057543sa5 =1.7475sa48 =1.74754689570642838213502205754303141190300520871 2.2.3 符

27、号表达式的基本操作collect合并同类项numden获取最小公分母和相应分子expand展开指定项simplify简化表达式，消除冗余项factor因式分解simple搜索寻找最简表达horner转换为嵌套形式pretty习惯方式显示syms x;f=(1/x3+6/x2+12/x+8)(1/3);g1=simple(f)g1 =(2*x + 1)3/x3)(1/3)验证:f2=g13f2 =(2*x + 1)3/x3expand(f2)ans =12/x + 6/x2 + 1/x3 + 8例2.2-2 简化2.2.4 表达式中的置换操作1.公共子表达式法简化表达 RS=subexpr(S)

28、; RS=subexpr(S,w); RS,w=subexpr(S,w) 其中，w为MATLAB自动寻找的子表达式syms a b c d W;V,D=eig(a,b;c,d);RVD,W=subexpr(V;D,W)【例2.2-3】 对符号矩阵 进行特征向量分解。RVD = (a/2 + d/2 - w/2)/c - d/c, (a/2 + d/2 + w/2)/c - d/c 1, 1 a/2 + d/2 - w/2, 0 0, a/2 + d/2 + w/2w =(a2 - 2*a*d + d2 + 4*b*c)(1/2) V =(a/2+d/2-(a2-2*a*d+d2+4*b*c)(

29、1/2)/2)/c-d/c, (a/2+d/2+(a2-2*a*d+ d2+4*b*c)(1/2)/2)/c-d/c 1, 1D =a/2 + d/2 - (a2 - 2*a*d + d2 + 4*b*c)(1/2)/2, 0 0, a/2 + d/2 + (a2 - 2*a*d + d2 + 4*b*c)(1/2)/2 2.2.4 表达式中的置换操作2. 通用置换指令RESsubs(ES,old,new)RESsubs(ES,new)例2.2-4 subs置换规则示例clearsyms a b x;f=a*sin(x)+b f1=subs(f,sin(x),log(y)f2=subs(f,a

30、,3.11)f3=subs(f,a,b,x,2,5,sym(pi/3)class(f3) f1=a* log(y) +b f2 =b + (311*sin(x)/100f3 =3(1/2) + 5ans =sym (5)例2.2-4 subs置换规则示例format %设为默认输出格式显示format compactf4=subs(f,a,b,x,2,5,pi/3)class(f4) f4 = 6.7321ans =double f5=subs(f,x,0:pi/2:pi)class(f5)f5 = b, a + b, bans =sym t=0:pi/10:2*pi;f6=subs(f,a,

31、b,x,2,3,t)plot(t,f6)k=(0.5:0.1:1);f6=subs(subs(f,a,b,k,2),x,t);size(f6) plot(t,f6) ans = 6 21 f6=2*sin(t)+3 f=a*sin(x)+b 符号限定假设对方程根不起作用。Matlab6.5 版本上机问题findsym(expression,n) 当n大于实际的基本变量数目时，按字母表顺序列出所有本符号变量；当n小于等于时，按与x距离顺序列出。symvar(expression,n) 多一个参数n在Matlab6.5中不能用syms x clearf=x3+4.75*x+2.5;rf=solve

32、(f,x)rf = -1/2 1/4 - (79(1/2)*i)/4 (79(1/2)*i)/4 + 1/4 syms x realrfr=solve(f,x) rfr =-1/2 evalin(symengine,getprop(x) ans =R_ evalin(symengine,anames(Properties)调用word问题。指令的输入、编写用M文件编辑器。习题2 (Page114) 除了第3题不能得到实根、正根外，112题, 2325题都可做，只是4题结果在差值有效数字上与2010版本有差别, 7题不能得到最简结果,9题结果 6.5版本错误,。 单周上机，鞋套带上。 clear

33、 all; clear mex syms x unreal通用置换指令RESsubs(ES,old,new)RESsubs(ES,new)digits(n) 设定 “十进制浮点”表示的有效数字位数； xs=vpa(x,n) 据表达式x得到n位有效数字的符号数字xsdouble(num_sym) % 把符号数字转换为双精度数字sym(num) %双精度数字转符号数字sym(num, r) 的简写形式sym(num) %定义符号数字， num是字符串数字 syms x flag para=sym(x, flag) syms a positive; a=sym(a, positive)reviewd

34、fdvn=diff(f,v,n) 求fjac=jacobian(f,v) 求多元向量函数f(v)的jacobian矩阵r=taylor(f,n,v,a) 把f(v)在v=a处进行泰勒展开2.3 符号微积分2.3.1 极限和导数的符号计算limit(f,v,a) 求极限limit(f,v,a,right) 求右极限limit(f,v,a,left) 求左极限【例2.3-1】试求syms x kf=(1-1/x)(k*x); Lf=limit(f,x,inf) Lf =1/exp(k) Lf1=vpa(subs(Lf,k,sym(-1),48) Lf1=2.71828182845904523536

35、02874713526624977572470937 syms a t x f=a, t3; t*cos(x), log(x); df=diff(f) dfdt2=diff(f,t,2) dfdxdt=diff( diff(f,x) ,t) 【例2.3-2】 求df = 0, 0 -t*sin(x), 1/xdfdt2 = 0, 6*t 0, 0dfdxdt = 0, 0 -sin(x), 0 例2.3-5:设cos(x+siny)=siny, 求dy/dx(隐函数求导).%将dy/dx表达式用x,y表达syms xg=sym(cos(x+sin(y(x)=sin(y(x)dgdx=diff(

36、g,x) g =cos(x + sin(y(x) = sin(y(x)dgdx = -sin(x + sin(y(x)*(cos(y(x)*diff(y(x), x) + 1) = cos(y(x)*diff(y(x), x) dgdx1=subs(dgdx,diff(y(x),x),dydx)dgdx1 =-sin(x + sin(y(x)*(dydx*cos(y(x) + 1) = dydx*cos(y(x) dydx=solve(dgdx1, dydx) dydx =-sin(x + sin(y(x)/(cos(y(x) + cos(y(x)*sin(x + sin(y(x)sym xr

37、=taylor(x*exp(x),9,x,0)pretty(r) 例2.3-6:求f(x)=xex在x=0处展开的8阶Maclaurin级数，即忽略9阶及以上小量的泰勒级数展开 。 r =x8/5040 + x7/720 + x6/120 + x5/24 + x4/6 + x3/2 + x2 + x 8 7 6 5 4 3 x x x x x x 2 + + + - + - + - + x + x 5040 720 120 24 6 2 MATLAB求解通式求和 问题的指令为：s=symsum(f,v,a,b) 求通式f在指定变量v取遍 a,b中所有整数时的和。 2.3.2 序列/级数的符号求

38、和说明： （1）f是矩阵时，求和对逐个元素进行，但自变量定义在整个矩阵上。（2）v省缺时，f中的自变量由findsym（symvar）自动辨认；b可以取有限整数也可以取无穷大。（3）a，b同时省缺时，默认的自变量求和区间为0,v-1syms n kf1=1/(k*(k+1);s1=symsum(f1,k,1,n)s1 =1 - 1/(n + 1) 【例2.3-8】求 ， 。f2=1/(2*k-1)2, (-1)k/k s2=symsum(f2,1,inf) s3 = pi2/8, -log(2) 2.3.3 符号积分【例2.3-9】求 。intf=int(f,v) 给出f 对指定变量v的不定积

39、分intf=int(f,v,a,b) 给出f对指定变量v的定积分syms xf= x*log(x) s=int(f,x)s=simple(s)f = x*log(x)s = (x2*(log(x) - 1/2)/2s = x2*(log(x)/2 - 1/4)2.3.3 符号积分intf=int(f,v) 给出f 对指定变量v的不定积分intf=int(f,v,a,b) 给出f对指定变量v的定积分syms a b x;f2=a*x, b*x2; 1/x, sin(x)disp(The integral of f is)pretty(int(f2) 【例2.3-10】求 。The integra

40、l of f is +- -+ | 2 3 | | a x b x | | , | | 2 3 | | | | log(x), -cos(x) | +- -+ f2 = a*x, b*x2 1/x, sin(x) 例2.3-11求积分syms x y z;f=x2+y2+z2;F2=int(int(int(f, z, sqrt(x*y), x2*y), y, sqrt(x), x2), x, 1, 2)第一重积分第二重积分Warning: Explicit integral could not be found. F2 =(14912*2(1/4)/4641 - (6072064*2(1/2)

41、/348075 + (64*2(3/4)/225 + 1610027357/6563700VF2=vpa(F2) %默认32位有效数字 VF2 =224.92153573331143159790710032805教材中48位有效数字 syms a r theta phi r=a*theta; x=r*cos(theta);y=r*sin(theta); dLdth=sqrt(diff(x,theta)2+diff(y,theta)2);L=simple(int(dLdth,theta,0,phi)例2.3-12 求阿基米德螺线r=a*(a0)在=0到 间的曲线长度函数，并求出a=1, 时的曲线

42、长度。解:L_2pi=subs(L,a,phi,sym(1,2*pi) L_2pi_vpa=vpa(L_2pi) L =(phi*(a2*phi2 + a2)(1/2) + log(phi + (phi2 + 1)(1/2)*(a2)(1/2)/2 L_2pi_vpa =21.256294148209098800702512272566例2.3-12 画阿基米德螺线r = a*及其曲线长度图形。L1=subs(L,a,sym(1);ezplot(L1*cos(phi),L1*sin(phi),0,2*pi)grid on; hold onx1=subs(x,a,sym(1);y1=subs(y

43、,a,sym(1);h1=ezplot(x1,y1,0,2*pi);set(h1,Color,r,LineWidth,5) title( ); legend(螺线长度-幅角曲线,阿基米德螺线)hold offL =(phi*(a2*phi2 + a2)(1/2) + log(phi + (phi2 + 1)(1/2)*(a2)(1/2)/2 见P1982.4 微分方程的符号解法S=dslove(eq1,eq2,eqn, cond1,cond2,condn,v)2.4.1 符号解法和数值解法有互补作用2.4.2 求微分方程符号解的一般指令从数值计算角度看，微分方程边值问题的求解比初值问题复杂困难

44、；而对符号计算来说，边值问题和初值问题的求解微分方程的指令形式相同且简单，不过计算时间较长且未必有封闭形式解。S=dslove(eq1,eq2,eqn,cond1,cond2,condn,v)默认独立变量t输入量为字符串必选可选条件：y(a)=b, Dy(c)=d若应变量为y ，用“Dny”表示“y的n阶导数”， Dy为一阶导数。解在S.y中。说明：（1）输入量包括三部分：微分方程、初始条件、指定独立变量（不指定时，默认为t）。输入量必须以字符串的形式给出。（2）微分方程的记述规定：当y是应变量时， 用“Dny”表示“y的n阶导数”。（3）关于初始条件或边界条件的规定：应写成y(a)=b, D

45、y(c)=d等。（4）如果y是应变量，则它的解在S.y中。【例2.4-1】求 的解。2.4.3 微分方程符号解示例clearS=dsolve(Dx=y,Dy=-x) disp(微分方程的解, blanks(2),x,blanks(22),y)disp(S.x,S.y) S = x: 1x1 sym y: 1x1 sym微分方程的解 x y C2*cos(t) + C1*sin(t), C1*cos(t) - C2*sin(t)写微分方程遵循原则：导数在前函数在后，导数阶数降阶x, y= dsolve(Dx=y,Dy=-x) %按字母表输出顺序不变例2.4-2.图示微分方程y=xy-(y)2通解

46、和奇解的关系clear ally=dsolve(Dy)2-x*Dy+y=0,x)y = x2/4 C3*x - C32 clf, hold onhy1=ezplot(y(1),-6,6,-4,8,1);set(hy1,Color,r,LineWidth,5)for k=-2:0.5:2 y2=subs(y(2),C3,k); ezplot(y2,-6,6,-4,8,1)endhold off;box onlegend(奇解,通解,Location,Best)ylabel(y)title(fontsize14微分方程, (y )2 xy + y = 0 ,的解) 见P199例2.4-3.求解两点

47、边值问题: xy-3y=x2,y(1)=0,y(5)=0y=dsolve(x*D2y-3*Dy=x2,y(1)=0,y(5)=0,x) y =(31*x4)/468 - x3/3 + 125/468 ezplot(y,-1,6)hold onplot(1,5,0,0,.r,MarkerSize,20)text(1,1,y(1)=0)text(4,1,y(5)=0)title(x*D2y-3*Dy=x2, y(1)=0,y(5)=0)hold off2.5 符号变换和符号卷积2.5.1 Fourier变换及其反变换Fw=fourier(ft,t,w) 求“时域”函数ft的Fourier变换ft=

48、ifourier(Fw,w,t) 求“频域”函数Fw的Fourier 反变换2.5.2 Laplace变换及其反变换Fs=laplace(ft,t,s) 求“时域”函数ft的Laplace变换ft=ilaplace(Fs,s,t) 求“频域”函数Fs的Laplace 反变换2.5.3 Z变换及其反变换FZ=ztrans(ft,n,z) 求“时域”序列ft的Z变换fn=itrans(FZ,z,n) 求“频域”序列FZ的Z反变换2.5.4 符号卷积 利用拉氏变换求取。【例 2.5-1】求 的Fourier变换。（1）求Fourier变换syms t wut=heaviside(t);UT=four

49、ier(ut) UT = pi*dirac(w)-i/w （2）求Fourier反变换Ut=ifourier(UT,w,t) Ut =heaviside(t)plot(t(kk),ut(kk),.r,MarkerSize,30)ut(kk)=NaN; hold onplot(t,ut,-r,LineWidth,3)plot(t(kk),t(kk),ut(kk-1), ut(kk+1),or,MarkerSize,10)hold offgrid onaxis(-2,2,-0.2,1.2)xlabel(fontsize14t),ylabel(fontsize14ut)title(fontsize1

50、4Heaviside(t) （3）单位阶跃曲线t=-2:0.01:2;ut=heaviside(t);kk=find(t=0);【例2.5-4】求 的Laplace变换。 syms t s; syms a b positive; %不限定无结果 Dt=dirac(t-a); Ut=heaviside(t-b); Mt=Dt,Ut;exp(-a*t)*sin(b*t),t2*exp(-t); MS=laplace(Mt,t,s) MS = exp(-s*a), exp(-s*b)/s 1/b/(s+a)2/b2+1), 2/(s+1)3 dfdvn=diff(f,v,n) 求r=taylor(f

51、,n,v,a) 把f(v)在v=a处进行泰勒展开review2.3.1 极限和导数的符号计算limit(f,v,a) 求极限limit(f,v,a,right) 求右极限limit(f,v,a,left) 求左极限2.3.2 序列/级数的符号求和s=symsum(f,v,a,b)S=dslove(eq1,eq2,eqn, cond1,cond2,condn,v)2.4.2 求微分方程符号解的一般指令S=dslove(eq1,eq2,eqn,cond1,cond2,condn,v)默认独立变量t输入量为字符串必选可选条件：y(a)=b, Dy(c)=d若应变量为y ，用“Dny”表示“y的n阶导

52、数”， Dy为一阶导数。解在S.y中。reviewintf=int(f,v) 给出f 对指定变量v的不定积分intf=int(f,v,a,b) 给出f对指定变量v的定积分2.3.3 符号积分例2.4-2.图示微分方程y=xy-(y)2通解和奇解的关系clear ally=dsolve(Dy)2-x*Dy+y=0,x)clf, hold onhy1=ezplot(y(1),-6,6,-4,8,1);set(hy1,Color,r,LineWidth,5)for k=-2:0.5:2 y2=subs(y(2),C3,k); ezplot(y2,-6,6,-4,8,1)endhold off;box

53、 onlegend(奇解,通解,Location,Best) %参数的设置有多种ylabel(y)title(fontsize14微分方程, (y )2 xy + y = 0 ,的解)legend(奇解,通解,Location,northwest)legend(奇解,通解,Location, 0.4,0.6,0.3,0.1)legend(.,Location,location) location can be either a 1-by-4 position vector (left bottom width height) or one of the following strings. 见

54、P199docsearch legend Location %通过帮助指令获取参数SpecifierLocation in AxesNorthInside plot box near topSouthInside bottomEastInside rightWestInside leftNorthEastInside top right (default for 2-D plots)NorthWestInside top leftSouthEastInside bottom rightSouthWestInside bottom leftNorthOutsideOutside plot box

55、 near topSouthOutsideOutside bottomEastOutsideOutside rightWestOutsideOutside leftNorthEastOutsideOutside top right (default for 3-D plots)NorthWestOutsideOutside top leftSouthEastOutsideOutside bottom rightSouthWestOutsideOutside bottom leftBestLeast conflict with data in plotBestOutsideLeast unuse

56、d space outside plot2.6 符号矩阵分析和代数方程的解2.6.1 符号矩阵分析最常用的矩阵分析指令如下：det(A)diag(A)V,D=eig(A)expm(A) 矩阵指数函数eAinv(A)V,J=jordan(A) Jordan分解，使AV=VJpoly(A)rank(A)2.6.1 符号矩阵分析syms a11 a12 a21 a22A=a11, a12; a21, a22DA=det(A)IA=inv(A)A = a11, a12 a21, a22DA =a11*a22 - a12*a21IA = a22/(a11*a22 - a12*a21), -a12/(a1

57、1*a22 - a12*a21) -a21/(a11*a22 - a12*a21), a11/(a11*a22 - a12*a21) EA=subexpr(eig(A),D) D = (a112 - 2*a11*a22 + a222 + 4*a12*a21)(1/2)EA = a11/2 + a22/2 - D/2 a11/2 + a22/2 + D/2 【例2.6-1】求矩阵 的行列式、逆和特征根。例2.6-2 Givens旋转 对矩阵 (两点)的旋转作用.解: syms t; A=sym(sqrt(3)/2,1/2;1/2,sqrt(3)/2); %符号变量G=cos(t),-sin(t)

58、;sin(t),cos(t); GA=G*AAn=subs(GA,t,110/180*pi); %旋转后数值阵Op=0;0; Ad=double(A); %旋转前数值阵v1=Op,Ad(:,1), v2=Op,Ad(:,2) %旋转前各点与坐标原点相连直线阵u1=Op,An(:,1), u2=Op,An(:,2) %旋转后plot(v1(:,1),v1(:,2),-k,v2(:,1),v2(:,2),b), hold onhu=plot(u1(:,1),u1(:,2),-k,u2(:,1),u2(:,2),b)set(hu,LineWidth,4)Lstr=旋转前的 v1;旋转前的 v2;旋转

59、后的 u1;旋转后的 u2; legend(Lstr,Location,South)axis(-1,1,-1,1), axis square, hold off, grid on2.6.2 线性方程组的符号解（参考4.2节）矩阵计算是求解线性方程组最简单有效的方法,在matlab中,不管数据对象是数值还是符号,实现矩阵运算的指令形式几乎完全相同. A=sym(1,1/2,1/2 , -1; 1, 1, -1,1; 1, -1/4, -1,1; -8, -1, 1, 1);b=sym(0; 10; 0; 1)x=Ab【例2.6-3】求线性方程组的解。X = 1 8 8 9 A=1,1/2,1/2

60、 , -1; 1, 1, -1,1; 1, -1/4, -1,1; -8, -1, 1, 1;b=0; 10; 0; 1;x=AbX=1.0000 8.0000 8.0000 9.0000 2.6.3 一般代数方程组的解一般代数方程包括线性 ( Linear )、非线性( Nonlinear )和超越方程( Transcendental equation ) 等。S =solve(eq1,eq2,.,eqN,v1,v2,.vn) (推荐)S =solve(exp1, exp2,., expN, v1, v2,. vn) (可用)只能是符号表达式符号变量【例2.6-3】求线性方程组的解。syms