第 2 章符号计算

第2章符号计算 教学目标教学重点教学内容

教学目标 • 一是讲述MATLAB符号计算基本知识，包括符号对象的创建、符号数字、符号表达式的操作； • 二是介绍符号微积分的计算； • 三是介绍符号矩阵分析和代数方程（组）的符号解法； • 四是介绍符号计算结果的可视化。

教学重点 • 熟悉符号对象的创建、符号数字、符号表达式的操作。 • 熟悉符号微积分的基本计算函数指令。 • 熟悉代数方程（组）的符号解法。 • 熟悉符号计算结果可视化的的基本指令。 • 了解符号计算帮助系统及其帮助指令。

教学内容 2.1符号对象和符号表达式 2.2 符号数字及表达式的操作 2.3 符号微积分 2.4 微分方程的符号解法 2.5 符号变换和符号卷积 2.6 符号矩阵分析和代数方程解 2.7 代数状态方程求符号传递函数 2.8 符号计算结果的可视化 2.9 符号计算资源深入利用

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

符号运算与数值运算的区别： 符号运算中，解算数学表达式、方程时，不是在离散化的数值点上进行，而是凭借一系列恒等式和数学定理，通过推理和演绎，获得解析结果。这种计算建立在数值完全准确表达和推演严格解析的基础上，所得结果是完全准确的。符号运算----代数运算，公式推导代值数值运算---算术运算

2.1 符号对象和符号表达式 在matlab中,数值和数值变量用于数值的存储和各种数值计算.而符号常量,符号变量,符号函数,符号操作等则是用来形成符号表达式,严格按照代数,微积分等课程中的规则,公式进行运算,并尽可能给出解析表达式. 2.1.1 符号对象的创建和衍生 ●数值计算---变量先赋值,再使用. ●符号计算---先定义基本的符号对象(可以是常量,变量,表达式),然后用这些基本符号对象去构成新的表达式,再进行所需的符号运算．

2.1.1 符号对象的创建和衍生 1. 生成符号对象的基本规则 ① 任何基本符号对象（数字、参数、变量、表达式）都必须借助专门的符号函数指令sym或syms定义。 ② 任何包含符号对象的表达式或方程，将继承符号对象的属性。即任何包含符号对象的表达式、方程也一定是符号对象。

2 符号数字的定义 格式：sc=sym('num')% sc为值为num的符号数字注意： i) 单引号必须在英文状态下输入，构成字符串 ii) num为一个具体的数字如： sc=sym('2/3') sb=sym('pi+sqrt(5)') sc=2/3 sb =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.3777 sa = pi + 5^(1/2) Ca = double Csa = sym ans =0.00000000000000001382237584108520004859354256418 本例表现符号数字总是被准确记录和运算，而数值数字并不总能保证被准确存储，运算时会引进截断误差。

无逗号 iii) syms a b c syms a b c flag 3. 基本符号变量：符号参数（表达式中的参数）, “待解符号变量”或“自由符号变量” （表达式中的自变量x，默认为x）定义格式： 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-----参数为限定的复数

4. 自由符号变量 • 解题结果是“用符号参数构成的表达式表述自由符号变量”。解题时自由符号变量可“人为指定”，也可“默认地自动认定”：与小写字母 x 的ASII码距离最小的变量。 symvar(expression) 列出表达式中所有基本符号变量 symvar(expression,n) 列出表达式中认定n个自由符号变量 expression是符号表达式或符号表达式矩阵，x是首选自由符号变量，认定优先次序为x, y, w, z, v等

【例2.1-2】 用符号计算研究方程的解。 syms u v w z a5 f=sym('3'); Eq=sin(f)*u*z^2+v*z+f*w-a5; symvar(Eq) %按字母表顺序列出基本符号变量, 无 f ans =[ a5, u, v, w, z] symvar(Eq,100) %按离x的距离列出所有自由符号变量 ans =[ w, z, v, u, a5] result_1=solve(Eq) result_1 =a5/3 - (v*z)/3 - (u*sin(3)*z^2)/3 result_2=solve(Eq,z) result_2 = -(v - (v^2 + 4*a5*u*sin(3) - 12*u*w*sin(3))^(1/2))/(2*u*sin(3)) -(v + (v^2 + 4*a5*u*sin(3) - 12*u*w*sin(3))^(1/2))/(2*u*sin(3))

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

4. 自由符号变量 • 【例2.1-4】symvar确定自由变量是对整个矩阵进行的。 syms a b t u v x y A=[a+b*x, sin(t)+u; x*exp(-t), log(y)+v] symvar(A,1) A = [ a + b*x, u + sin(t)] [ x/exp(t), v + log(y)] ans = x

2.1.2 符号计算中的算符 • 由于新版matlab采用了重载(Overload)技术, 使得用来构成符号计算表达式的算符和基本函数,无论在形式,名称,还是使用方法上,都与数值计算中的算符和基本函数几乎完全相同,这给编程带来极大的方便. (1) 基本运算符算符 ”+”, ”-”, ”*”, ”\”, “/”,“^” 分别构成矩阵的加, 减, 乘,左除,右除,求幂运算. 算符 ”.*”, “./”, “.\”, “.^” 分别实现元素对元素的数组乘,除,求幂运算.算符” ' ”, “ .' ” 分别实现矩阵的共轭转置,非共轭转置

2.1.2 符号计算中的算符 (2) 关系运算符在符号对象的比较中,没有大于,大于等于,小于,小于等于的概念,而只有是否等于的概念。 ”==“ “~=“分别用来对算符两边的对象进行相等和不等的比较,返回为逻辑量。事实为真时，比较结果为1，事实为假时，结果为0.

2.1.3 符号计算中的函数指令 符号计算中的函数分成三个层次： • 与数值类函数和指令对应的同名符号类函数和指令。 • 约50个经典特殊函数（误差函数erf、贝塞尔函数besselj、椭圆积分ElliptiK等），借助mfun调用，用mfunlist可列出。 • 数量较大的MuPAD库函数，借助evalin和feval调用。

与数值类函数和指令对应的同名符号类函数和指令与数值类函数和指令对应的同名符号类函数和指令 (1) 三角函数,双曲线函数以及他们的反函数除atan2只能用于数值计算外,另外的在两种运算中使用方法相同. (2) 指数,对数函数 sqrt, exp, expm在两者中用法相同．符号计算中只有自然对数log，而没有数值计算中的log2, log10 (3) 复数函数 conj, imag, real, abs在两者中用法相同.但在符号计算中没有求相角的指令angle.

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

2.1.4 符号对象的识别 • 数值计算对象,符号计算对象,字符串是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) % 符号变量，与前面的各变量无关 Mn = 1 2 3 4 Ms = [ a, b] [ c, d] Mc = [a,b;c,d]

Mn = 1 2 3 4 Mc = [a,b;c,d] Ms = [ a, b] [ c, d] (3) 用class获得每种矩阵的类别 CMn=class(Mn), CMc=class(Mc), CMs=class(Ms) CMn = double CMc = char CMs = sym (4) 用isa判断每种矩阵的类别 isa(Mn,‘double’), isa(Mc,'char'), isa(Ms,'sym') ans = 1 ans = 1 ans = 1 (2) 三个矩阵的大小不同 SizeMn=size(Mn), SizeMc=size(Mc), SizeMs=size(Ms) SizeMn = 2 2 SizeMs = 2 2 SizeMc = 1 9

theta = 30*pi/180*a a与数组a无关 alfa = [ 0, 1/6*pi, 1/3*pi, 1/2*pi, 2/3*pi, 5/6*pi, pi] (5) 利用whos观察内存变量的类别和其他属性 whos P26 Name Size Bytes Class ans 1x1 8 double t 1x201 1608 double y 1x201 1608 double Name Size Bytes Class Mc 1x9 18 char array Mn 2x2 32 double array Ms 2x2 312 sym object a=0:1:6; theta=sym('30*pi/180*a') alfa=sym('30*pi/180')*a

2.1.5 符号运算机理和变量假设 1.符号运算的工作机理 • Matlab借助sym或syms定义符号变量时，启动MuPAD引擎并启动一个专供MuPAD使用的内存工作空间执行符号运算； • matlab内存空间只保存该符号变量和计算结果。 • 该定义变量的限定性假设（assumption）被保存在MuPAD的内存空间。 • 若变量不带限定性假设，则MuPAD默认为复数。

2.1.5 符号运算机理和变量假设 2. 对符号变量的限定性假设 i) syms x  para=sym('x') ii) syms x flag para=sym('x', 'flag') syms a positive; a=sym('a', 'positive') flag为参数属性： positive----参数取正实数 real-----参数为实数 unreal-----参数为限定的复数 iii) syms a b c syms a b c flag

2.1.5 符号运算机理和变量假设 3. 清除变量和撤销假设 • 符号变量和其假设存放在不同的内存空间，因此删除符号变量和撤销关于变量的假设是两件需要分别处理的事情。 clear all 可同时删除 clear x 删除matlab内存中的x变量 syms x clear 撤销MuPAD内存中对变量x的假设 clear all删除matlab及 MuPAD内存中的所有变量 evalin(symengine,'getprop(x) ') 获取x的限定性假设 evalin(symengine,'anames (Properties)') 列出MuPAD内存中带限定性假设的符号变量 reset(symengine) 重启MuPAD引擎，清空MuPAD内存

【例2.1-6】syms 对符号变量限定性假设的影响 clear x syms x g=x^2+x+5; rg=solve(g,x) Warning: Explicit solution could not be found. > In solve at 98 rg =[ empty sym ] syms x clear f=x^3+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(x)') ans =C_ evalin(symengine,'anames(Properties)') ans =x syms x real rfr=solve(f,x) rfr =-1/2 evalin(symengine,'getprop(x)') ans =R_ syms x clear rg=solve(g,x) rg = - (19^(1/2)*i)/2 - 1/2 (19^(1/2)*i)/2 - 1/2 reset(symengine) clear all

2.1.6 符号帮助体系 • Matlab符号指令帮助系统与2008年前不同，因为用MuPAD取代maple作为符号计算的内核。 help SymName helpwin SymName doc SymName doc mfunlist doc(symengine)

2.1.6 符号帮助体系 • Matlab符号指令帮助系统与2008年前不同，因为用MuPAD取代maple作为符号计算的内核。 help SymName helpwin SymName docsearch laplace doc mfunlist doc(symengine)

2.1.6 符号帮助体系 • Matlab符号指令帮助系统与2008年前不同，因为用MuPAD取代maple作为符号计算的内核。 help SymName helpwin SymName docsearch laplace doc mfunlist doc(symengine) • 开启MuPAD帮助浏览器。

sym(num, ‘f’) %以数值 表示“浮点数”，N, e为整数 sym(1/10,'f') is 3602879701896397/36028797018963968 2.2 符号数字及表达式的操作 2.2.1双精度数字与符号数字之间的转换 1.双精度数字向符号数字的转换 sym(num, ‘r’) % “有理分数”表达的符号数字: p/q, n^(p/q) sym(num) % sym(num, ‘r’) 的简写形式 sym(num, ‘e’) % “有理分数”表达+机器实际浮点表达误差e sym(3*pi/4,'e') is 3*pi/4 - 103*eps/249 sym(num, ‘d’) %十进制小数近似表示，有效数字位数受 digits指令控制。默认为digits(32)情况.

2.2.1双精度数字与符号数字之间的转换 1.双精度数字向符号数字的转换 sym(num) % sym(num, ‘r’) 的简写形式 sym(‘num’) % ‘num’是字符串数字 • sym(‘num’)中‘num’ 用作符号计算函数的输入时，体现其理论真值， • 对sym(num)中的num用作符号计算函数的输入时，体现其字面数子的双精度近似。

2.2.1双精度数字与符号数字之间的转换 2.符号数字向双精度数字转换 double(num_sym) % 把符号数字转换为双精度数字注意：double(‘num’)把字符串数字转换为字符的ASCII码 double(‘0.1’) ans = 48 46 49 f=sym('10/3') Df=double(f) class(Df) Df = 3.3333 ans = double

2.2.2 符号数字的任意精度表达形式 • 数值计算与符号计算的最重要区别：数值计算存在截断误差，且在计算中不断传播，产生累计误差；符号计算完全准确，没有累计误差，但降低计算速度，增加内存。为兼顾精度和速度，采用“变精度”算法。 digits 显示当前环境下符号数字“十进制浮点” 表示的有效数字位数； digits(n) 设定 “十进制浮点”表示的有效数字位数； xs=vpa(x) 据表达式x得到digits指定精度下的符号数字xs xs=vpa(x,n) 据表达式x得到n位有效数字的符号数字xs

例2.2-1 digits, vpa, symengine指令演示 sa32=vpa(sa) digits(48) sa5=vpa(sa,5) sa48=vpa(sa) sa32 =1.747546895706428382135022057543 sa5 =1.7475 sa48 =1.74754689570642838213502205754303141190300520871 reset(symengine) sa=sym('1/3+sqrt(2)') sa =2^(1/2) + 1/3 digits Digits = 32 format long a=1/3+sqrt(2) sa_Plus_a=vpa(sa+a,20) sa_Minus_a=vpa(sa-a,20) a = 1.747546895706428 sa_Plus_a = 3.4950937914128567869 sa_Minus_a = -0.000000000000000022658064826339973669

例2.2-2 简化 2.2.3符号表达式的基本操作验证: f2=g1^3 f2 =(2*x + 1)^3/x^3 expand(f2) ans = 12/x + 6/x^2 + 1/x^3 + 8 syms x; f=(1/x^3+6/x^2+12/x+8)^(1/3); g1=simple(f) g1 =((2*x + 1)^3/x^3)^(1/3)

【例2.2-3】对符号矩阵 进行特征向量分解。 1.公共子表达式法简化表达 RS=subexpr(S); RS=subexpr(S,’w’); [RS,w]=subexpr(S,’w’) 其中，w为MATLAB自动寻找的子表达式 2.2.4 表达式中的置换操作 syms a b c d W; [V,D]=eig([a,b;c,d]); [RVD,W]=subexpr([V;D],W) 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/2] w =(a^2 - 2*a*d + d^2 + 4*b*c)^(1/2) V = [(a/2+d/2-(a^2-2*a*d+d^2+4*b*c)^(1/2)/2)/c-d/c, (a/2+d/2+(a^2-2*a*d+ d^2+4*b*c)^(1/2)/2)/c-d/c] [ 1, 1] D = [a/2 + d/2 - (a^2 - 2*a*d + d^2 + 4*b*c)^(1/2)/2, 0] [ 0, a/2 + d/2 + (a^2 - 2*a*d + d^2 + 4*b*c)^(1/2)/2]

2.通用置换指令 • RES＝subs(ES,old,new) • RES＝subs(ES,new) 2.2.4 表达式中的置换操作例2.2-4 subs置换规则示例 clear syms a b x; f=a*sin(x)+b f1=subs(f,sin(x),'log(y)') f2=subs(f,a,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))/100 f3 =3^(1/2) + 5 ans =sym

例2.2-4 subs置换规则示例 f=a*sin(x)+b (5) k=(0.5:0.1:1)'; f6=subs(subs(f,{a,b},{k,2}),x,t); size(f6) plot(t,f6) ans = 6 21 format %设为默认输出格式显示 format compact f4=subs(f,{a,b,x},{2,5,pi/3}) class(f4) f4 = 6.7321 ans =double f5=subs(f,x,0:pi/2:pi) class(f5) f5 =[ b, a + b, b] ans =sym t=0:pi/10:2*pi; f6=subs(f,{a,b,x},{2,3,t}) plot(t,f6) f6=2*sin(t)+3

单周上机，鞋套带上。 clear all; clear mex syms x unreal • 符号限定假设对方程根不起作用。 syms x clear f=x^3+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 syms x real rfr=solve(f,x) rfr =-1/2 evalin(symengine,'getprop(x)') ans =R_ evalin(symengine,'anames(Properties)') Matlab6.5 版本上机问题 findsym(expression,n) 当n大于实际的基本变量数目时，按字母表顺序列出所有本符号变量；当n小于等于时，按与x距离顺序列出。symvar(expression,n) 多一个参数n在Matlab6.5中不能用 • 调用word问题。指令的输入、编写用M文件编辑器。 • 习题2 (Page114)除了第3题不能得到实根、正根外，1~12题, 23~25题都可做，只是4题结果在差值有效数字上与2010版本有差别, 7题不能得到最简结果,9题结果 6.5版本错误,。

syms x flag para=sym('x', 'flag') syms a positive; a=sym('a', 'positive') sym(num) %双精度数字转符号数字sym(num, ‘r’) 的简写形式 sym(‘num’) %定义符号数字， ‘num’是字符串数字 review double(num_sym) % 把符号数字转换为双精度数字 digits(n)设定 “十进制浮点”表示的有效数字位数； xs=vpa(x,n)据表达式x得到n位有效数字的符号数字xs 通用置换指令 • RES＝subs(ES,old,new) • RES＝subs(ES,new)

limit(f,v,a) 求极限 limit(f,v,a,’right’) 求右极限 limit(f,v,a,’left’) 求左极限 2.3 符号微积分 2.3.1极限和导数的符号计算 dfdvn=diff(f,v,n) 求 fjac=jacobian(f,v) 求多元向量函数f(v)的jacobian矩阵 r=taylor(f,n,v,a) 把f(v)在v=a处进行泰勒展开

syms x k f=(1-1/x)^(k*x); Lf=limit(f,x,inf) Lf =1/exp(k) Lf1=vpa(subs(Lf,k,sym('-1')),48) Lf1=2.7182818284590452353602874713526624977572470937 【例2.3-1】试求

【例2.3-2】求 syms a t x f=[a, t^3; t*cos(x), log(x)]; df=diff(f) dfdt2=diff(f,t,2) dfdxdt=diff( diff(f,x) ,t) df =[ 0, 0] [-t*sin(x), 1/x] dfdt2 =[ 0, 6*t] [ 0, 0] dfdxdt =[ 0, 0] [ -sin(x), 0]

例2.3-5:设cos(x+siny)=siny, 求dy/dx(隐函数求导). syms x g=sym('cos(x+sin(y(x)))=sin(y(x))') dgdx=diff(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)))) %将dy/dx表达式用x,y表达

例2.3-6:求f(x)=xex在x=0处展开的8阶Maclaurin级数，即忽略9阶及以上小量的泰勒级数展开 。 sym x r=taylor(x*exp(x),9,x,0) pretty(r) r = x^8/5040 + x^7/720 + x^6/120 + x^5/24 + x^4/6 + x^3/2 + x^2 + 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 序列/级数的符号求和说明：（1）f是矩阵时，求和对逐个元素进行，但自变量定义在整个矩阵上。（2）v省缺时，f中的自变量由findsym（symvar）自动辨认；b可以取有限整数也可以取无穷大。（3）a，b同时省缺时，默认的自变量求和区间为[0,v-1]

【例2.3-8】求， 。 syms n k f1=1/(k*(k+1)); s1=symsum(f1,k,1,n) s1 =1 - 1/(n + 1) f2=[1/(2*k-1)^2, (-1)^k/k ] s2=symsum(f2,1,inf) s3 = [ pi^2/8, -log(2)]

intf=int(f,v) 给出f 对指定变量v的不定积分 • intf=int(f,v,a,b) 给出f对指定变量v的定积分 2.3.3 符号积分【例2.3-9】求。 syms x f= x*log(x) s=int(f,x) s=simple(s) f = x*log(x) s = (x^2*(log(x) - 1/2))/2 s = x^2*(log(x)/2 - 1/4)

【例2.3-10】求。 • intf=int(f,v) 给出f 对指定变量v的不定积分 • intf=int(f,v,a,b) 给出f对指定变量v的定积分 2.3.3 符号积分 syms a b x; f2=[a*x, b*x^2; 1/x, sin(x)] disp('The integral of f is') pretty(int(f2)) The integral of f is +- -+ | 2 3 | | a x b x | | ---- , ---- | | 2 3 | | | | log(x), -cos(x) | +- -+ f2 = [ a*x, b*x^2] [ 1/x, sin(x)]

第 2 章符号计算