Transcript 2. 第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
本例表现符号数字总是被准确记录和运算，而数值数字并不
总能保证被准确存储，运算时会引进截断误差。
3. 基本符号变量： sin(3)uz  xz  3w  a5  0
符号参数（表达式中的参数）, “待解符号变量”或
“自由符号变量” （表达式中的自变量x，默认为
x）
定义格式：
2
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')
无逗号
iii) syms a b c
syms a b c flag
flag为参数属性：
positive----参数取正实数
real-----参数为实数
unreal-----参数为限定的复数
4. 自由符号变量
解题结果是“用符号参数构成的表达式表述自由符号
变量”。解题时自由符号变量可“人为指定”，也可
“默认地自动认定”：与小写字母 x 的ASII码距离最
小的变量。
symvar(expression)
列出表达式中所有基本符号变量
symvar(expression,n) 列出表达式中认定n个自由符号变量
expression是符号表达式或符号表达式矩阵，x是
首选自由符号变量，认定优先次序为x, y, w, z, v等
【例2.1-2】
用符号计算研究方程 sin(3)uz 2  vz  3w  a5  0 的解。
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】元符号表达式、衍生符号表达式定义，基
本符号变量、自由符号变量的机器认定。
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]
无kz
symvar(EXPR,10)
ans =[ x, y, t, d, c, b, a, X, Y ]
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 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 符号计算中的函数指令
符号计算中的函数分成三个层次：
1. 与数值类函数和指令对应的同名符号类函数和指令。
2. 约50个经典特殊函数（误差函数erf、贝塞尔函数
besselj、椭圆积分ElliptiK等），借助mfun调用，
用mfunlist可列出。
3. 数量较大的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
3
2
4
Mc =
[a,b;c,d]
Ms =
[ a, b]
Mn =
1
3
2
4
Mc =
[a,b;c,d]
Ms =
[ a, b]
[ c, d]
(2) 三个矩阵的大小不同
SizeMn=size(Mn),
SizeMn =
2 2
SizeMc=size(Mc),
SizeMc =
1 9
SizeMs=size(Ms)
SizeMs =
2 2
(3) 用class获得每种矩阵的类别
CMn=class(Mn),
CMc=class(Mc),
CMs=class(Ms)
CMn =
double
CMc =
char
CMs =
sym
(4) 用isa判断每种矩阵的类别
isa(Mn,‘double’),
ans =
1
isa(Mc,'char'),
ans =
1
isa(Ms,'sym')
ans =
1
(5) 利用whos观察内存变量的类别和其他属性
whos
Name
Size
Name
Size
Mc
1x9
ans
1x1
Mn
2x2
t
1x201
Ms
2x2
y
1x201
Bytes Class P26
Bytes Class
18 char array
8
double
32 double array
1608
double
312
sym object
1608 double
a=0:1:6;
theta=sym('30*pi/180*a')
alfa=sym('30*pi/180')*a
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]
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')
iii) syms a b c
syms a b c flag
flag为参数属性：
positive----参数取正实数
real-----参数为实数
unreal-----参数为限定的复数
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 clear
syms x
f=x^3+4.75*x+2.5;
g=x^2+x+5;
rf=solve(f,x)
rg=solve(g,x)
rf = -1/2
Warning: Explicit solution could
1/4 - (79^(1/2)*i)/4
not be found. > In solve at 98
(79^(1/2)*i)/4 + 1/4
rg =[ empty sym ]
evalin(symengine,'getprop(x)')
ans =C_
evalin(symengine,'anames(Properties)')
ans =x
syms x real
syms x clear
rfr=solve(f,x)
rg=solve(g,x)
rfr =-1/2
rg = - (19^(1/2)*i)/2 - 1/2
evalin(symengine,'getprop(x)')
(19^(1/2)*i)/2 - 1/2
ans =R_
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帮
助浏览器。
2.2 符号数字及表达式的操作
2.2.1双精度数字与符号数字之间的转换
1. 双精度数字向符号数字的转换
sym(num, ‘r’) % “有理分数”表达的符号数字 : p/q, n^(p/q)
sym(num)
% sym(num, ‘r’) 的简写形式
sym(num, ‘f’) %以数值 N *2e 表示“浮点数”，N, e为整
数
[2e  N  2( 52e) ]
sym(1/10,'f') is 3602879701896397/36028797018963968
sym(num, ‘e’) % “有理分数”表达+机器实际浮点表达误差
e
sym(3*pi/4,'e') is 3*pi/4 - 103*eps/249
sym(num, ‘d’)
%十进制小数近似表示，有效数字位数受
2.2.1双精度数字与符号数字之间的转换
1. 双精度数字向符号数字的转换
sym(num)
sym(‘num’)
% sym(num, ‘r’) 的简写形式
% ‘num’是字符串数字
sym(‘num’) 中‘num’ 用作符号计算函数的输入
时，体现其理论真值，
对sym(num)中的num用作符号计算函数的输入
时，体现其字面数子的双精度近似。
2.2.1双精度数字与符号数字之间的转换
2. 符号数字向双精度数字转换
double(num_sym)
% 把符号数字转换为双精度数字
注意：double(‘num’) 把字符串数字转换为字符的ASCII码
double(‘0.1’)
ans = 48 46
f=sym('10/3')
Df=double(f)
class(Df)
49
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)
reset(symengine)
sa=sym('1/3+sqrt(2)') digits(48)
sa5=vpa(sa,5)
sa =2^(1/2) + 1/3
sa48=vpa(sa)
digits
sa32 =1.747546895706428382135022057543
Digits = 32
sa5 =1.7475
format long
sa48
a=1/3+sqrt(2)
=1.74754689570642838213502205754303
sa_Plus_a=vpa(sa+a,20) 141190300520871
sa_Minus_a=vpa(sa-a,20)
a=
1.747546895706428
sa_Plus_a =
3.4950937914128567869
sa_Minus_a =
-0.000000000000000022658064826339973669
2.2.3 符号表达式的基本操作
collect
expand
合并同类项
展开指定项
factor
horner
simple
因式分解
转换为嵌套形式 pretty
numden
simplify
获取最小公分母和相应分子
简化表达式，消除冗余项
搜索寻找最简表达
习惯方式显示
例2.2-2 简化
f 
syms x;
3 1
x3
 x62  12x  8
f=(1/x^3+6/x^2+12/x+8)^(1/3);
g1=simple(f)
g1 =((2*x + 1)^3/x^3)^(1/3)
验证:
f2=g1^3
f2 =(2*x + 1)^3/x^3
expand(f2)
ans =
12/x + 6/x^2 + 1/x^3 + 8
2.2.4 表达式中的置换操作
1.公共子表达式法简化表达
RS=subexpr(S);
RS=subexpr(S,’w’); [RS,w]=subexpr(S,’w’)
其中，w为MATLAB自动寻找的子表达式
【例2.2-3】 对符号矩阵
进行特征向量分解。
syms a b c d W;
A  a b 
 c d 
[V,D]=eig([a,b;c,d]);
[RVD,W]=subexpr([V;D],W)
V RVD
=
= [ (a/2 + d/2 - w/2)/c - d/c, (a/2 + d/2 + w/2)/c - d/c]
[(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,1,
1]
D=
[
a/2 + d/2 - w/2,
0]
[a/2 + d/2 - (a^2 - 2*a*d + d^2 + 4*b*c)^(1/2)/2,
0]
[
0,
a/2
+
d/2
+
w/2]
[
0, a/2 + d/2 + (a^2 - 2*a*d + d^2 + 4*b*c)^(1/2)/2]
w =(a^2 - 2*a*d + d^2 + 4*b*c)^(1/2)
2.2.4 表达式中的置换操作
2. 通用置换指令
• RES＝subs(ES,old,new)
• RES＝subs(ES,new)
例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置换规则示例
(5)
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)
5
4.5
4
3.5
3
2.5
f=a*sin(x)+b
k=(0.5:0.1:1)';
f6=subs(subs(f,{a,b},{k,2}),x,t);
size(f6)
plot(t,f6)
ans =
6 21
2
3
1.5
1
0
1
2
3
4
5
6
7
2.8
2.6
2.4
2.2
2
1.8
1.6
1.4
1.2
f6=2*sin(t)+3
1
0
1
2
3
4
5
6
7
Matlab6.5 版本上机问题
 单周上机，鞋套带上。
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)')
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版本错误,。
review
syms x flag para=sym('x', 'flag')
syms a positive; a=sym('a', 'positive')
sym(num) %双精度数字转符号数字sym(num, ‘r’) 的简写形式
sym(‘num’) %定义符号数字， ‘num’是字符串数字
double(num_sym)
% 把符号数字转换为双精度数字
digits(n) 设定 “十进制浮点”表示的有效数字位数；
xs=vpa(x,n) 据表达式x得到n位有效数字的符号数字xs
通用置换指令
• RES＝subs(ES,old,new)
• RES＝subs(ES,new)
2.3 符号微积分
2.3.1 极限和导数的符号计算
limit(f,v,a)
求极限
limit(f,v,a,’right’)
limit(f,v,a,’left’)
dfdvn=diff(f,v,n)
lim f (v )
v a
求右极限
lim f (v )
va
求左极限 lim f (v )
va
n
d
f
(
v
)
求
n
dv
fjac=jacobian(f,v)
求多元向量函数f(v)的jacobian矩阵
r=taylor(f,n,v,a)
把f(v)在v=a处进行泰勒展开
n 1

k 0
f
(k )
(a)
k
( x  a)
k!
1

【例2.3-1】试求 lim1  
x 
x

syms x k
kx
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.718281828459045235360287471352662
4977572470937
 a
t 
【例2.3-2】 f  
求
t cos x ln x 
syms a t x
3
df
dx
d2 f
dt 2
d2 f
dtdx
f=[a, t^3; t*cos(x), log(x)];
df=diff(f)
dfdt2=diff(f,t,2)
df =[
0, 0]
[-t*sin(x), 1/x]
dfdxdt=diff( diff(f,x) ,t)
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
2.3.2 序列/级数的符号求和
b
MATLAB求解通式求和  f (v)问题的指令为：
va
• s=symsum(f,v,a,b) 求通式f在指定变量v取遍
[a,b]中所有整数时的和。
说明：
（1）f是矩阵时，求和对逐个元素进行，但自变量定义在
整个矩阵上。
（2）v省缺时，f中的自变量由findsym（symvar）自动辨
认；b可以取有限整数也可以取无穷大。
（3）a，b同时省缺时，默认的自变量求和区间为[0,v-1]
n
1
【例2.3-8】求 
，
k 1 k (k  1)
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)]

1
( 1) k 

 ( 2k  1) 2 , k。
k 1 


2.3.3 符号积分
• intf=int(f,v)
给出f 对指定变量v的不定积分
• intf=int(f,v,a,b) 给出f对指定变量v的定积分
【例2.3-9】求
syms x
f= x*log(x)
s=int(f,x)
s=simple(s)
。
x
ln
xdx

f=
x*log(x)
s=
(x^2*(log(x) - 1/2))/2
s=
x^2*(log(x)/2 - 1/4)
2.3.3 符号积分
• intf=int(f,v) 给出f 对指定变量v的不定积分
• intf=int(f,v,a,b) 给出f对指定变量v的定积分
 ax
bx 2 
dx
【例2.3-10】求   1
。

sin x 


x

syms a b x;
The integral of f is
f2=[a*x, b*x^2; 1/x, sin(x)]
+-+
disp('The integral of f is')
|
2
3 |
| ax
bx |
pretty(int(f2))
f2 =
[ a*x, b*x^2]
[ 1/x, sin(x)]
| ---- , ---|
|
2
3
|
|
|
| log(x), -cos(x) |
+-+
例2.3-11求积分
2
x2
x2 y
x
xy
 
1
( x  y  z )dzdydx
2
2
2
syms x y z;
f=x^2+y^2+z^2;
F2=int(int(int(f, z, sqrt(x*y), x^2*y), y, sqrt(x), x^2), x, 1, 2)
第一重积分
第二重积分
Warning: Explicit integral could not be found.
F2 =(14912*2^(1/4))/4641 - (6072064*2^(1/2))/348075 +
(64*2^(3/4))/225 + 1610027357/6563700
VF2=vpa(F2)
%默认32位有效数字
VF2 =
224.92153573331143159790710032805
教材中48位有效数字

例2.3-12 求阿基米德螺线r=a*θ(a>0)在θ=0到 间的曲线长
度函数，并求出a=1,
 时的曲线长度。
2

2
2
L
(
x
'
)

(
y
'
)
  ( x ' ) 2  ( y ' ) 2 d
  ,ddl
x 0r cos , y  r sin


解: 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))
L =(phi*(a^2*phi^2 + a^2)^(1/2) + log(phi + (phi^2 + 1)^(1/2))*(a^2)^(1/2))/2
L_2pi=subs(L,[a,phi],sym('[1,2*pi]'))
L_2pi_vpa=vpa(L_2pi)
L_2pi_vpa =21.256294148209098800702512272566
x = cos( ) (log( + ( 2 + 1)1/2)/2 + ( ( 2 + 1)1/2)/2), y = sin( ) (log( + ( 2 + 1)1/2)/2 + ( ( 2 + 1)1/2)/2
例2.3-12 画阿基米德螺线r = a*θ及其曲线长度图形。
4
2
x  r cos , y  r sin 
0
-2
y
-4
-6
L =(phi*(a^2*phi^2 + a^2)^(1/2) + log(phi + (phi^2 +
1)^(1/2))*(a^2)^(1/2))/2
-8
-10
-12
-14
-16
-5
0
5
10
15
20
L1=subs(L,a,sym('1'));
ezplot(L1*cos(phi),L1*sin(phi),[0,2*pi])
grid on; hold on
x1=subs(x,a,sym('1'));
y1=subs(y,a,sym('1'));
h1=ezplot(x1,y1,[0,2*pi]);
set(h1,'Color','r','LineWidth',5)
title(' '); legend('螺线长度-幅角曲线','阿基米德螺线')
hold off
x
见P198
2.4 微分方程的符号解法
2.4.1 符号解法和数值解法有互补作用
从数值计算角度看，微分方程边值问题的求解比初值问题复杂困
难；而对符号计算来说，边值问题和初值问题的求解微分方程的
指令形式相同且简单，不过计算时间较长且未必有封闭形式解。
2.4.2 求微分方程符号解的一般指令
S=dslove(‘eq1,eq2,…,eqn’, ’cond1,cond2,…,condn’,’v’)
输入量为字符串
必选
可选
默认独立变量t
若应变量为y ，用“Dny”表示“y的n
条件：y(a)=b, Dy(c)=d
阶导数”， Dy为一阶导数。解在S.y中。
S=dslove(‘eq1’,’eq2’,…,’eqn’,’cond1’,‘cond2’,
说明：
（1）输入量包括三部分：微分方程、初始条件、指
定独立变量（不指定时，默认为t）。输入量必须
以字符串的形式给出。
（2）微分方程的记述规定：当y是应变量时，
用“Dny”表示“y的n阶导数”。
（3）关于初始条件或边界条件的规定：应写成
y(a)=b, Dy(c)=d等。
（4）如果y是应变量，则它的解在S.y中。
2.4.3 微分方程符号解示例
dx
dy
【例2.4-1】求
 y,
  x 的解。
dt
dt
clear
S=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通解和奇解的关系
clear all
y=dsolve('(Dy)^2-x*Dy+y=0','x')
y = x^2/4
C3*x - C3^2
clf, hold on
hy1=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)
end
hold off;box on
legend('奇解','通解','Location','Best')
ylabel('y')
title(['\fontsize{14}微分方程',' (y '')^2 – xy '' + y = 0 ','的解'])
见P199
例2.4-3.求解两点边值问题:
xy’’-3y’=x2,y(1)=0,y(5)=0
y=dsolve('x*D2y-3*Dy=x^2','y(1)=0,y(5)=0','x')
y =(31*x^4)/468 - x^3/3 + 125/468
ezplot(y,[-1,6])
hold on
plot([1,5],[0,0],'.r','MarkerSize',20)
text(1,1,'y(1)=0')
text(4,1,'y(5)=0')
title(['x*D2y-3*Dy=x^2',', y(1)=0,y(5)=0'])
hold off
2.5 符号变换和符号卷积
2.5.1 Fourier变换及其反变换
• Fw=fourier(ft,t,w) 求“时域”函数ft的Fourier变
换
• ft=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 符号卷积
利用拉氏变换求取。
（1）求Fourier变换
syms t w
ut=heaviside(t);
UT=fourier(ut)
UT = pi*dirac(w)-i/w
（2）求Fourier反变换
Ut=ifourier(UT,w,t)
Ut =heaviside(t)
（3）单位阶跃曲线
t=-2:0.01:2;
ut=heaviside(t);
kk=find(t==0);
Heaviside(t)
的Fourier变换。
1
t0
plot(t(kk),ut(kk),'.r','MarkerSize',30)
ut(kk)=NaN;
hold on
plot(t,ut,'-r','LineWidth',3)
plot([t(kk),t(kk)],[ut(kk-1),
ut(kk+1)],'or','MarkerSize',10)
hold off
t
grid on
axis([-2,2,-0.2,1.2])
xlabel('\fontsize{14}t'),ylabel('\fontsize{
14}ut')
title('\fontsize{14}Heaviside(t)')
0.8
0.6
ut
【例 2.5-1】求
t0
1
f (t )  
0
0.4
0.2
0
-0.2
-2
-1.5
-1
-0.5
0
0.5
1
1.5
2
【例2.5-4】求  at(t  a) u(t2 bt ) 的Laplace变换。
e sin bt
t e

>> syms t s; syms a b positive; %不限定无结果
>> Dt=dirac(t-a);
>> Ut=heaviside(t-b);
>> Mt=[Dt,Ut;exp(-a*t)*sin(b*t),t^2*exp(-t)];
>> MS=laplace(Mt,t,s)
MS =
[
exp(-s*a),
[ 1/b/((s+a)^2/b^2+1),
exp(-s*b)/s]
2/(s+1)^3]
review
2.3.1 极限和导数的符号计算
limit(f,v,a)
求极限
limit(f,v,a,’right’)
limit(f,v,a,’left’)
dfdvn=diff(f,v,n)
r=taylor(f,n,v,a)
lim f (v )
v a
求右极限
lim f (v )
va
求左极限 lim f (v )
va
n
d
f
(
v
)
求
dvn
把f(v)在v=a处进行泰勒展开
2.3.2 序列/级数的符号求和
b
 f (v )
va
s=symsum(f,v,a,b)
review
2.3.3 符号积分
• intf=int(f,v)
给出f 对指定变量v的不定积分
• intf=int(f,v,a,b) 给出f对指定变量v的定积分
2.4.2 求微分方程符号解的一般指令
S=dslove(‘eq1,eq2,…,eqn’, ’cond1,cond2,…,condn’,’v’)
输入量为字符串
必选
可选
默认独立变量t
若应变量为y ，用“Dny”表示“y的n
条件：y(a)=b, Dy(c)=d
阶导数”， Dy为一阶导数。解在S.y中。
S=dslove(‘eq1’,’eq2’,…,’eqn’,’cond1’,‘cond2’,
例2.4-2.图示微分方程y=xy'-(y')2通解和奇解的关系
clear all
y=dsolve('(Dy)^2-x*Dy+y=0','x')
clf, hold on
hy1=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)
end
hold off;box on
legend(‘奇解’,‘通解’,‘Location’,‘Best’) %参数的
设置有多种
ylabel('y')
title(['\fontsize{14}微分方程',' (y '')^2 – xy '' + y = 0 ','的解'])
legend('奇解','通解','Location','northwest')
docsearch legend Location %通过帮助指令获取参数
legend(...,'Location','location')
location can be either a 1-by-4 position vector ([left bottom width
height]) or one of the following strings.
Specifier
Location in Axes
North
Inside plot box near top
South
Inside bottom
East
Inside right
West
Inside left
NorthEast
Inside top right (default for 2-D plots)
NorthWest
Inside top left
SouthEast
Inside bottom right
SouthWest
Inside bottom left
NorthOutside
Outside plot box near top
SouthOutside
Outside bottom
EastOutside
Outside right
WestOutside
Outside left
NorthEastOutside
Outside top right (default for 3-D plots)
NorthWestOutside
Outside top left
SouthEastOutside
Outside bottom right
SouthWestOutside
Outside bottom left
Best
Least conflict with data in plot
见P199
BestOutside
Least unused space outside plot
2.6 符号矩阵分析和代数方程的解
2.6.1 符号矩阵分析
最常用的矩阵分析指令如下：
• det(A)
• diag(A)
• [V,D]=eig(A)
• expm(A) 矩阵指数函数e^A
• inv(A)
• [V,J]=jordan(A) Jordan分解，使AV=VJ
• poly(A)
• rank(A)
2.6.1 符号矩阵分析
a11 a12 

【例2.6-1】求矩阵 A   a a 的行列式、逆和特征根。
 21
22 

syms a11 a12 a21 a22
A=[a11, a12; a21, a22]
DA=det(A)
IA=inv(A)
A =[ a11, a12]
[ a21, a22]
DA =a11*a22 - a12*a21
IA =[ a22/(a11*a22 - a12*a21), -a12/(a11*a22 - a12*a21)]
[ -a21/(a11*a22 - a12*a21), a11/(a11*a22 - a12*a21)]
EA=subexpr(eig(A),'D')
D = (a11^2 - 2*a11*a22 + a22^2 + 4*a12*a21)^(1/2)
EA = a11/2 + a22/2 - D/2
a11/2 + a22/2 + D/2
 sin t对矩阵

cos t 
X  x cos t  y sin t
1/ 2 
(两点)的旋转作用.
Y  x sin t  y cos t
3 / 2
例2.6-2 Givens旋转
 3/2
A
 1/ 2
G  cos t
 sin t

解: syms t; A=sym([sqrt(3)/2,1/2;1/2,sqrt(3)/2]); %符号变量
G=[cos(t),-sin(t);sin(t),cos(t)]; GA=G*A
An=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 on
hu=plot(u1(:,1),u1(:,2),'--k',u2(:,1),u2(:,2),'b')
set(hu,'LineWidth',4)
Lstr=['旋转前的 v1';'旋转前的 v2';'旋转后的 u1';'旋转后的 u2'];
legend(Lstr,'Location','South')
axis([-1,1,-1,1]), axis square, hold off, grid on
2.6.2 线性方程组的符号解（参考4.2节）
矩阵计算是求解线性方程组最简单有效的方法,在
matlab中,不管数据对象是数值还是符号,实现矩阵
运算的指令形式几乎完全相同.  d  n  p  q
【例2.6-3】求线性方程组的解。

2 2
n  d  q  p  10

n
 qd   p
4

 q  p - n - 8d  1
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= 1
8
x=A\b
8
9
A=[1,1/2,1/2 , -1; 1, 1, -1,1; 1, -1/4, -1,1; -8, -1, 1, 1];
X=1.0000
b=[0; 10; 0; 1];
8.0000
8.0000
x=A\b
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) (可用)
只能是符号表达式 符号变量  d  n  p  q
【例2.6-3】求线性方程组的解。
syms d n p q
eq1= d+n/2+p/2-q ; eq2= d+n-p+q-10 ;
eq3= d-n/4-p+q ; eq4= -8*d-n+p+q-1 ;
S=solve(eq1,eq2,eq3,eq4, d , n , p , q );
disp([' d',' n',' p',' q'])
disp([S.d,S.n,S.p,S.q])

2 2
n  d  q  p  10

n
 qd   p
4

 q  p - n - 8d  1
d n pq
[ 1, 8, 8, 9]
2.6.3 一般代数方程组的解
【例2.6-3】求线性方程组的解。
eq1=sym('d+n/2+p/2-q');
eq2=sym('d+n-p+q-10');
eq3=sym('d-n/4-p+q');
eq4=sym('-8*d-n+p+q-1');
S=solve(eq1,eq2,eq3,eq4,'d','n','p','q');
disp([' d',' n',' p',' q'])
disp([S.d,S.n,S.p,S.q])
n p

 q

d

2 2
n  d  q  p  10

n
 qd   p
4

 q  p - n - 8d  1
d n pq
[ 1, 8, 8, 9]

uy
^
2

vz

w

0
【例2.6-4】求方程组
关于y, z 的解。
yzw0
S=solve('u*y^2+v*z+w=0','y+z+w=0','y','z')
disp('S.y'), disp(S.y), disp('S.z'), disp(S.z)
S=
y: [2x1 sym]
z: [2x1 sym]
S.y
(v + 2*u*w + (v^2 + 4*u*w*v - 4*u*w)^(1/2))/(2*u) - w
(v + 2*u*w - (v^2 + 4*u*w*v - 4*u*w)^(1/2))/(2*u) - w
S.z
-(v + 2*u*w + (v^2 + 4*u*w*v - 4*u*w)^(1/2))/(2*u)
-(v + 2*u*w - (v^2 + 4*u*w*v - 4*u*w)^(1/2))/(2*u)
不能是等式
其他形式：
S=solve('u*y^2+v*z+w=0','y+z+w=0', 'z', 'y')
S=solve('u*y^2+v*z+w','y+z+w', 'z,y')
syms y z u v w; S=solve(u*y^2+v*z+w,y+z+w, z,y)
[y, z]=solve(‘u*y^2+v*z+w’,‘y+z+w’, ‘z,y’) %按字母表输出顺序不变
n p

例2.6-5 求欠定方程的解。  d    q
2 2

n  d  q  p  10
syms d n p q
 qd  n  p
eq1=d+n/2+p/2-q;

4
eq2=n+d+q-p-10;
eq3=q+d-n/4-p;
S=solve(eq1,eq2,eq3,d,n,p,q)
disp([' S.d',' S.n',' S.p',' S.q'])
disp([S.d,S.n,S.p,S.q])
S = d: [1x1 sym]
n: [1x1 sym]
p: [1x1 sym]
q: [1x1 sym]
S.d S.n
S.p
S.q
[z/3 - 2, 8, (4*z)/3 - 4, z ]
n p
n
例2.6-5 求 d    q, n  d  q  p  10, q  d   p
4
欠定方程的解 2 2
解: syms d n p q
eq1=d+n/2+p/2-q;
eq2=n+d+q-p-10;
eq3=q+d-n/4-p;
S=solve(eq1,eq2,eq3)
disp([' S.n',' S.p','
disp([S.n,S.p,S.q])
S = n: [1x1 sym]
p: [1x1 sym]
q: [1x1 sym]
S.n S.p S.q
[8, 4*d+4, 3*d+6]
S.q'])
说明:解中自由变量选择规则:方程中所有变量按排序
最前的变量作为自由变量.
例2.6-6 求(x+2)x=2的解
syms x;
s=solve('(x+2)^x=2', 'x')
xs=(s(1)+2)^s(1)
s=
0.69829942170241042826920133106081
xs =
2.0
s 是MuPAD格式的矩阵，尽管只有一个元素。在
MuPAD中，标量不能与矩阵相加，并且矩阵的指数矩
阵没有定义。所以xs=(s+2)^s 格式会出错。
在matlab中是
2.7.代数状态方程求符号传递函数
2.7.1 结构框图的代数状态方程解法
2.8 符号计算结果的可视化
 符号计算结果的可视化两条途径：
1. 利用符号表达式直接绘图
2. 转换成数值数据，利用数值绘图指令绘图。
2.8.1 直接可视化符号表达式
Maltlab有一组专门实现函数可视化的指令,例ezplot,
冠以"ez"字头的指令为函数绘图(ez意为简易Easy to)
常用简捷作图指令见表2.8-1
2.8 符号计算结果的可视化
2.8.1 直接可视化符号表达式(ez意为简易Easy to)
指令名
ezcontour
含
义
画等位线
可执行实例
ezcontour('cos(x+sin(y))-sin(y)')
ezcontourf
画填色等位线
ezcontourf('sin(x)*sin(y)')
ezmesh
ezmeshc
画网线图
ezmesh('exp(-s)*cos(t)','exp(-s)*sin(t)','t',[0,8,0,4*pi])
ezplot
ezplot3
ezpolar
画二维曲线
画三维曲线
画极坐标曲线
ezplot('1/y-log(y)+log(-1+y)+x-1')
ezsurf
画曲面图
ezsurf('(x+8)*((y)^2/((x+8)^2+(y)^4+eps))','circ')
ezsurfc
画带等位线的曲面图
ezsurfc('sin(x)*sin(y)')
画带等位线的网线图 ezmeshc('y/(1+x^2+y^2)',[-5,5,-2*pi,2*pi])
ezplot3('sin(3*t)*cos(t)','sin(3*t)*sin(t)','t','animate')
ezpolar('sin(tan(t))')
2.8.1 直接可视化符号表达式
1. 单独立变量符号函数的可视化
ezplot(Fx) — 在x=[-2 2]*pi内绘制f(x)的函数图
ezplot(Fx,[xmin,xmax,ymin,ymax])
— 给定区间绘y=f(x)图形
ezplot(Fxy,[xmin,xmax,ymin,ymax],fig)
— 指定绘图窗口给定区间绘f(x,y)=0图形
ezplot(xt,yt,[tmin,tmax])
— 给定区间绘x=x(t),y=y(t)图形
ezplot(xt,yt, zt,[tmin,tmax])
— 给定区间绘x=x(t),y=y(t),z=x(t)三维曲形
Fx, Fxy,xt,yt,zt可以是字符表达式, 符号函数, 内联函
数,匿名函数句柄,函数M文件句柄.
注意:
plot(x, y)
%x,y为数值型数组
(2 cos((31/2 t)/2))/(3 exp(t/2))
例2.8-1：
0.6
0.4
2
绘制y  e
3

t
2
t
3
cos
t和它的积分s(t )   y (t )dt在[0,4 ]的图形.
2
0
0.2
0
-0.2
0
2
4
6
8
10
12
8
10
12
t
syms t tao;
s =  y(t)dt
0.5
y=2/3*exp(-t/2)*cos(sqrt(3)/2*t);
0.4
sint=int(y,t,0,tao);
s=subs(sint,tao,t)
0.3
0.2
0
2
4
subplot(2,1,1),ezplot(y,[0,4*pi]);grid on
subplot(2,1, 2),ezplot(s,[0,4*pi]);grid on
title('s= \int y(t)dt')
见P199
6
t
2. 双独立变量符号函数的可视化
ezsurf(Fxy) 绘制由表达式f(x,y)定义的表面图，自变量x
和y的变化范围为[-2π,2π]；
ezsurf(Fxy,domain) 绘制由表达式f(x,y)定义的表面图，
自变量 x和 y的变化范围由domain确定，domain可以
是[xmin,xmax,ymin,ymax]，也可是[min,max]，此时，
min<x<max,min<y<max。
ezsurf(x,y,z) 绘制由表达式x=x(s,t),y=y(s,t)和z=z(s,t)定义
的表面图,自变量s和t的变化范围均为[-2π,2π]
ezsurf(x,y,z,[smin,smax,tmin,tmax]) ；
ezsurf(…,ngrid) 绘表面图时按n×n的网格密度绘图，
ngrid 的缺省值为60；
ezsurf(…,’circ’) 以圆盘为自变量域绘制表面图，
2. 双独立变量符号函数的可视化
【例2.8-2】使用球坐标参量画部分球壳。
clf
x='cos(s)*cos(t)';
y='cos(s)*sin(t)';
z='sin(s)';
ezsurf(x,y,z,[0,pi/2,0,3*pi/2])
view(17,40)
shading interp
colormap(spring)
light('position',[0,0,-10],'style','local')
light('position',[-1,-0.5,2],'style','local')
material([0.5,0.5,0.5,10,0.3])
2.8.2 符号计算结果的数值化绘图
2
【例2.8-3】符号求函数 y  f ( x)  1 
x的积分，用
1 e
函数及其积分函数
plot指令绘制函数极其积分函数的曲线；对反函数积
分的两种算法进行可视化比较。
0.9
0.8
0.7
clear; syms x y real
fx=1-2/(1+exp(x));
fxint=int(fx,x,0,x)
0.6
0.5
fxint =
log((exp(x) + 1)^2/4) - x
0.4
f(x)
0.3
x0 f(x) dx
0.2
xk=0:0.1:2;
fxk=subs(fx,x,xk);
fxintk=subs(fxint,x,xk);
plot(xk,fxk,'g',xk,fxintk,'r','LineWidth',2.5)
title('函数极其积分函数' )
xlabel('x')
legend('f(x)','\int^x_0 f(x) dx','Location','best')
0.1
0
0
0.2
0.4
0.6
0.8
1
x
1.2
1.4
1.6
1.8
2
见P199
2
【例2.8-3】符号求函数 y  f ( x)  1 直 接 法 计 算 反的积分，对
x函 数 积 分
1

e
互补法求反函数积分
0.7
反函数积分的两种算法进行可视化比较。
clear; syms x y real;fx=1-2/(1+exp(x));
0.6
0.5
0.4
gy=subs(finverse(fx),x,y) % （3）求反函数
gyint=int(gy,y,0,y)
0.3
gy =log(-(y + 1)/(y - 1))
gyint =piecewise([y < 1, log(1 - y^2) +0.2 y*log(y + 1) - y*log(1 - y)], [1 <= y,
log(y^2 - 1) + y*log(-(y
+ 1)/(y - 1)) - pi*i])
0.1
gf=simplify(subs(gy,y,fx)) % （4）函数与反函数的互补性
0
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
gf =x
y
yk=subs(fx,x,xk); % （5）函数与反函数的积分的互补性
gyintk=subs(gyint,y,yk);
GYintk=xk.*fxk-fxintk; % fxk, fxintk是x取0:0.1:2时fx及其积分的值
plot(yk,gyintk,‘r’,‘LineWidth’,5); hold on
plot(yk,GYintk,'+k'); hold off; xlabel('y')
legend('直接法计算反函数积分', '互补法求反函数积分
','location','best')
2
【例2.8-3】符号求函数 y  f ( x)  1 
的积分，
x
1 e
对反函数积分的两种算法进行可视化比较。
画图验证函数积分与反函数积分的互补性
plot(xk,fxk,'-og',xk,fxintk,'r','LineWidth',2.5)
hold on
plot(xk,gyintk,'b','LineWidth',5)
hold on
plot(xk,xk.*fxk,'k','LineWidth',5)
hold on
plot(xk,fxintk+gyintk,'ok','LineWidth',5)
legend('f(x)','\int^x_0 f(x) dx','反函数积分', '函数积分+反函
数积分','矩形面积','location','best')
%分段函数的写法
syms y x
z=int(x*y,x,0,1);
g = evalin(symengine,[' piecewise([y > 1/2,' char(z) '], [y <= 1/2, 1])'])
G=int(g,y)
g=
piecewise([1/2 < y, y/2], [y <= 1/2, 1])
G=
piecewise([1/2 < y, y^2/4], [y <= 1/2, y])
习题2 (Page114)
3,4,5,6,7,8,11,12, 23,25题