第二章 matlab符号计算.docx
《第二章 matlab符号计算.docx》由会员分享，可在线阅读，更多相关《第二章 matlab符号计算.docx（31页珍藏版）》请在冰豆网上搜索。
第二章matlab符号计算
Matlab符号计算
MATLAB自产生之日起就在数值计算应用方面独占鳌头。
但仍有缺少符号运算的缺陷，1993年MathWorks公司从加拿大滑铁卢大学购入了Maple函数库的使用权，在此基础上开发了MATLAB自己的符号计算工具箱（SymbolicToolbox），且保留了与Maple的接口。
从此，MATLAB便集数值计算、符号计算和图形可视化三大基本功能于一体，成为在数值计算领域中功能强大、操作简单、最受用户喜爱的语言。
在MATLAB中实现符号计算功能主要有以下3种途径：
（1）通过调用MATLAB自己开发的各种功能函数进行常用符号运算。
这些功能主要包括符号矩阵的基本操作、符号矩阵的运算、符号微积分运算、符号线性方程求解、符号微分方程求解、特殊数学符号函数、符号函数图形等，这些内容将在本章中作详细的介绍。
对于众多喜爱和熟悉MATLAB的用户来说，这些操作十分简单，很容易学习和掌握。
（2）MATLAB语言中的符号计算功能已径很强大了，但为了一些特殊专业人员提供方便，MATLAB中还保留着与Maple的接口，以实现更多功能。
在应用过程中需了解一些Maple的操作语法。
本章中将对此作简单介绍。
（3）对那些习惯了计算器的人来说，MATLAB同样是最佳的选择，因为MATLAB还提供了符号运算计算器（Functioncalculator）功能。
§2.1符号表达式的生成
MATLAB将数值计算中的变量认为是被赋值的数值变量，将符号运算中数字也视为符号。
符号表达式包括符号函数和符号方程，两者的区别是：
前者不包括等号，而后者必须带有符号。
但这两者的创建方式是相同的，并且和MATLAB中的字符串变量的生成方法相同。
（1）创建符号函数
>>f='log（x）'
f=
log（x）
（2）创建符号方程
>>eqation='a*x^2+b*x+c=0'
eqation=
a*x^2+b*x+c=0
（3）创建符号微分方程
>>diffeq='Dy-y=x'
diffeq=
Dy-y=x
由此种方法创建的符号表达式对空格是很敏感的。
因此不要在字符间乱加空格符，
否则在其他地方调用此表达式的时候会出错。
由于符号表达式在MATLAB中被看成是非1阶的符号矩阵，因此它也可用sym命令来创建。
如：
>>f=sym（'sin（x）'）
f=
sin（x）
>>f=sym（'sin（x）^2=0'）
f=
sin（x）^2=0
另外一种创建符号函数的方法是用syms命令来创建。
如：
>>symsx%此法不能用来创建符号方程
>>f=sin（x）+cos（x）
f=
sin（x）+cos（x）
§2.2符号矩阵的生成
在MATLAB中创建符号矩阵的方法和创建数值矩阵的形式很相似，只不过要用到符号定义函数sym,下面介绍使用此函数创建符号函数的几种形式。
2.2.1使用sym函数直接创建符号矩阵
>>a=sym（'[1/s+x,sin（x）,cos（x）^2/（b+x）;9,exp（x^2+y^2）,log（tanh（y））]'）%注意’的使用及位置
a=
[1/s+x,sin（x）,cos（x）^2/（b+x）]
[9,exp（x^2+y^2）,log（tanh（y））]
2.2.2用创建子阵的方法创建符号矩阵
这种创建方式需用逗号将同行各元素分隔，列元素对不齐时要用空格调整。
>>a=['[1/s+xsin（x）cos（x）^2/（b+x）]';'[9exp（x^2+y^2）log（tanh（y））]']
a=
[1/s+xsin（x）cos（x）^2/（b+x）]
[9exp（x^2+y^2）log（tanh（y））]
>>a=['[1/s+x,sin（x）,cos（x）^2/（b+x）]';'[9,exp（x^2+y^2）,log（tanh（y））]']
a=
[1/s+x,sin（x）,cos（x）^2/（b+x）]
[9,exp（x^2+y^2）,log（tanh（y））]
>>b=[a;'[exp（-1）,3,x^3+y^9]']%？
报错
2.2.3将数值矩阵转化为符号矩阵
在MATLAB中，数值矩阵不能直接参与符号计算，必须先转化为符号矩阵。
注意不论数值矩阵的元素原先是用分数还是用浮点数表示的，转化后的符号矩阵都将以最接近的精确有理数给出。
如：
>>a=[2/3,sqrt
（2）,0.222;1.4,1/0.23,log（3）]
a=
0.66671.41420.2220
1.40004.34781.0986
>>b=sym（a）
b=
[2/3,sqrt
（2）,111/500]
[7/5,100/23,4947709893870346*2^（-52）]
2.2.4符号矩阵的索引和修改
MATLAB的矩阵索引和修改同数值矩阵的索引和修改完全相同，即用矩阵的坐标括号表达式实现。
如：
>>b（1,3）%对上矩阵
ans=
111/500
>>b（2,3）='log（9）'
b=
[2/3,sqrt
（2）,111/500]
[7/5,100/23,log（9）]
§2.3符号矩阵的基本运算
作为MATLAB语言核心的矩阵在符号运算中也发挥着重要作用，对符号矩阵本身的操作也显得十分重要，MATLAB为符号矩阵提供了与数值矩阵相应的操作方式和函数。
2.3.1四则运算
1．矩阵的加（+）、减（-）法
>>a=sym（'[1/x,1/（x+1）;1/（x+2）,1/（x+3）]'）
a=
[1/x,1/（x+1）]
[1/（x+2）,1/（x+3）]
>>b=sym（'[x,1;x+2,0]'）
b=
[x,1]
[x+2,0]
>>a+b
ans=
[1/x+x,1/（x+1）+1]
[1/（x+2）+x+2,1/（x+3）]
>>b-a
ans=
[x-1/x,1-1/（x+1）]
[x+2-1/（x+2）,-1/（x+3）]
2．矩阵的乘（*）、除（/、\）法
>>a*b%矩阵a与b的线性代数乘法
ans=
[1+1/（x+1）*（x+2）,1/x]
[1/（x+2）*x+1/（x+3）*（x+2）,1/（x+2）]
>>a/b%相当于a*inv（b）
ans=
[1/（x+1）,-（x^2-x-1）/（x^2+3*x+2）/x]
[1/（x+3）,-（x^2+x-3）/（x^3+7*x^2+16*x+12）]
>>a\b%相当于inv（a）*b
ans=
[-6*x-2*x^3-7*x^2,3/2*x^2+x+1/2*x^3]
[6+2*x^3+10*x^2+14*x,-1/2*x^3-2*x^2-3/2*x]
3．矩阵的转置（’、.’）
>>a=sym（'[1/x,1/（x+1）;1/（x+2）,1/（x+3）]'）
a=
[1/x,1/（x+1）]
[1/（x+2）,1/（x+3）]
>>a'%结果中用到求复数共轭函数conj（）
ans=
[1/conj（x）,1/（2+conj（x））]
[1/（1+conj（x））,1/（3+conj（x））]
>>a.'
ans=
[1/x,1/（x+2）]
[1/（x+1）,1/（x+3）]
2.3.2符号矩阵的其他基本运算
1．行列式
>>det（a）
ans=
2/x/（x+3）/（x+1）/（x+2）
2．符号矩阵的逆（运算）
>>inv（b）
ans=
[0,1/（x+2）]
[1,-x/（x+2）]
3．符号矩阵的秩（运算）
>>rank（a）
ans=
2
4．符号矩阵的幂运算
>>a^2
ans=
[1/x^2+1/（x+1）/（x+2）,1/x/（x+1）+1/（x+1）/（x+3）]
[1/（x+2）/x+1/（x+3）/（x+2）,1/（x+1）/（x+2）+1/（x+3）^2]
5．符号矩阵的指数运算
>>b
b=[x,1]
[x+2,0]
>>exp（b）
ans=
[exp（x）,exp
（1）]
[exp（x+2）,1]
>>c=sym（'[0,a;b,0]'）
c=
[0,a]
[b,0]
>>c=sym（'[0,a;b,0]'）
c=
[0,a]
[b,0]
>>expm（c）%[V,D]=EIG（X）andEXPM（X）=V*diag（exp（diag（D）））/V
ans=
[1/2*exp（（a*b）^（1/2））+1/2*exp（-（a*b）^（1/2））,1/2*a*（-exp（-（a*b）^（1/2））+exp（（a*b）^（1/2）））/（a*b）^（1/2）]
[1/2*b*（-exp（-（a*b）^（1/2））+exp（（a*b）^（1/2）））/（a*b）^（1/2）,1/2*exp（（a*b）^（1/2））+1/2*exp（-（a*b）^（1/2））]
§2.4矩阵分解
这是MATLAB符号矩阵操作的重要组成部分之一，其中包括特征值解、奇异值分解等，以及约当型、三角抽取和矩阵空间等。
2.4.1符号矩阵的特征值分解——由eig函数实现
>>a='a';b='b';c='c';%清除变量a、b、c——同Maple
>>d=sym（'[a,b,c;b,c,a;c,a,b]'）
>>[v,lambda]=eig（d）
v=
[1,-（a+（b^2-b*a-c*b-c*a+a^2+c^2）^（1/2）-b）/（a-c）,-（a-（b^2-b*a-c*b-c*a+a^2+c^2）^（1/2）-b）/（a-c）]
[1,-（b-（b^2-b*a-c*b-c*a+a^2+c^2）^（1/2）-c）/（a-c）,-（b+（b^2-b*a-c*b-c*a+a^2+c^2）^（1/2）-c）/（a-c）]
[1,1,1]
lambda=
[b+a+c,0,0]
[0,（b^2-b*a-c*b-c*a+a^2+c^2）^（1/2）,0]
[0,0,-（b^2-b*a-c*b-c*a+a^2+c^2）^（1/2）]
2.4.2符号矩阵的奇异值分解——由svd函数实现
>>symstreal
>>A=[0,1;-1,0];
>>E=expm（t*A）
E=
[cos（t）,sin（t）]
[-sin（t）,cos（t）]
>>sigma=svd（E）
sigma=
[（cos（t）^2+sin（t）^2）^（1/2）]
[（cos（t）^2+sin（t）^2）^（1/2）]
>>simplify（sigma）
ans=
[1]
[1]
2.4.3符号矩阵的约当标准型——由jordan函数计算得到
>>a=sym（'[112;013;002]'）
a=
[1,1,2]
[0,1,3]
[0,0,2]
>>[x,y]=jordan（a）
x=
[5,-5,-5]
[3,0,-5]
[1,0,0]
y=
[2,0,0]
[0,1,1]
[0,0,1]
2.4.4符号矩阵的三角抽取——由diag、tril、triu函数实现
diag——抽主对角线元素
tril、triu——分别抽左、右三角形
>>z=sym（'[x*yx^asin（y）;t^alog（y）b;yexp（t）x]'）
z=
[x*y,x^a,sin（y）]
[t^a,log（y）,b]
[y,exp（t）,x]
>>diag（z）
ans=
[x*y]
[log（y）]
[x]
>>tril（z,-1）
ans=
[0,0,0]
[t^a,0,0]
[y,exp（t）,0]
>>triu（z,1）
ans=
[0,x^a,sin（y）]
[0,0,b]
[0,0,0]
>>triu（z,1）
ans=
[0,x^a,sin（y）]
[0,0,b]
[0,0,0]
2.4.5符号矩阵的列空间——求解函数为colspace
>>a=sym（'[112;013;002]'）
a=
[1,1,2]
[0,1,3]
[0,0,2]
>>b=colspace（a）%求解出来的列向量构成原矩阵列空间的基
b=
[1,0,0]
[0,0,1]
[0,1,0]
>>size（b,2）%给出a的秩（rank）
2.4.6符号矩阵的零空间——函数为null
Z=null（A）得到由奇异分解所得的零空间（AB=0时，B的列向量组的基）的正交基
Z=null（A,’r’）得到零空间的有理基。
AZ为零。
如果A为一个整数的小矩阵，元素R为小整数的比。
如：
A=sym（'[10,-1,1;-1,10,-2;-2,20,-4]'）
>>x=null（a）
x=
[-8/19]
[1]
[99/19]
>>x=null（a,'r'）%这种格式对符号矩阵不成立
?
?
?
Errorusing==>null
Toomanyinputarguments.
>>A=[1:
3;1:
3;1:
3]%对于数值矩阵
>>rank（a）
ans=
1
>>null（A）
ans=
0.96360
-0.1482-0.8321
-0.22240.5547
>>Z=null（A,'r'）
Z=
-2-3
10
11
>>A=[1:
3;012;002]
A=
123
012
002
>>rank（a）
ans=
3
>>Z=null（A）
Z=
Emptymatrix:
3-by-0
§2.5符号矩阵的简化
符号工具箱中还提供了符号矩阵因式分解、展开、合并、简化及通分等符号操作作函数。
下面将一一进行介绍。
2.5.1因式分解——factor
格式：
factor（S）——当S为整数时，则最佳因式分解将被计算
factor（sym（‘N’））——当S>2^52时，要使用此格式
>>symsx
>>factor（x^9-1）%这两行等价于“factor（sym（'x^9-1'））”
ans=
（x-1）*（x^2+x+1）*（x^6+x^3+1）
>>factor（sym（'12345678901234567890'））%12345678901234567890-2^52=1.2e+019>0
ans=
（2）*（3）^2*（5）*（101）*（3803）*（3607）*（27961）*（3541）
练习：
分解A=[x^2-1,x^3-1;1-2*x^2+x^4,x^2*y^2-y^2]，得B
2.5.2符号矩阵的展开——expand
常用于多项式、三角函数、指数函数和对数函数中。
>>symsxy
>>expand（（x+1）^3）
ans=
x^3+3*x^2+3*x+1
>>expand（sin（x+y））
ans=
sin（x）*cos（y）+cos（x）*sin（y）
练习：
将上面练习的结果B展开得到C，并比较C与A
2.5.3同类项合并——collect
格式：
collect（S,v）——按S中变量v的按降幂合并系数
collect（S）——对由findsym函数返回的缺省变量进行同类项合并
>>symsxy;
>>collect（x^2*y+y*x-x^2-2*x）
ans=
（y-1）*x^2+（y-2）*x
练习：
接着上面的练习完成：
collect（B,x）,collect（B,y）,collect（B）。
对C也完成同样的工作，比较结果
2.5.4符号简化——simple,simplify
1.simple寻找符号矩阵或符号表达式的最简型
格式：
simple（S）——对表达式S尝试多种不同的算法简化，以显示S表达式的长度
最短的简化形式。
若S为一矩阵，则结果是全矩阵的最短型，而可能非每个元素的最短型。
[R,HOW]=simple（S）——返回的R为简化型，HOW为简化过程中使用的主要方法。
>>symsx
>>[R,HOW]=simple（cos（x）^2+sin（x）^2）
R=
1
HOW=
combine%
>>[R,HOW]=simple（2*cos（x）^2-sin（x）^2）
R=
3*cos（x）^2-1
HOW=
simplify
>>[R,HOW]=simple（cos（x）^2-sin（x）^2）
R=
cos（2*x）
HOW=
combine
>>[R,HOW]=simple（cos（x）+（-sin（x）^2）^（1/2））
R=
cos（x）+i*sin（x）
HOW=
radsimp%
练习：
[R,HOW]=simple（（x+1）*x*（x-1）），simple（x^3+3*x^2+3*x+1）
simple（3*acos（x））
2．simplify符号简化函数符号矩阵的每一个元素
格式：
simplify（S）——对简化函数符号矩阵的每一个元素
>>symsx
>>simplify（sin（x）^2+cos（x）^2）
ans=
1
>>symscalphabeta
>>simplify（exp（c*log（sqrt（alpha+beta））））
ans=
（alpha+beta）^（1/2*c）
2.5.5分式通分——numden
格式：
[N,D]=numden（A）——转换A的各元素为分子和分母都是整系数的最佳多项式型
>>symsxy
>>[n,d]=numden（x/y+y/x）
n=
x^2+y^2
d=
y*x
练习：
观察“>>numden（x/y+y/x）”的结果，并说明结果的含义。
2.5.6特殊表达式的“秦九昭型”——horner
格式:
horner（P）——将符号多项式转换成嵌套形式表达式
>>horner（x^3-6*x^2+11*x-6）
ans=
-6+（11+（-6+x）*x）*x
§2.6符号函数
2.6.1复合函数的运算——compose
格式：
compose（f,g）——f=f（x）,g=g（y）,返回f（g（y））.
compose（f,g,z）——返回自变量z的复合函数f（g（g））,z可以是原来的自变量也可不是.
compose（f,g,x,z）——返回f（g（z））且使得x为f的独立变量.如f=cos（x/t）,则compose（f,g,x,z）返回cos（g（z）/t）,而compose（f,g,t,z）返回cos（x/g（z））.
compose（f,g,x,y,z）——返回f（g（z））并使得x为f的独立变量,y为g的独立变量。
即f=cos（x/t）,g=sin（y/u）,compose（f,g,x,y,z）返回cos（sin（z/u/t）,而compose（f,g,x,u,z）返回cos（sin（y/z）/t.
设：
>>symsxyztu;
>>f=1/（1+x^2）;g=sin（y）;h=x^t;p=exp（-y/u）;
（1）f（x）,g（y）:
>>compose（f,g）
ans=
1/（1+sin（y）^2）
>>compose（f,g,x）
ans=
1/（1+sin（x）^2）
>>compose（h,g,x,z）
ans=
sin（z）^t
>>compose（g,h,x,z）
ans=
sin（y）
>>compose（g,h,t,z）
ans=
sin（y）
>>compose（h,p,x,y,z）
ans=
exp（-z/u）^t
>>compose（h,p,t,u,z）
ans=
x^exp（-y/z）
练习：
>>compose（g,f），compose（f,g,y），compose（f,g,t），compose（h,g,t,z），
compose（g,h,t,z）
2.6.2反函数的运算——finverse
g=finvers（f）——解出反函数后，仍用原自变量作为反函数的自变量
g=finvers（f,v）——返回的自变量为v,这里v为一符号，是表达式的向量变量。
用于当f中变量是多个情表
>>symsxy
>>f=x^3;
>>finverse（f）
Warning:
finverse（x^3）isnotunique.
>InC:
\MATLAB6P1\toolbox\symbolic\@sym\finverse.matline43
ans=
x^（1/3）
>>f=x^2+y;
>>finverse（f,y）%（视x为参数）解出f中的y后，仍用y作自变量
ans=
-x^2+y
>>finverse（f）%多元函数求反函数，不指明自变量时，求的是什么？
Warning:
finverse（x^2+y）isnotunique.
>InC:
\MATLAB6P1\toolbox\symbolic\@sym\finverse.matline43
ans=
（-y+x）^（1/2）
练习：
（1）f=x+y^2，>>finverse（f）；
（2）f=s+t^2,>>finverse（f）;
（3）f=t+s^2,>>finverse（f）.
看有没有明确的结果。
2.6.3符号函数的二维图形——ezplot,fplot
1．符号函数的简易绘图函数
ezplot（f）,范围：
约[-2pi,2pi]
ezplot（f,xmin,xmax）
ezplot（f,[xmin,xmax]）
Ex.1绘制误差函数的图形：
>>ezplot（'erf（x）'）
>>ezploterf（x）
注:
在示图窗口菜单选:
Edit/CopyFigure可实现图形的复制。
2．绘制函数图函数fplot
fplot（fun,lims）:
fun必须为M文件函数或对变量x可执行字符串，此字符串被送入运行字符串命令函数eval后被执行。
函数fun（x）必须返回
升级会员 