第三讲MATLAB的符号运算

Download Report

Transcript 第三讲MATLAB的符号运算

第三讲 MATLAB的符号运算
• 科学与工程技术中的数值运算固然重要,但自
然科学理论分析中各种各样的公式、关系式及
其推导就是符号运算要解决的问题。
• 在Matlab7.0中,符号计算虽以数值运算的补充
身份出现,但它们都是科学计算研究的重要内
容。
• Matlab开发了实现符号计算的工具包Symbolic
Math Toolbox 。
• 符号数学工具箱中的工具是建立在功能强
大的Maple的基础上。
• 它最初是由加拿大的滑铁卢(Waterloo)大
学开发出来的。
• 如果要求Matlab7.0进行符号运算,那么首
先由Maple计算并将结果返回到Matlab7.0
命令窗口。
两个数学分析的可视化界面
• 图示化符号计算器
(由命令funtool引出)
• 泰勒级数逼近分析界面
(由命令taylortool引出)
图示化符号计算器
• 由三个独立的窗口构成,通过函数运算
控制窗口来演示另外两个图形窗口,任
何时候,只有一个窗口属于激活状态。
而被激活的函数图像可随运算控制窗口
的操作而做相应的变化。
• 下面给出运算控制窗口的键位功能。
• 前两行是函数 f 和 g 的具体解析式,第三
•
•
•
•
行是自变量 x 的取值范围和常数 a 的值。
第四行只对 f 起作用,如求导、积分、简
化、提取分子和分母、倒数、反函数。
第五行是处理 f 和 a 的加减乘除等运算。
第六行前四个进行 f 和 g 之间的运算,后
三个分别是:求复合函数;把 f 传递给 ;
swap是实现 f 和 g 功能的交换。
最后一行是对计算器自身进行操作。
• Funtool计算器存有一张函数列表fxlist
•
•
•
•
•
•
•
这7个功能键分别是:
Insert:把当前激活窗的函数写入列表
Cycle:依次循环显示fxlist中的函数
Delete:从fxlist列表中删除激活窗的函数
Reset:使计算器恢复到初始调用状态
Help:获得关于界面的在线提示说明
Demo:自动演示
Close:关闭整个计算器
泰勒级数逼近分析
• 该界面用于观察函数f(x)在给定区间被N
阶泰勒多项式Tn(x)逼近的情况。
• f(x)的输入可由命令taylortool(fx)引入,
或者在栏中直接输入表达式,回车确定。
• N默认值为7,a是级数的展开点。
• 函数的观察区间默认为(-2pi,2pi)。
符号运算的功能
•
•
•
•
•
•
符号表达式、符号矩阵的创建
符号线性代数
因式分解、展开和简化
符号代数方程求解
符号微积分
符号微分方程
一、符号运算的基本操作
1. 什么是符号运算
• 与数值运算的区别
※ 数值运算中必须先对变量赋值,
然后才能参与运算。
※ 符号运算无须事先对独立变量
赋值,运算结果以标准的符号形式
表达。
• 特点:
 运算对象可以是没赋值的符号变量,以推理解析的方
式进行,因此不受计算误差累积所带来的困扰。
 可以给出完全正确的封闭解或任意精度的数值解(当封
闭解不存在时)。
③符号计算指令的调用简单,和经典教科书公式相近。
④计算所需的时间较长。
• Symbolic Math Toolbox——符号运算工具包通过调用
Maple软件实现符号计算的。
• Maple软件——主要功能是符号运算,它占据符号软件
的主导地位。
2. 字符串与符号变量、符号常量
字符串对象
f = 'sin(x)+5x'
f —— 字符串名
sin(x)+5x—— 函数表达式
'
'—— 字符串标识
字符串表达式一定要用' '单引
号括起来Matlab才能识别。
用class( )来返回对象的数据类型。
‘
’ 里的内容可以是函数表达式,也
可以是方程。
例:
f1='a*x^2+b*x+c' —— 二次三项式
f2= 'a*x^2+b*x+c=0' —— 方程
f3='Dy+y^2=1' ——微分方程
※函数表达式或方程可以赋给字符串
或符号变量,以后方便调用。
符号变量
• 符号变量是内容可变的符号对象。
• 符号变量通常是指一个或几个特定的字
符,不是指符号表达式,甚至可以将一
个符号表达式赋值给一个符号变量。
• 符号变量有时也称自由变量,它的命名
规则和数值变量的命名规则相同。
• 相关指令为: sym( ) 和 syms( )
(symbolic的缩写)
例:用函数命令sym( )和syms( )来创建符号对
象并检测数据类型。
a=sym('a')
注意两个 a的区别
b=sym('c')
classa=class(a)
classb=class(b) 可看出两个变量均为符号对象
syms a b c d e f g h
whos
也可以查看所有变量类型
从上述比较来看:当需要同时定义多个符号
变量时,使用syms( )更简洁一些。
符号常量
• 当数值常量作为sym( )的输入参量时,就
建立了一个符号对象——符号常量。
• 虽然看上去是一个数值量,但已经是一
个符号对象了。
例:a=3/4; b='3/4'; c=sym(3/4); d=sym('3/4');
whos 查看变量类型
a为实双精度浮点数值类型;b为实字符类
型;c和d都是符号对象类型。
由符号变量构成的符号函数和
符号方程
• 符号表达式是由符号常量、符号变量、符号函
数运算符以及专用函数连接起来的符号对象。
• 包括:符号函数和符号方程。判断看带不带等
号。
例:syms x y z; f1=x*y/z;
f2=x^2+y^2+z^2; f3=f1/f2;
e1=sym('a*x^2+b*x+c')
e2=sym('sin(x)^2+2*cos(x)=1')
e3=sym('Dy-y=x')
3.符号矩阵的创建
数值矩阵 clear clc A=[1,2;3,4]
A=[a,b;c,d] —— 不识别
用Matlab函数sym创建矩阵
命令格式:A=sym('[
]')
※ 符号矩阵内容同数值矩阵
※ 需用sym指令定义,需用'
※ 注意与'[a,b;c,d]'的区别
'标识
例如:A = sym('[a , 2*b ; 3*a , 0]')
A=
[
a, 2*b]
[3*a,
0]
这就完成了一个符号矩阵的创建。
注意:符号矩阵的每一行的两端都有方
括号,这是与 Matlab数值矩阵的
一个重要区别。
用字符串直接创建矩阵
 模仿Matlab数值矩阵的创建方法
 需保证同一列中各元素字符串有相
同的长度。
例:A =['[ a,2*b]'; '[3*a, 0]']
A=
[ a, 2*b]
[3*a, 0]
 符号矩阵的修改
a.直接修改
可用、 键找到所要修改的矩阵,
直接修改
b.指令修改
 用A1=subs(A, 'new', 'old')来修改
例如:A =[ a, 2*b]
[3*a, 0]
A(2,2)='4*b'
A1 = [ a, 2*b]
[3*a, 4*b]
A2=subs(A1, 'c', 'b')
A2 =[ a, 2*c]
[3*a, 4*c]
 符号矩阵与数值矩阵的转换
将数值矩阵转化为符号矩阵
函数调用格式:sym(A)
clear A=[1/3,2.5;1/0.7,2/5]
A=
0.3333 2.5000
1.4286 0.4000
sym(A)
ans =
[ 1/3, 5/2]
[10/7, 2/5]
将符号矩阵转化为数值矩阵
函数调用格式: numeric(A)???
A=
[ 1/3, 5/2]
[10/7, 2/5]
numeric(A)
ans =
0.3333 2.5000
1.4286 0.4000
二、符号运算
由于Matlab7.0采用了重载技术,使得符号计
算表达式的运算符和基本函数,无论在
形状、名称上,还是在使用方法上,都
与数值计算中的运算符和基本函数几乎
完全相同。这无疑给用户带来了极大的
方便。
例外:在符号对象的比较中,没有”大
于”、 ”大于等于”、 ”小于”、 ”
小于等于”的概念,而只有是否“等于”
1.符号矩阵运算
数值运算中,所有矩阵运算操作指令都比
较直观、简单。例如:a=b+c; a=a*b ;
A=2*a^2+3*a-5等。
符号运算中,很多方面在形式上同数值计
算都是相同的,没必要重新学习新的规则。
2. 任意精度的数学运算
在symbolic中有三种不同的算术运算:
1. 数值类型
matlab的浮点算术运算
2. 有理数类型 maple的精确符号运算
3. vpa类型
maple的任意精度算术
运算
• 浮点算术运算
format long --(定义输出格式) 1/2+1/3
ans =
0.83333333333333
• 符号运算
sym(1/2)+(1/3)
或sym(1/2+1/3)
ans =
5/6 --精确解
• 任意精度算术运算
digits(n) —— 设置近似解的精读为n位有效数字,
默认32位有效数字。
vpa(x,n) —— 求符号解的近似解,该近似解的有
效位数由n来决定。
digits(25)
vpa(1/2+1/3)
ans =
.8333333333333333333333333
vpa(5/6,40)
ans =
.8333333333333333333333333333333333333333
a=sym('[1/4,exp(1);log(3),3/7]')
a=
[ 1/4,exp(1)]
[log(3), 3/7]
vpa(a,10)
ans =
[.2500000000, 2.718281828]
[1.098612289, .4285714286]
3.符号表达式的化简
• 可以对符号计算结果进行简化,诸如因
式分解、同类项合并、符号表达式的展
开、符号表达式的化简和通分等等。
• 合并同类项 collect(v) ----将表达式v的相
同次幂的项合并。
例:syms x t
% 定义基本变量
f=(x-1)*(x-2)*(x-3) %定义符号表达式
collect(f)
%合并f中x的同类项
• expand(s)
将s中的各项进行展开,用于
多项式,三角函数、指数函数、对数函数。
例:syms x y;
f=(x+y)^3;
f1=expand(f)
f1 =
x^3+3*x^2*y+3*x*y^2+y^3
例:h=cos(x-y)
expand(h)
• factor(S)
将系数为有理数的多项式(矩
阵)S,表示成低阶多项式相乘的形式,
如果不能分解,则返回S本身。
例:syms x y
factor(x^3-y^3)
• simplify( ) 该函数是一个强有力的具有
普遍意义的工具,它利用Maple化简规则
对表达式进行简化。
例:S=sym('[(x^2+5*x+6)/(x+2);sqrt(16)]')
simplify(S)
• simple( ) 用几种不同的算术简化规则对
符号表达式进行简化,使其用最少的字
符来表示。
• 虽然并非表达式中的字符越少,表达式
就越简单,但采用这个标准往往能够得
到满意的结果,尤其是对于包含三角函
数的表达式。
例:sym x
simple(cos(x)^2+sin(x)^2)
• 从结果看出,simple比较这些不同函数的
结果,最终把最少字符作为标准。
4. 符号微积分与积分变换
• diff(f) — 对缺省变量求f的微分
• diff(f,v) — 对指定变量v求微分
• diff(f,n) — 对默认变量求n阶微分
• diff(f,v,n) —对指定变量v求f的n阶微分
例:syms a x
f=sin(a*x)
df=diff(f)
dfa=diff(f,a,2)
符号表达式的极限
•
•
•
•
•
limit(F,x,a) 求当x→a时,表达式F的极限
limit(F, a)
默认自变量时,趋于a的极限
limit(F)
默认自变量,默认a=0
limit(F,x,a, 'left')
取F的左极限
limit(F,x,a, 'right')
取F的右极限
例:syms h n x
dc=limit((sin(x+h)-sin(x))/h,h,0)
%按照导数的定义求sin的导数
注意:对于极限不存在,返回NaN
例: limit(1/x,x,0)
limit(1/x,x,0, 'left')
limit(1/x,x,0, 'right')
结果分别为:
• ans =
NaN
• ans =
-Inf
• ans =
Inf
符号表达式的积分
• int(f) — 对f表达式的缺省变量求不定积分
• int(f,v) — 对f表达式的v变量求不定积分
• int(f,v,a,b) — 对f表达式的v变量在(a,b)
区间求定积分
findsym(f) —可以找出f中的每个变量
注意:当函数的积分不存在时,Matlab7.0将简单地
返回原来的积分表达式。
int('被积表达式','积分变量','积分上限',
'积分下限')—— 定积分
——缺省时为不定积分
例:int(-2*x/(1+x^2)^2)
ans =
1/(1+x^2)
int(log(x))
int(log10(x))
int(sin(x),x,-pi,pi)
• taylor(f,n,v) —— n阶泰勒级数展开
例:syms x
f=1/(2+cos(x))
r=taylor(f,8)
• symsum(f,v,a,b) —表达式f中变量 v从a变
到b时的有限和
例:syms x k
s1=symsum(1/k^2,1,inf)
s2=symsum(x^k,k,0,inf)
上述都是求无穷级数的和
符号积分变换
•
•
•
•
•
•
ztrans(f) —— Z变换
Invztrans(f) —— 反Z变换
Laplace(f) —— 拉氏变换
Invlaplace(f) —— 反拉氏变换
fourier(f) —— 付氏变换
Invfourier(f) —— 反付氏变换
注意 :上述函数均缺省了部分参数
符号积分的例子
例1.计算二重不定积分
xe
dxdy

 xy
F=int(int('x*exp(-x*y)','x'),'y')
例2.计算 f='x*exp(-x*10)'的Z变换
F=ztrans(f)
F=
z*exp(-10)/(z-exp(-10))^2
>> syms x y
>> F=int(int(x*exp(-x*y),x),y)
F=
1/y*exp(-x*y)
>> syms x
>> f=x*exp(-x*10);
>> F=ztrans(f)
>> F=ztrans(x*exp(-x*10);
F=
z*exp(-10)/(z-exp(-10))^2
例3. 计算指数函数eAt。
用拉氏反变换法计算eAt的公式为:
eAt = L-1[(SI-A)-1]
0
1
系统矩阵A=
2 3
结果:
2t
t
At
e =
2t
t
 e  2e  e  e
 2t
t
 2t
t
2e  2e 2e  e
>> a=[0 1;-2 -3];
>> syms s
>> b=(s*eye(2)-a)
b=
[ s, -1]
[ 2, s+3]
>> B=inv(b)
[ (s+3)/(s^2+3*s+2),
[ -2/(s^2+3*s+2),
1/(s^2+3*s+2)]
s/(s^2+3*s+2)]
>> b11=ilaplace(sym(b,1,1));b(1,1)=b11;
>> b12=ilaplace(sym(b,1,2));b(1,2)=b12;
>> b21=ilaplace(sym(b,2,1));b(2,1)=b21;
>> b22=ilaplace(sym(b,2,2));b(2,2)=b22;
>> b
b=
[ -exp(-2*t)+2*exp(-t),
exp(-t)-exp(-2*t)]
[ -2*exp(-t)+2*exp(-2*t), 2*exp(-2*t)-exp(-t)]
5.符号代数方程求解
Matlab符号运算能够解一般的线性
方程、非线性方程、超越方程。当方
程组不存在符号解时,又无其他自由
参数,则给出数值解。
命令格式:
solve(f,v) —— 求一个方程f=0的解
Solve(f1,f2, …fn) —— 求n个方程的解
例1. f = ax2+bx+c 求解
f='a*x^2+b*x+c';
• solve(f) —— 对缺省变量x求解
ans =
计算机
[1/2/a*(-b+(b^2-4*a*c)^(1/2))]
格式
[1/2/a*(-b-(b^2-4*a*c)^(1/2))]
 b  b  4ac
2a
2
一般格式
• solve(f , 'b' ) —— 对指定变量b求解
' '加与不加效果一样
例2. 符号方程cos(x)=sin(x)
tan(2*x)=sin(x)求解
f1=solve('cos(x)=sin(x)'),
f1 =
1/4*pi
f2=solve('tan(2*x)=sin(x)')
例3. 解方程组
x+y+z=1
x-y+z=2
2x-y-z=1
g1='x+y+z=1',g2='x-y+z=2',g3='2*x-y-z=1'
f=solve(g1,g2,g3)
f=solve('x+y+z=1','x-y+z=2','2*x-y-z=1')
f=
z = 5/6, y = -1/2, x = 2/3
f=solve('x+y+z=1','x-y+z=2','2*x-y-z=1')
f=
x: [1x1 sym] f.x
ans =2/3
y: [1x1 sym] f.y
ans =-1/2
z: [1x1 sym] f.z
ans =5/6
[x,y,z]=solve('x+y+z=1','x-y+z=2','2*x-y-z=1')
x = 2/3
y =-1/2
z =5/6
6. 符号微分方程求解
—— 用一个函数可以方便地得到微
分方程的符号解
符号微分方程求解指令:dsolve
命令格式:dsolve(f,g)
• f —— 微分方程,可多至12个微分方程的求
解;g为初始条件
• 默认自变量为 'x',可任意指定自变量't', 'u'等
• 微分方程的各阶导数项以大写字母D表示
dy
dy
或
y的一阶导数—— Dy
dt
dx
d2y
d2y
2 或
2 y的二阶导数—— D2y
dt
dx
n
dny
d y
y的 n 阶导数—— Dny
n 或
n
dt
dx
[y1,y2…]=dsolve(x1,x2,…xn) —— 返回
微分方程的解
一阶微分方程
dsolve('Dx=y','Dy=x','x(0)=0','y(0)=1')
ans =
x(t) = sin(t), y(t) = cos(t)
二阶微分方程
dsolve('D2y=-a^2*y','y(0)=1','Dy(pi/a)=0')
ans =
cos(a*x)
2
d y
dy
2
 2y  0
例.
2
dx
dx
dy
y (0)  1, (0) 0 求该方程的解
dx
y=dsolve('D2y+2*Dy+2*y=0','y(0)=1','Dy(0)=0')
ans =
exp(-x)*cos(x)+exp(-x)*sin(x)
ezplot(y) —— 方程解y(t)的时间曲线图
exp(-x)*cos(x)+exp(-x)*sin(x)
300
250
200
150
100
50
0
-6
-4
-2
0
x
2
4
exp(-t) sin(t)+exp(-t) cos(t)
250
200
150
100
50
0
-6
-4
-2
0
t
2
4
三、maple函数——符号运算的扩展
maple——是专门进行数学运算的软件工具,
具有超强的符号运算能力,提供了
几乎包括所有数学领域的专用函数
matlab——依赖于maple的内核与函数库,扩
展了自己的符号运算功能。
matlab还设计了对maple库函数的调用功能
使得已有的maple数学功能,可以扩充matlab
中,作为自身符号运算能力的扩展。
1. maple内核访问函数
可以访问maple内核的matlab函数:
maple ——— 访问maple内核函数
mapleinit —— maple函数初始化
mpa ———— maple函数定义
mhelp ——— maple函数帮助命令
procread —— maple函数程序安装
. maple 的调用格式
maple('表达式')
—— 将表达式送至maple内核,
返回符号表达式结果。
maple ('函数',变量1,变量2)
——调用maple函数,传递给定
变量。
例1. 展开5阶 bernoulli 多项式,计
算 x=3 时bernoulli 数。
a=maple('bernoulli(5,x)')
a=
-1/6*x+5/3*x^3+x^5-5/2*x^4
a=maple('bernoulli(5,3)')
a=
85
例2. 化简三角函数式sin2x+cos2x
a=maple('simplify(sin(x)^2+cos(x)^2);')
a=
1
例4. 求f(t)=e-3tsint的拉式变换
f=maple('laplace(exp(-3*t)*sin(t),t,s);')
f=
1/((s+3)^2+1)
例4. 寻找二次多项式的完全平方
f (x) = x2+2x+2
a=maple('completesquare(x^2+2*x+2)')
a=
completesquare(x^2+2*x+2) 将工具包
装入内存
maple('with(student);')
a=maple('completesquare(x^2+2*x+2)')
a=
(x+1)^2+1
• maple软件中的所有函数,在初
始化时并没有完全装入内存,可
用readlib指令把库函数读入内存,
或用with指令将应用工具包装入
内存。
• 调用格式
maple('readlib(函数名);')
maple('with(工具包名);')
例5.求sin(x2+y2)在x=0,y=0处泰勒级数
展开式,8阶截断。
maple('mtaylor(sin(x^2+y^2),[x=0,y=0],8)')
ans =
mtaylor(sin(x^2+y^2),[x = 0, y = 0],8)
maple('readlib(mtaylor);')
maple('mtaylor(sin(x^2+y^2),[x=0,y=0],8)')
ans =
x^2+y^2-1/6*x^6-1/2*y^2*x^4-1/2*y^4*x^21/6*y^6
2. mpa —— maple变量定义
• 任何一个matlab定义的函数f,可使用mpa
语句直接调用,还可把 f 定义成maple变
量v。
• maple的工作空间与matlab工作空间是相
互独立的, 所以f 与v是属于不同工作空
间中的变量
• mpa的调用格式:
mpa('v',f)
f为matlab工作空间中已存在的变量
mpa v f
S 0 ( NI )
例. 电磁力计算公式为 F ( I , x) 
2
4x
试I=0.5,x=0.1邻域展开泰勒级数,3阶截
断,令常数 k  S 0 N , F ( I , x)  kI
2
2
4
2
x
1.直接调用
maple('readlib(mtaylor);')
maple('mtaylor(k*I^2/x^2,[I=0.5,x=0.1],3);')
2
2.定义符号函数f(matlab6.1无map函数)
f='k*I^2/x^2';
maple('mtaylor(f,[I=0.5,x=0.1],3);')
ans =
mtaylor(f,[I = .5, x = .1],3)
mpa('u',f)
maple('mtaylor(u,[I=0.5,x=0.1],3);')
ans =
25.*k-.50e3*k*(x-.1)+.10e3*k*(I.5)+7500.000000000000*k*(x.1)^2+.1e3*k*(I-.5)^2-.20e4*k*(I-.5)*(x-.1)
注意:matlab符号运算时,可以识
别matlab定义的符号变量,但在调
用 maple 函数时,需将matlab变量
定义为maple变量后,所调用的函
数方可识别和执行
3.mhelp —— maple函数帮助命令
 mhelp 是协助检索maple库函数的专
用命令
工具词条
调用格式:mhelp 相关词条
函数词条
例如:
mhelp intro — maple介绍
mhelp maple — maple命令格式
mhelp tutorial —maple入门
mhelp index —maple检索
mhelp index 用于工具包检索
library ——maple标准库函数
packages —— 应用工具包
libmisc —— 其它库函数
statements —— maple语句描述
expressions —— maple表达式
datatypes —— maple数据格式
tables —— maple表格和阵列
procedures —— maple程序
misc —— maple其它应用
一般帮助文本主要包括以下部分
• FUNCTION—— 函数功能说明
• CALLING SEQUENCE—— 调用格式
• PARAMETERS —— 调用参数说明
• SYNOPSIS —— 语法说明
• EXAMPLES —— 应用举例
• SEE ALSO —— 相关词条
4.maple库函数
maple库函数共分四类
maple内部函数:驻留函数任何条件
下都可调用
mhelp index[internal]
 maple的外部函数—读库定义部分:
调用时先执行读库命令,因此与内部函
数一样可直接调用
mhelp index[external]
 maple的外部函数—读库装入部分
maple其余外部函数需要在使用前执
行maple('readlib(函数名);')命令将其
装入内存
mhelp index[libmisc]
 maple的惰性函数—不能直接调用,
还需一些函数如mod,evala,evalf 等才
能调用
mhelp index[intert]
小
结
本节介绍了matlab语言的符号运算
功能,通过学习应该掌握:
• 掌握如何创建、修改符号矩阵
• 掌握符号运算功能
• maple函数调用
• mhelp检索