MATLAB入门北航 张志涌版ch2符号计算.docx

MATLAB入门北航 张志涌版ch2符号计算.docx

《MATLAB入门北航 张志涌版ch2符号计算.docx》由会员分享，可在线阅读，更多相关《MATLAB入门北航 张志涌版ch2符号计算.docx（31页珍藏版）》请在冰豆网上搜索。

MATLAB入门北航张志涌版ch2符号计算

第2章符号计算

符号计算：

解算数学表达式、方程不是在离散化的数值点上进行，而是凭借一系列恒等式，数学定理，通过推理和演绎，获得解析结果。

特点：

一，相对于MATLAB的数值计算“引擎”和“函数库”而言，符号计算的“引擎”和“函数库”是独立的。

二，在相当一些场合，符号计算解算问题的指令和过程，显得比数值计算更自然、更简明。

三，大多数理工科的本科学生在学过高等数学和其他专业基础课以后，比较习惯符号计算的解题理念和模式。

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）。

●Num为整数、有理数、预定义常数（pi）

●如果定义十进制小数，可能会造成符号精度误差

formatlong

a=sym（'sin（0.3）'）;

b=sym（'sin（3e-1）'）;

c=sym（'sin（3/10）'）;

a==b

vpa（b-a,40）

【例2.1-1】符号（类）数字与数值（类）数字之间的差异。

a=pi+sqrt（5）%创建方式

sa=sym（'pi+sqrt（5）'）

Ca=class（a）%类别判断

Csa=class（sa）

vpa（sa-a）

10三基本符号变量

表达式e-axsinbx中的a,b称为参数，x为变量。

定义格式：

symsPara定义符号参数Para

Para=sym（'Para'）

symsParaFlag定义具有Flag指定属性的符号参数Para

Para=sym（'Para','Flag'）

symsPara1Para2ParaN定义Para1Para2ParaN为符号参数

symsPara1Para2ParaNFlag定义Para1Para2ParaN为具有Flag指定属性的符号参数

●符号参数名不要用处于“字母表中小写字母x及其两侧的英文字母”开头。

●Flag表示参数属性，可具体取以下词条：

positive表示那些符号参数取正实数；

real表示那些符号参数限定为实时；

unreal和clear作用相同，撤销定义，恢复复数域

默认值是复数域符号变量

sym指令只能对单变量作用，syms不能对数值、常数相关的定义。

symsxab

int（1/（x）,a,b）

Var=sym（'x','positive'）;

Upp=sym（'a','real'）;

Low=sym（'b','real'）;

Intergral=int（1/（x）,a,b）

Var=sym（'x','positive'）;

Upp=sym（'a','positive'）;

Low=sym（'b','positive'）;

Intergral=int（1/（x）,a,b）

10四符号变量

e-axsinbx中的x称为变量，符号变量的定义同符号参数。

确定自由符号变量的规则：

●在专门指定变量名的符号运算中，解题一定围绕指定变量名进行。

●自动识别符号变量时，字母的优先次序为x,y,w,z,v等。

自动识别表达式中自由、独立的符号变量的指令：

symvar（EXPR）确认表达式EXPR中所有自由符号变量

symvar（EXPR,N）确认表达式EXPR中距离x最近的N个自由符号变量

原来版本是

findsym（EXPR）确认表达式EXPR中所有自由符号变量

findsym（EXPR,N）确认表达式EXPR中距离x最近的N个自由符号变量

【例2.1-2】用符号计算研究方程

的解。

（1）不指定变量情况

symsuvwz%定义符号参数/变量

f=sym（'3'）;%f是符号常数

Eq=sin（f）*u*z^2+v*z+w;

result_1=solve（Eq）

symvar（Eq）%

symvar（Eq,100）

（2）指定变量情况

result_2=solve（Eq,z）

【例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）

symvar（EXPR）

（3）

symvar（EXPR,1）

（4）

symvar（EXPR,2）,findsym（EXPR,9）

【例2.1-4】symvar确定自由变量是对整个矩阵进行的。

symsabtuvxy

A=[a+b*x,sin（t）+u;x*exp（-t）,log（y）+v]

symvar（A,5）

.1.2符号计算中的算符

●与数值计算中的算符在形状、名称和使用方法上几乎完全相同。

●仅注意：

在符号对象的关系运算符中，只有算符“==”，“~=”

比较结果为“真”时，用1表示;否则用0表示。

a=sym（'4'）

b=sym（'5'）

a~=b

a

.1.3符号计算中的函数指令

表2.1-1MATLAB中可调用的符号计算函数指令

类别

情况描述

与数值计算对应关系

基本函数

三角函数、双曲函数及反函数；除atan2外

名称和使用方法相同

指数、对数函数（如exp，expm）只有log，无log2和log10

symsx

log2（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可以看到所有经典函数名

部分

MuPAD库函数

MuPAD库函数在符号计算的扩展目录上；可通过mhelpindex看到各子函数库的名称；函数的数量很大；使用库函数，需要具备MuPAD语言知识

无

注意：

使用函数注意数据类型。

就数字而言，有双精度和符号类数字之分。

plot

.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是一个符号矩阵，它与前面各变量无关

（2）

SizeMn=size（Mn）

SizeMc=size（Mc）

SizeMs=size（Ms）

（3）

CMn=class（Mn）

CMc=class（Mc）

CMs=class（Ms）

（4）

isa（Mn,'double'）

isa（Mc,'char'）

isa（Ms,'sym'）

（5）

whosMnMcMs

.1.5符号运算的工作机理

工作机理

●sym或syms启动MuPAD引擎，开启MuPAD单独的内存空间

●定义变量保存至Matlab内存空间

●对变量的限定性假设，保存在MuPAD空间中

对该定义变量的限定性假设

x=sym（'x','real'）

symsxreal

清除变量、消除假设

clearx

symsxclear

assumptions（x）

assumptions

reset（symengine）

【例2.1-7】

clearall

reset（symengine）

Da=1.2;Dw=1/3;

symssaswsxsysz

symsABpositive

（2）

whos

（3）

assumptions

（4）

clearA

who

assumptions

（5）

symsBclear

who

assumptions

.1.6符号帮助体系

1、不依赖mfun，evalin，feval指令调用

docsymname

doclaplace

helpsymname

helplaplace

2、借助mfun调用的特殊函数

docmfunlist

helpmfunlist

3、借助evalin，feval指令调用MuPAD指令

doc（symengine）

【例2.1-8】

（1）

图2.1-1

（2）

图2.1-2

●（3）

图2.1-3

symsabcx

p=a*x^2+b*x+c;

feval（symengine,'polylib:

:

discrim',p,x）

evalin（symengine,'polylib:

:

discrim（a*x^2+b*x+c,x）'）

.2符号数字及表达式的操作

.2.1数值数字与符号数字之间的转换

10一数值数字向符号数字的转换

在符号运算中，“数值类数字”会自动地转换为符号数字。

亦可借助sym函数：

sym（Num,'r'）或sym（Num）数值类数字Num的广义有理表达

sym（Num,'d'）数值类数字Num的“十进制浮点”近似表达

sym（Num,'e'）数值类数字Num的带eps误差的理性近似表达

sym（Num,'f'）'f'standsfor"floating-point."AllvaluesarerepresentedintheformN*2^eor-N*2^e

a=pi;

a1=sym（a,'r'）

a2=sym（a,'d'）

a3=sym（a,'e'）

a4=sym（a,'f'）

log2（281474976710656）

a=pi+0.1;

a1=sym（a,'r'）

a2=sym（a,'d'）

a3=sym（a,'e'）

a4=sym（a,'f'）

1、Num只能是数值类数字

2、sym（'Num'）和sym（Num）区别问题

'Num'数字字符串，理论值，Num表示数字，近似（双精度值）

3、在符号运算中，“数值数字”自动转换为符号数字运算sym（Num，'r'）

10二符号数字向双精度数字转换

double（Num_sym）把符号数字Num_sym转换为双精度数字

区别double（'Num'）

a=double（sym（5455454））

b=double（'34'）

.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））

e64=vpa（abs（p0-pd）,64）

digits

注意e32r等的类型，操作

e32r

.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）

[r,how]=simple（f）

r为返回的简化形式，how为化简过程中使用的一种方法。

how有以下几种形式：

（1）simplify函数对表达式进行化简；

（2）radsimp函数对含根式的表达式进行化简；

（3）combine函数将表达式中以求和、乘积、幂运算等形式出现的项进行合并；

（4）collet合并同类项

（5）factor函数实现因式分解

（6）convert函数完成表达式形式的转换

170的阶乘

factorial（170）

factorial（171）

sym（'171!

'）

.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

（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（f6）

.3符号微积分

.3.1极限和导数的符号计算

大学本科高等数学中的大多数微积分问题，都能用符号计算解决，手工笔算演绎的烦劳都可以由计算机完成。

limit（f,v,a）求极限

limit（f,v,a,'right'）求极限

limit（f,v,a,'left'）求极限

【例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>

（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

【例2.3-6】求f（x）=xex在x=0处8阶Maclaurin级数

symsx

r=taylor（x*exp（x）,9,x,0）%忽略9阶以及以上小量展开

pretty（r）%

R=evalin（symengine,'series（x*exp（x）,x=0,8）'）%<4>

pretty（R）

.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

.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）

s1=simple（simple（s））

【例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）

symsarthetaphipositive

x=r*cos（theta）;x=subs（x,r,a*theta）;

y=r*sin（theta）;y=subs（y,r,a*theta）;

dLdth=sqrt（diff（x,theta）^2+diff（y,theta）^2）;

L=simple（int（dLdth,theta,0,phi））

（2）

L_2pi=subs（L,[a,phi],sym（'[1,2*pi]'））

L_2pi_vpa=vpa（L_2pi）

（3）

L1=subs（L,a,sym（'1'））;

ezplot（L1*cos（phi）,L1*sin（phi）,[0,2*pi]）

gridon

holdon

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（'螺线长度-幅角曲线','阿基米德螺线'）

holdoff

.4微分方程的符号解法

.4.1符号解法和数值解法的互补作用

优点：

从数值计算角度，与初值问题求解相比，微分方程边值问题的求解显得复杂和困难。

对于运用数学工具去解决实际问题的科研人员来说，此时，不妨通过符号计算指令进行求解尝试。

对于符号计算来说，不论是初值问题，还是边值问题，其求解微分方程的指令形式相同，且相当简单。

缺点：

符号计算可能花费较长的时间，可能得不到简单的解析解，也可能得不到封闭形式的解，甚至可能无法求解.

求解微分方程的符号法和数值法有很好