ch2符号计算.docx
《ch2符号计算.docx》由会员分享，可在线阅读，更多相关《ch2符号计算.docx（26页珍藏版）》请在冰豆网上搜索。
ch2符号计算
第2章符号计算
符号计算：
解算数学表达式、方程不是在离散化的数值点上进行，而是凭借一系列恒等式，数学定理，通过推理和演绎，获得解析结果。
特点：
一，相对于MATLAB的数值计算“引擎”和“函数库”而言，符号计算的“引擎”和“函数库”是独立的。
二，在相当一些场合，符号计算解算问题的指令和过程，显得比数值计算更自然、更简明。
三，大多数理工科的本科学生在学过高等数学和其他专业基础课以后，比较习惯符号计算的解题理念和模式。
SymbolicMathToolbox™ThecomputationalengineunderlyingthetoolboxesisthekernelofMaplesoftware,asystemdevelopedprimarilyattheUniversityofWaterloo,Canadaand,morerecently,attheEidgen?
ssicheTechnischeHochschule,Zürich,Switzerland.MaplesoftwareismarketedandsupportedbyWaterlooMaple,Inc.
运算引擎MuPAD
MuPAD作为MATLAB7.8的符号计算工具箱，是一具有人工智能的数学软件。
方程式可以处理复数计算，完美的绘图功能，图型输入，输出,可以输入多个2-D函数或极坐标函数或3-D函数,选择所要绘图参数，就可以完成图形，以及图形的动画制作也是非常方便。
数值计算结果并不是MATLAB命令行窗口所得的类似代码形式，而是规范数学格式。
并拥有一内建的程序语言，帮助文档以及文本操作，文本操作在一定程度上可以取代word。
MathWorks自从2008年10开始，在Matlab的新版本（Matlab2008a，即7.6之后）中使用MuPAD内核替换原来的Maple符号计算内核！
.1符号对象和符号表达式
MATLAB依靠基本符号对象（包括数字、参数、变量）、运算符及一些预定义函数来构造和衍生符号表达式和符号方程。
.1.1符号对象的创建和衍生
10一生成符号对象的基本规则
●任何基本符号对象都必须借助专门的符号函数指令sym或syms定义。
●任何包含符号对象的表达式或方程，将继承符号对象的属性。
10二符号数字
符号（类）数字的定义：
sym（'Num'）创建一个符号数字Num
sc=sym（'Num'）创建一个符号常数sc，该常数值准确等于Num
说明：
Num代表一个具体的数字
Num必须处于（英文状态下的）单引号内，构成字符串（关于字符串参见附录A.1）。
【例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+sqrt（5）
Ca=
double
Csa=
sym
ans=
.138223758410852e-16
10三符号参数
表达式e-axsinbx中的a,b称为参数。
定义格式：
symsPara定义符号参数Para
Para=sym（'Para'）
symsParaFlag定义具有Flag指定属性的符号参数Para
Para=sym（'Para','Flag'）
symsPara1Para2ParaN定义Para1Para2ParaN为符号参数
symsPara1Para2ParaNFlag定义Para1Para2ParaN为具有Flag指定属性的符号参数
●符号参数名不要用处于“字母表中小写字母x及其两侧的英文字母”开头。
●Flag表示参数属性，可具体取以下词条：
positive表示那些符号参数取正实数；
real表示那些符号参数限定为实时；
unreal表示那些符号参数为不限定的复数。
symsxab
int（1/（x）,a,b）
Warning:
Warning,unabletodetermineif0isbetweenaand
b;trytouseassumptionsorset_EnvAllSolutionstotrue
Warning:
Explicitintegralcouldnotbefound.
>Insym.intat58
ans=
int（1/x,x=a..b）
Var=sym（'x'）;
Upp=sym（'a','real'）;
Low=sym（'b','real'）;
Intergral=int（1/（x）,a,b）
Warning:
Warning,unabletodetermineif0isbetweenaand
b;trytouseassumptionsorsetoption_EnvAllSolutionsto
true
Warning:
Explicitintegralcouldnotbefound.
>Insym.intat58
Intergral=
int（1/x,x=a..b）
Var=sym（'x'）;
Upp=sym（'a','positive'）;
Low=sym（'b','positive'）;
Intergral=int（1/（x）,a,b）
Intergral=
-log（a）+log（b）
10四符号变量
e-axsinbx中的x称为变量，符号变量的定义同符号参数。
确定自由符号变量的规则：
●在专门指定变量名的符号运算中，解题一定围绕指定变量名进行。
●自动识别符号变量时，字母的优先次序为x,y,w,z,v等。
自动识别表达式中自由、独立的符号变量的指令：
findsym（EXPR）确认表达式EXPR中所有自由符号变量
findsym（EXPR,N）确认表达式EXPR中距离x最近的N个自由符号变量
【例2.1-2】用符号计算研究方程
的解。
（1）不指定变量情况
symsuvwz%定义符号参数/变量
Eq=u*z^2+v*z+w;
result_1=solve（Eq）%
findsym（Eq,1）
result_1=
-u*z^2-v*z
ans=
w
（2）指定变量情况
result_2=solve（Eq,z）
result_2=
-1/2*（v-（v^2-4*u*w）^（1/2））/u
-1/2*（v+（v^2-4*u*w）^（1/2））/u
【例2.1-3】对独立自由符号变量的自动辨认。
（1）
symsabxXY%定义符号参数/变量
k=sym（'3'）;%符号常数
z=sym（'c*sqrt（delta）+y*sin（theta1）'）;%直接定义符号表达式
EXPR=a*z*X+（b*x^2+k）*Y;%构成衍生符号表达式
（2）
findsym（EXPR）
ans=
X,Y,a,b,c,delta,theta1,x,y
（3）
findsym（EXPR,1）
ans=
x
（4）
findsym（EXPR,2）,findsym（EXPR,9）
ans=
x,y
ans=
x,y,theta1,delta,c,b,a,X,Y
【例2.1-4】findsym确定自由变量是对整个矩阵进行的。
symsabtuvxy
A=[a+b*x,sin（t）+u;x*exp（-t）,log（y）+v]
findsym（A,5）
A=
[a+b*x,sin（t）+u]
[x*exp（-t）,log（y）+v]
ans=
x,y,v,u,t
.1.2符号计算中的算符
●与数值计算中的算符在形状、名称和使用方法上几乎完全相同。
●仅注意：
在符号对象的关系运算符中，只有算符“==”，“~=”
比较结果为“真”时，用1表示;否则用0表示。
.1.3符号计算中的函数指令
表2.1-1MATLAB中可调用的符号计算函数指令
类别
情况描述
与数值计算对应关系
基本函数
三角函数、双曲函数及反函数；除atan2外
名称和使用方法相同
指数、对数函数（如exp，expm）只有log，无log2和log10
symsx
log10（x）
名称和使用方法相同
复数函数（注意：
没有幅角函数angle）
z=1+i;
angle（z）
a=sym（'1+i'）;
abs（a）
angle（a）
名称和使用方法相同
矩阵分解函数（如eig等）
名称和使用方法相同
方程求解函数solve
不同
微积分函数（如diff，int）
不完全相同
积分变换和反变换函数（如laplace，ilaplace）
只有离散Fourier变换
绘图函数（如ezplot，ezsurf）
数值绘图指令更丰富
经典特殊函数
如误差函数erf、贝塞尔函数besselj、第一类完全椭圆积分EllipticK等；通过mfunlist可以看到所有经典函数名
部分
Maple库函数
Maple库函数在符号计算的扩展目录上；可通过mhelpindex看到各子函数库的名称；函数的数量很大；使用库函数，需要具备Maple语言知识
注意：
使用函数注意数据类型。
就数字而言，有双精度和符号类数字之分。
.1.4符号对象的识别
为了函数指令与数据对象的适配，MATLAB提供了用于识别数据对象属性的指令：
class（var）给出变量var的数据类别（如double，sym等）
isa（var,'Obj'）若变量var是Obj代表的类型，给出1，表示“真”
whos给出所有MATLAB内存变量的属性
【例2.1-5】数据对象及其识别指令的使用。
（1）
clear
a=1;b=2;c=3;d=4;%产生4个数值变量
Mn=[a,b;c,d]%利用已赋值变量构成数值矩阵
Mc='[a,b;c,d]'%字符串中的a,b,c,d与前面输入的数值变量无关
Ms=sym（Mc）%Ms是一个符号矩阵，它与前面各变量无关
Mn=
12
34
Mc=
[a,b;c,d]
Ms=
[a,b]
[c,d]
（2）
SizeMn=size（Mn）
SizeMc=size（Mc）
SizeMs=size（Ms）
SizeMn=
22
SizeMc=
19
SizeMs=
22
（3）
CMn=class（Mn）
CMc=class（Mc）
CMs=class（Ms）
CMn=
double
CMc=
char
CMs=
sym
（4）
isa（Mn,'double'）
isa（Mc,'char'）
isa（Ms,'sym'）
ans=
1
ans=
1
ans=
1
（5）
whosMnMcMs
NameSizeBytesClassAttributes
Mc1x918char
Mn2x232double
Ms2x2312sym
.2符号数字及表达式的操作
.2.1数值数字与符号数字之间的转换
10一数值数字向符号数字的转换
在符号运算中，“数值类数字”会自动地转换为符号数字。
亦可借助sym函数：
sym（Num,'r'）或sym（Num）数值类数字Num的广义有理表达
sym（Num,'d'）数值类数字Num的“十进制浮点”近似表达
sym（Num,'e'）数值类数字Num的带eps误差的理性近似表达
sym（Num,'f'）数值类数字Num的“十六进制浮点”近似表达
a=pi
a1=sym（a,'r'）
a2=sym（a,'d'）
a3=sym（a,'e'）
a4=sym（a,'f'）
a=
3.1416
a1=
pi
a2=
3.1415926535897931159979634685442
a3=
pi-198*eps/359
a4=
'1.921fb54442d18'*2^
（1）
a=pi+0.1;
a1=sym（a,'r'）
a2=sym（a,'d'）
a3=sym（a,'e'）
a4=sym（a,'f'）
a1=
7299417733396965*2^（-51）
a2=
3.2415926535897932048158054385567
a3=
7299417733396965*2^（-51）
a4=
'1.9eec82110f9e5'*2^
（1）
1、Num只能是数值类数字
2、sym（'Num'）和sym（Num）区别问题
'Num'数字字符串，理论值，Num表示数字，近似（双精度值）
3、在符号运算中，“数值数字”自动转换为符号运算
10二符号数字向双精度数字转换
double（Num_sym）把符号数字Num_sym转换为双精度数字
.2.2符号数字的任意精度计算
digits显示当前环境下符号数字“十进制浮点”表示的有效数字位数
digits（n）设定符号数字“十进制浮点”表示的有效数字位数（默认32位）
xs=vpa（x）据表达式x得到digits指定精度下的符号数字xs
xs=vpa（x,n）据表达式x得到n位有效数字的符号数字xs
【例2.2-1】digits,vpa指令的使用。
digits
p0=sym（'（1+sqrt（5））/2'）
pr=sym（（1+sqrt（5））/2）
pd=sym（（1+sqrt（5））/2,'d'）
e32r=vpa（abs（p0-pr））
e16=vpa（abs（p0-pd）,16）%16位变精度运算计算误差
e32d=vpa（abs（p0-pd））
e16=vpa（abs（p0-pd）,64）
digits
.2.3符号表达式的基本操作
collect（合并同类项）
factor（进行因式分解）
numden（提取公因式）等
最常用：
simple（EXPR）把EXPR转换成最简形式
【例2.2-2】简化
。
symsx
f=（1/x^3+6/x^2+12/x+8）^（1/3）;
g1=simple（f）
g2=simple（g1）
有时需要多次使用simple
simple（f）%直接看化简过程
simple（ans）
pretty（f）
.2.4表达式中的置换操作
10一子表达式置换操作
[RS,ssub]=subexpr（S,ssub）运用符号变量ssub置换子表达式，并重写S为RS
【例2.2-3】对符号矩阵
进行特征向量分解。
clearall
symsabcdW
[V,D]=eig（[ab;cd]）%V：
特征向量阵D：
特征值阵
[RVD,W]=subexpr（[V;D],W）%对矩阵元素中的公共子表达式进行置换表达
只置换较长的表达式
10二通用置换指令
RES=subs（ES,old,new）用new置换ES中的old后产生RES
RES=subs（ES,new）用new置换ES中的自由变量后产生RES
【例2.2-4】用简单算例演示subs的置换规则。
（1）产生符号函数
symsax;f=a*sin（x）+5
f=
a*sin（x）+5
（2）符号表达式置换
f1=subs（f,'sin（x）',sym（'y'））%<2>
class（f1）
（3）符号常数置换
f2=subs（f,{a,x},{2,sym（'pi/3'）}）%<3>a被双精度数字置换，x被符号数字置换
class（f2）
（4）双精度数值置换
f3=subs（f,{a,x},{2,pi/3}）%<4>
class（f3）
（5）数值数组置换之一
f4=subs（subs（f,a,2）,x,0:
pi/6:
pi）%<5>
class（f4）
（6）数值数组置换之二
f5=subs（f,{a,x},{0:
6,0:
pi/6:
pi}）%<6>
class（f5）
f6=subs（f,{a,x},{0:
7,0:
pi/6:
pi}）%<6>
class（f5）
.3符号微积分
.3.1极限和导数的符号计算
大学本科高等数学中的大多数微积分问题，都能用符号计算解决，手工笔算演绎的烦劳都可以由计算机完成。
limit（f,v,a）求极限
【例2.3-1】试求
。
symsxk
Lim_f=limit（（1-1/x）^（k*x）,x,inf）
diff（f,v,n）求
（n缺省时，默认n=1）
【例2.3-2】求
求
，
，
。
symsatx;f=[a,t^3;t*cos（x）,log（x）];
df=diff（f）%f对x的导数
dfdt2=diff（f,t,2）%f对x的二阶导数
dfdxdt=diff（diff（f,x）,t）%二阶混合导数
【例2.3-3】求
的Jacobian（雅可比）矩阵
。
symsx1x2;f=[x1*exp（x2）;x2;cos（x1）*sin（x2）];
v=[x1x2];fjac=jacobian（f,v）
【例2.3-4】
，求
，
。
（1）
clear
symsx
symsdpositive
f_p=sin（x）;%
df_p=limit（（subs（f_p,x,x+d）-f_p）/d,d,0）%<4>
df_p0=limit（（subs（f_p,x,d）-subs（f_p,x,0））/d,d,0）%<5>
（2）
f_n=sin（-x）;
df_n=limit（（f_n-subs（f_n,x,x-d））/d,d,0）%<7>
df_n0=limit（（subs（f_n,x,0）-subs（f_n,x,-d））/d,d,0）%<8>
（3）
f=sin（abs（x））;
dfdx=diff（f,x）%<10>
dfdx=
cos（abs（x））*abs（1,x）%abs（1,x）表示
一阶导数
（4）
xn=-3/2*pi:
pi/50:
0;xp=0:
pi/50:
3/2*pi;xnp=[xn,xp（2:
end）];
holdon
plot（xnp,subs（f,x,xnp）,'k','LineWidth',2.5）%<13>
plot（xn,subs（df_n,x,xn）,'--r','LineWidth',2.5）
plot（xp,subs（df_p,x,xp）,':
r','LineWidth',2.5）
legend（char（f）,char（df_n）,char（df_p）,'Location','NorthEast'）%<16>
gridon
xlabel（'x'）
holdoff
.3.2序列/级数的符号求和
symsum（f,v,a,b）求f在变量v取遍[a,b]中所有整数时的和。
a,b缺省时默认求和区间[0,v-1]。
【例2.3-8】求
，
。
symskt;f1=[tk^3];f2=[1/（2*k-1）^2,（-1）^k/k];
s1=simple（symsum（f1））%f1的自变量被确认为t
s2=simple（symsum（f2,1,inf））%f1的自变量被确认为k
symsxy;f1=[yx^3];f2=[1/（2*y-1）^2,（-1）^y/y];
s1=simple（symsum（f1））%f1的自变量被确认为x
.3.3符号积分
int（f,v）求f对变量v的不定积分
int（f,v,a,b）求f对变量v的定积分
【例2.3-9】求
。
clear
symsx
f=sqrt（（1+x）/x）/x
s=int（f,x）
s=simple（simple（s））
f=
（（x+1）/x）^（1/2）/x
s=
-2*（1/x+1）^（1/2）-2*atan（（1/x+1）^（1/2）*i）*i
s=
-2*（1/x+1）^（1/2）-2*atan（（1/x+1）^（1/2）*i）*i
【例2.3-10】求
。
symsabx;f=[a*x,b*x^2;1/x,sin（x）];
disp（'Theintegraloffis'）;pretty（int（f））
【例2.3-11】求积分
。
symsxyz
F2=int（int（int（x^2+y^2+z^2,z,sqrt（x*y）,x^2*y）,y,sqrt（x）,x^2）,x,1,2）
VF2=vpa（F2）%积分结果用32位数字表示
【例2.3-12】求阿基米德（Archimedes）螺线
在
到
间的曲线长度函数，并求出
时的曲线长度。
（1）
symsartheta1phi1positive
x=r*cos（theta1）;x=subs（x,r,a*theta1）;
y=r*sin（theta1）;y=subs（y,r,a*theta1）;
dLdth=sqrt（diff（x,theta1）^2+diff（y,theta1）^2）;
L=simple（int（dLdth,theta1,0,phi1））
（2）
L_2pi=subs（L,[a,phi1],sym（'[1,2*pi]'））
L_2pi_vpa=vpa（L_2pi）
L_2pi=
asinh（2*pi）/2+pi*（4*pi^2+1）^（1/2）
L_2pi_vpa=
21.256294148209098800702512272566
（3）
L1=subs（L,a,sym（'1'））;
ezplot（L1*cos（phi1）,L1*sin（phi1）,[0,2*pi]）
gr
升级会员 