09 符号计算

Report
Matlab 符号计算
一、创建符号对象
二、符号对象的基本操作
三、符号对象的基本运算
四、符号运算在数学中的应用
五、符号绘图
Matlab符号计算功能简介
• 符号计算以推理解析的方式进行,不受计算误差积累问题困
扰:
3
1 4
1 4
3

10
•
•
•
•

2
?
5

2

5
?
10
符号计算给出解析解,当解析解不存在时,会给出数值解;
符号计算指令的调用比较简单,与教科书上的公式相近;
计算所需时间比数值计算要多很多;
Matlab 的符号计算功能通过Symbolic Math Toolbox(符号数
学工具箱)实现。
• Matlab 2008b以前的版本,符号计算引擎为 Maple
• Matlab 2008b以后的版本(包括Maltab 2008b),符号计算引
擎为MuPAD
• 常用的符号计算软件还有Maple、Mathmatica、MathCAD等。
一、创建符号对象
1、用sym( )创建符号变量
2、用syms创建符号变量
3、符号变量的假设(assumptions)
4、创建符号表达式
5、创建符号函数
6、创建符号数组
1、 sym( )
每次只能创建 一个 符号量。
var = sym(‘var’)
x = sym('x')
y = sym('y')
z = sym('alpha')
创建符号变量var
x=
x
y=
y
z=
alpha
var = sym(Num) 由数(组)Num创建符号常量
以最接近的“有理”表
示的形式存储符号常数
a = sym(1/3 + sqrt(2))
a=
3935125774157969/2251799813685248
p , q为 正 整 数 时 ,
p
q
,
p
 ,
p ,2
q
, 10
q
称为“ 有理” 表示。
q
b = sym('1/3+sqrt(2)')
b=
2^(1/2) + 1/3
用单引号括起来,则以最精确
的形式存储符号常数
2、 syms 可以同时建立 多个 符号变量
syms var1 ... varN 创建符号变量var1 … varN
等效
syms a b c
a=sym('a');
b=sym('b');
c=sym('c');
a b c中间一定要用空格隔开,不能用逗号
syms不能用来创建符号常量
3、符号变量的假设(Assumptions)
符号变量默认都(假)设定为复(数)变量
(1) assumptions( ) 返回符号变量的假设(assumptions)
assumptions(var) 返回变量var的所有 assumptions
assumptions 返回matlab workspace中所有变量的assumptions
syms z
assumptions(z)
ans =
Empty sym: 1-by-0
当符号变量为复变量时,assumptions( )函数返回一个
空符号对象(Empty sym)
(2) sym( ) 可以假设变量为实数或正数
var = sym(‘var’,set)创建变量var并设定假设(set可取real或positive)
sym(‘var’,‘clear’)
删除变量var的假设
A = sym(A,set)
为已存在的符号数组A的元素设定假设
sym(A,‘clear’)
删除已存在的符号数组A的元素的假设
clear all
a = sym('a', 'real');
b = sym('b', 'real');
c = sym('c', 'positive');
A = sym('A%d', [1,3]);
A = sym(A, 'positive');
assumptions
sym(b, 'clear');
sym(A, 'clear');
assumptions
ans =
[ 0 < A1, 0 < A2, 0 < A3,
a in R_, b in R_, 0 < c]
ans =
[ a in R_, 0 < c]
(3) syms 也可以假设变量为实数或正数
syms var1 ... varN set
创建变量并设定假设
syms var1 ... varN clear 创建变量并删除以前的假设
clear all
syms a b real;
syms c positive;
assumptions
syms a clear;
assumptions
ans =
[ a in R_, b in R_, 0 < c]
ans =
[ b in R_, 0 < c]
(4) assume( ) 设定假设(assumption)
assume(condition):设定假设,以前的假设删除
assume(expr,set):设定set,删除以前的假定
set :‘integer’(整数)、‘rational’(有理数) 、 ‘real’(实数)
clear all
syms z;
assume(z ~= 0);
assumptions(z)
assume(z, 'integer');
assumptions(z)
ans =
z ~= 0
ans =
z in Z_
syms n
simplify(sin(2*n*pi))
ans =
sin(2*pi*n)
assume(n,'integer')
simplify(sin(2*n*pi))
ans =
0
(5) assumeAlso : 添加假设(assumption)
assumeAlso(condition):添加假设,以前的假设保留
assumeAlso(expr,set):设定set,保留以前的假定
set :‘integer’(整数)、‘rational’(有理数) 、 ‘real’(实数)
clear all
syms y
assume(y >= 0)
assumeAlso(y,'integer')
assumptions(y)
ans =
[ y in Z_, 0 <= y]
syms x
solve(x^4 == 1, x)
syms x real
solve(x^4 == 1, x)
ans =
1
-1
i
-i
ans =
1
-1
assumeAlso(x > 0);
solve(x^4 == 1, x)
ans =
1
(6) 删除符号变量及其假设
符号变量、假设是分开存储的。
syms x
assume(x,'real')
使用上面的语句创建符号变量x及假设,变量x存储在
matlab的workspace中,x的假设则存储在symbolic engine
(符号计算引擎)中
clear x
上面的语句仅删除了符号变量x,x的假设仍然存储在
symbolic engine(符号计算引擎)中。
若再创建符号变量x,则x会继承以前的假设。
syms x real
clear x
% 删除变量x
syms x
% 重新定义x,会继承之前的假设
solve(x^2 + 1 == 0, x)
Warning: Explicit solution could not be found.
> In solve at 81
ans =
[ empty sym ]
正确删除符号变量及其假设的方法:
syms x clear
clear x
% 先删除变量x的假设
% 然后再删除变量x
或:
clear all
% 删除所有变量,重启MuPAD引擎(这会
% 删除所有的假设)。
4、创建符号表达式
(1) 用 sym( ) 创建符号表达式:
f = sym('a*x^2 + b*x + c')
f=
a*x^2 + b*x + c
注意:sym()函数只创建了符号量f,并没有创建符号量a,b,c,x
(2)推荐先定义符号变量,再按普通书写形式创建符号表达式:
syms a b c x
f = a*x^2 + b*x + c;
5、创建符号函数
(1) 用sym( )创建符号函数
f(arg1,...,argN) = sym('f(arg1,...,argN)')
sym()只创建符号函数f(),但不会定义符号变量arg1,,…,argN,所
以必须先定义符号变量arg1,…,argN。
syms x y
f(x, y) = sym('f(x, y)');
(2) 用syms定义符号函数
syms f(x, y)
创建符号函数f,同时创建了符号变量x,y
f(x, y) = x + 2*y
指定函数f(x,y)的具体形式
f(1, 2)
求函数值
(3) 直接使用赋值符号 = 创建符号函数
syms x y
f(x, y) = x^3*y^3
(4) 用symfun函数创建符号函数:
f = symfun(formula,inputs)
syms x y
f = symfun(x + y, [x y])
(5) argnames 函数返回符号函数的参数(按定义时的顺序输出):
argnames(f)
(6) formula函数返回符号函数的定义式
formula(f)
例:创建符号函数数组
syms x
f(x) = [x, x^2; x^3, x^4];
f(2)
ans =
[ 2, 4]
[ 8, 16]
y = f([1 2; 3 4])
y=
[2x2 sym] [2x2 sym]
[2x2 sym] [2x2 sym]
y{1}
ans =
[ 1, 2]
[ 3, 4]
6、创建符号数组
(1) 用sym( )创建符号数组
A = sym(‘A’,dim)
A = sym('A', [3 4])
A(2, 3)
ans =
A2_3
A=
[ A1_1, A1_2, A1_3, A1_4]
[ A2_1, A2_2, A2_3, A2_4]
[ A3_1, A3_2, A3_3, A3_4]
注意:sym()只创建符号变量A,并不创
建独立的符号变量A1_1,A1_2,…
B = sym('x%d%d', [4 4])
可以修改数组元素的标记形式
B=
[ x11, x12, x13, x14]
[ x21, x22, x23, x24]
[ x31, x32, x33, x34]
[ x41, x42, x43, x44]
(2) 将数值数组转化成符号数组:
B = [2/3,sqrt(2);5.2,log(3)];
C = sym(B)
(3) 使用 sym 函数直接生成:
A=sym('[1+x,sin(x); 5,exp(x)]')
f=sym('[1,ab;c,d]')
(4) 先定义符号变量,再按普通书写 形式创建符号数组:
syms a b x y
C = [a,b;x,y]
二、符号对象的基本操作
1、symvar() 找出符号表达式(函数)中的符号变量
2、children() 返回符号表达式中的子表达式
3、pretty() 美化输出符号表达式
4、符号量与数值量之间的转换
5、符号计算精度控制
6、funtool 可视化单变量函数计算器
1、symvar( ) 找出符号表达式(函数)中的符号变量
按字母顺序列出 s 中的所有符号变量
列出 s 中离 x 最近的 n 个符号变量
symvar(s)
symvar(s,n)
若两个符号变量与 x 的距离相等,则ASCII 码大者优先。
若s为符号函数,则symvar(s,n)会先输出函数的输入变量,
然后再输出其它变量
f = sym('2*w-3*y+z^2+5*a') f =
z^2 + 5*a + 2*w - 3*y
symvar(f)
ans =
[ a, w, y, z]
查找默认变量:
symvar(f,1)
ans =
y
2、children() 返回符号表达式中的子表达式
syms x y
children(x^2 + x*y + y^2)
children(x^2 + x*y == y^2 + 1)
符号方程两边各看作是一个子表达式
ans =
[ x*y, x^2, y^2]
ans =
[ x^2 + y*x, y^2 + 1]
syms x y
s = children([x + y, sin(x)*cos(y); x^3 - y^3,
exp(x*y^2)])
s=
返回符号数组每个元素的子表达式,返回值为元胞数组
[1x2 sym] [1x2 sym]
[1x2 sym] [1x1 sym]
查看元胞数组的
具体内容:
s{1:4}
3、pretty() 美化输出符号表达式
A=sym('[3/2 - 5^(1/2)/2; 5^(1/2)/2 + 3/2]')
A=
3/2 - 5^(1/2)/2
5^(1/2)/2 + 3/2
pretty(A)
/
|
|
|
|
|
\
3
sqrt(5) \
- - ------- |
2
2
|
|
sqrt(5)
3 |
------- + - |
2
2 /
4、符号量与数值量之间的转换
(1) double( )把符号常数转化为16位相对精度的浮点数值对象
a = sym('1/123') + sym(sqrt(2)) + sym('pi')
b = double(a)
a = 1/123+2^(1/2)+pi
b = 4.5639
(2) s = poly2sym(p) 把多项式转化相应的符号表达式
p = [1,2,3]
s = poly2sym(p,x)
s=
x^2+2*x+3
(3) p = sym2poly(s) 把多项式的符号表达式转换为数值多项式
syms x
s = x + 2*x^3 - 6
p = sym2poly(s1)
p=
2
0
1
-6
5、符号计算的精度控制
digits(n) 设置符号数值计算以n位相对精度进行,默认为32位
xs = vpa(x) 在digits指定的精度下,给出x的数值型符号结果xs
xs = vpa(x,n) 在n位相对精度下,给出x的数值型符号结果xs
a = sym(pi)
digits(20)
b = vpa(a*a)
digits(32)
c = vpa(a*a)
d = vpa(a*a,40)
a = pi
b = 9.8696044010893586191
c = 9.8696044010893586188344909998761
d = 9.869604401089358618834490999876151135313
6、funtool 可视化单变量函数计算器
三、符号对象的基本运算
1、matlab符号计算中,数学运算符、关系运算符、逻辑运算
符在名称和使用上,与数值计算相同。
2、matlab符号计算中,基本的数学函数、线性代数函数、数
组操作函数,在名称和使用上,与数值计算相同。
X=sym('[a,b;c,d]')
Y=sym('[a,1;b,0]')
Z1=X*Y
Z2=X.'.*Y
Z2 =
Z1 =
[ a^2+b^2,
a]
[ a^2, c]
[ c*a+d*b,
c]
[ b^2, 0]
数学运算举例 1:
syms x
A = [x sin(x) 3];
A + x
ans =
syms x
M = [x x^2;Inf 0];
M + eye(2)
ans =
syms f(x) g(x)
f(x) = x^2 + 5*x + 6;
g(x) = 3*x - 2;
h = f + g
syms
f(x)
expr
f(x)
f(x)
= x^2 + 3*x + 2;
= x^2 - 2;
= f(x) + expr
[ 2*x, x + sin(x), x + 3]
[ x + 1, x^2]
[ Inf, 1]
h(x) =
x^2 + 8*x + 4
f(x) =
2*x^2 + 3*x
数学运算举例 2:坐标系平面转动
t = sym('t')
G = [cos(t),sin(t); -sin(t),cos(t)]
G2 = G*G
G2 =
[ cos(t)^2 - sin(t)^2, 2*cos(t)*sin(t)]
[ -2*cos(t)*sin(t), cos(t)^2 - sin(t)^2]
G2 =
[ cos(2*t), sin(2*t)]
[ -sin(2*t), cos(2*t)]
G2 = simplify(G2)
IG = inv(G)
G=
[ cos(t), sin(t)]
[ -sin(t), cos(t)]
IG =
[ cos(t)/(cos(t)^2 + sin(t)^2), -sin(t)/(cos(t)^2 + sin(t)^2)]
[ sin(t)/(cos(t)^2 + sin(t)^2), cos(t)/(cos(t)^2 + sin(t)^2)]
IG = simplify(IG)
IG =
[ cos(t), -sin(t)]
[ sin(t), cos(t)]
E = simplify(IG*G)
E=
[ 1, 0]
[ 0, 1]
矩阵运算举例:
syms a b c d
A = [a b; c d];
A*A/A
A*A-A^2
ans =
[ a, b]
[ c, d]
syms b1 b2
A = sym('a%d%d', [2 2]);
B = [b1 b2];
X = B/A;
x1 = X(1)
x2 = X(2)
x1 =
-(a21*b2 - a22*b1)/(a11*a22 - a12*a21)
x2 =
(a11*b2 - a12*b1)/(a11*a22 - a12*a21)
ans =
[ 0, 0]
[ 0, 0]
练习:用det()函数计算下列符号行列式:
a
1
0
0
x
y
x y
y
x y
x
1
b
1
0
x y
x
y
0
1
c
1
0
0
1
d
关系运算举例 1 : ==
A == B
创建符号方程
syms x
solve(sin(x) == cos(x), x)
ans =
pi/4
当 == 两边均为符号常量时,则比较两边是否相等,
返回逻辑真(1)或假(0)
A = sym(hilb(3));
B = sym([1,1/2,5; 1/2,2,1/4; 1/3,1/8,1/5]);
A == B
ans =
1 1
1 0
1 0
0
1
1
关系运算举例 2 :Isequaln()
isequaln(A,B) 若A、B的维数及内容完全相同,则返回逻辑真 (1),否则返回
逻辑假 (0)。 NaN (not a number) 值被认为彼此相等
isequaln(A1,A2,...,An) 若A1、A2、……、An完全相同,则返回逻辑真(1)
syms x
isequaln(abs(x), x)
ans =
0
用 > 设置条件
assume(x > 0)
isequaln(abs(x), x)
syms x
A1 = [x NaN NaN];
A2 = [x NaN NaN];
A3 = [x NaN NaN];
isequaln(A1, A2, A3)
ans =
1
ans =
1
关系运算举例 3 : <=
clear all
syms x
cond = (abs(sin(x)) <= 1/2);
for i = 0:sym(pi/12):sym(pi)
if subs(cond, x, i)
disp(i)
end
end
0
pi/12
pi/6
(5*pi)/6
(11*pi)/12
pi
用<=设置条件
四、符号运算在数学中的应用
(一)公式推导和化简
(二)求解符号方程(组)
(三)微积分
(一)公式推导和化简
化简
1、simplify
2、simplifyFraction
公式重排和重写
3、coeffs
4、collect
5、combine
6、compose
7、divisors
8、expand
9、factor
10、horner
11、numden
12、rewrite
替代
13、subexpr
14、subs
1. simplify()
符号对象化简
simplify(S): 对 S 进行代数化简
simplify(S,Name,Value) 用指定的选项(Name,Nalue)对S进行化简
syms x a b c
simplify(sin(x)^2 + cos(x)^2)
simplify(exp(c*log(sqrt(a+b))))
ans =
1
ans =
(a + b)^(c/2)
syms x
simplify([(x^2 + 5*x + 6)/(x + 2),...
sin(x)*sin(2*x) + cos(x)*cos(2*x);
(exp(-x*i)*i)/2 - (exp(x*i)*i)/2, sqrt(16)])
ans =
[ x + 3, cos(x)]
[ sin(x), 4]
Name-Value Pair Arguments
' Criterion ' — Simplification criterion 指定化简结果是否倾向于实
数(复数在表达式内部的几率减少)
'default' (default) | 'preferReal'
'IgnoreAnalyticConstraints' — Simplification rules 指定是否进行纯
粹的解析化简(结果可能更简单,但可能与原表达式不相等)
false(default) | true
'Seconds' — Time limit for the simplification 指定化简进程的时间
限制
processInf (default) | positive number
'Steps' — Number of simplification steps 指定化简的步骤
1 (default) | positive number
例:指定 IgnoreAnalyticConstraints 选项获得更简单的化简结果
syms x
s = (log(x^2 + 2*x + 1) - log(x + 1))*sqrt(x^2);
simplify(s)
ans =
-(log(x + 1) - log((x + 1)^2))*(x^2)^(1/2)
simplify(s, 'IgnoreAnalyticConstraints', true)
ans =
x*log(x + 1)
例:指定 Steps 选项获得更简单的化简结果
syms x
f = ((exp(-x*i)*i)/2 - (exp(x*i)*i)/2)/(exp(x*i)/2 + exp(x*i)/2);
simplify(f)
ans =
-(exp(x*2*i)*i - i)/(exp(x*2*i) + 1)
simplify(f, 'Steps', 10)
simplify(f, 'Steps', 30)
simplify(f, 'Steps', 50)
ans =
(2*i)/(exp(x*2*i) + 1) - i
ans =
((cos(x) - sin(x)*i)*i)/cos(x) - i
ans =
tan(x)
2. simplifyFraction():分式化简
simplifyFraction(expr) :化简结果的分子、分母都为多项式,且最大公因式为1
simplifyFraction(expr,Name,Value) :用指定的选项(Name,Nalue)进行化简
Name-Value Pair Arguments
'Expand'
Expand the numerator and denominator of the resulting fraction
(将化简结果的分子、分母展开)
Default: false
syms x y
simplifyFraction((x^2 - 1)/(x + 1))
simplifyFraction(((y+1)^3*x)/((x^3-x*(x+1)*(x-1))*y))
ans =
x-1
ans =
(y + 1)^3/y
syms x y
simplifyFraction(((y + 1)^3*x)/((x^3 - x*(x + 1)*(x 1))*y),'Expand', true)
ans =
(y^3 + 3*y^2 + 3*y + 1)/y
3. coeffs()
提取多项式系数
C = coeffs(p) :返回单变量多项式p的系数组成的向量C,变量为默认变量
C = coeffs(p,var):返回单变量多项式p的系数组成的向量C ,变量为var
C = coeffs(p,vars):返回多变量多项式p的系数组成的向量C ,多变量为vars
[C,T] = coeffs(___):返回多项式p的系数组成的向量C及对应的项T
syms x
c = coeffs(16*x^2 + 19*x + 11)
c=
[ 11, 19, 16]
注意:这里c是升幂排列的
syms x
[c,t] = coeffs(16*x^2 + 19*x + 11)
c=
[ 16, 19, 11]
注意:这里c是降幂排列的
t=
[ x^2, x, 1]
4. collect() Collect coefficients(合并同类项)
collect(P)
按 默认变量 降幂顺序重写符号表达式P
collect(P,var) 按 变量 var降幂顺序重写符号表达式P
syms x
collect((exp(x) + x)*(x + 2))
ans =
x^2 + (exp(x) + 2)*x + 2*exp(x)
syms x y
collect(x^2*y + y*x - x^2 - 2*x, x)
collect(x^2*y + y*x - x^2 - 2*x, y)
ans =
(y - 1)*x^2 + (y - 2)*x
ans =
(x^2 + x)*y - x^2 - 2*x
5. combine()
Y = combine(S)
Y = combine(S,T)
合并指定的代数结构
将S中出现的多个幂的乘积(底数相同)合并为一个
将S中多次出现的函数T合并为一个
syms x y z
combine(x^y*x^z)
ans =
x^(y + z)
syms a b
assume(abs(a*b) < 1)
combine(atan(a) + atan(b),'atan')
ans =
-atan((a + b)/(a*b - 1))
求复合函数
6、compose()
com pose  f , g  :返 回 f
 g  y  , 其 中 f
 f
x, g
 g  y 。
这 里 的 x , y 分 别 是 由 sym var 函 数 确 定 的 f , g 中 的 默 认 变 量 。
com pose  f , g , z  :返 回 f
 g  z  , 其 中 f
 f
x, g
 g  y 。
这 里 的 x , y 分 别 是 由 sym var 函 数 确 定 的 f , g 中 的 默 认 变 量 。
 g  z  , 指 定 x为 f 的 变 量 。
com pose  f , g , x , y , z  :返 回 f  g  z   , 指 定 x 为 f 的 变 量 , y 为 g 的 变 量 。
clear all
e=
syms x y z t u
exp(-z/u)^t
h = x^t;
p = exp(-y/u);
f=
e = compose(h,p,x,y,z)
x^exp(-y/z)
f = compose(h,p,t,u,z)
com pose  f , g , x , z  :返 回 f
可用symfun函数将符号表达式转为符号函数:
f(x, y, z) =
f = symfun(f,symvar(f))
x^exp(-y/z)
7、Divisors Divisors of integer or expression
(R2014b中的新函数)
divisors(n) :返回整数n的所有非负数因子
divisors(expr,vars) :返回多项式expr的因式,vars为变量
divisors(sym(42))
ans = [ 1, 2, 3, 6, 7, 14, 21, 42]
syms x
divisors(x^4 - 1, x)
ans =
[ 1, x - 1, x + 1, (x - 1)*(x + 1), x^2 + 1, (x^2 + 1)*(x - 1),...
(x^2 + 1)*(x + 1), (x^2 + 1)*(x - 1)*(x + 1)]
8. expand()
展开符号表达式
expand(S)把S展开为多项式形式,也可展开三角函数、指数函数、对数函数
expand(S,Name,Value) 指定选项(Name,Nalue)进行展开
Name-Value Pair Arguments
'ArithmeticOnly‘
Default: false
若为 true,只展开算术部分,不展开三角函数、双曲
函数、对数函数及特殊函数。
'IgnoreAnalyticConstraints‘
Default: false
若为true,进行纯代数化简,结果可能更简单,但可
能与原表达式不相等
syms x
expand((x-2)*(x-4))
syms x y
expand(cos(x+y))
ans =
x^2 - 6*x + 8
ans =
cos(x)*cos(y) - sin(x)*sin(y)
expand(log((a*b/c)^2), 'IgnoreAnalyticConstraints', true)
ans =
2*log(a) + 2*log(b) - 2*log(c)
9. factor() 因式分解
factor(f)
对 f 进行因式分解,也可用于正整数的分解
ans =
syms x y
factor(x^3 - y^3)
(x-y)*(x^2+x*y+y^2)
a = 1234567890
b1 = factor(a)
b2 = factor(sym(a))
输入数值量,返回一个行数组
输入符号量,返回符号表达式
b1 =
2
3
3
b2 =
(2)*(3)^2*(5)*(3803)*(3607)
5
3607
3803
10、horner()
输出horner多项式
horner(P)
将多项式 P 改写为horner形式(数值计算中,这
种形式可以快速计算多项式的值)
syms x
horner(x^3 - 6*x^2 + 11*x - 6)
ans =
x*(x*(x - 6) + 11) - 6
syms x y
horner([x^2 + x; y^3 - 2*y])
ans =
x*(x + 1)
y*(y^2 - 2)
11. numden()
Numerator(分子)and denominator(分母)
[N,D] = numden(A) 将符号或数值数组A转为有理形式N/D,N、D
为互质多形式(系数为整数)
syms x y
[n,d] = numden(x/y + y/x)
n=
x^2 + y^2
d=
x*y
syms a b
A = [a, 1/b]
[n,d] = numden(A)
n=
[ a, 1]
d=
[ 1, b]
12. rewrite() 用指定函数重写符号表达式
rewrite(expr,target 用target重写符号表达式expr
rewrite(A,target)
用target重写符号数组A中的每个元素
target
可以是: exp, log, sincos, sin, cos, tan, sqrt, heaviside.
syms x
rewrite(sin(x), 'exp')
ans =
(exp(-x*i)*i)/2 - (exp(x*i)*i)/2
syms x
rewrite(tanh(x), 'sin')
ans =
(sin(x*i)*i)/(2*sin((x*i)/2)^2 - 1)
syms x
rewrite(acos(x), 'log')
ans =
-log(x + (1 - x^2)^(1/2)*i)*I
13、subexpr()
重写字符串
[r,sigma] = subexpr(expr) 重写符号表达式expr为r,expr中的共同的
子表达式用符号sigma代替(即使输出变量不是sigma,r里出现的
也是sigma)。
[r,var] = subexpr(expr,var) 共同的子表达式用符号变量var代替(必
须先定义符号变量var)。
[r,var] = subexpr(expr,‘var’) 与[r,var] = subexpr(expr,var)功能相同。
syms a b c d x
solutions = solve(a*x^3 + b*x^2 + c*x + d == 0, x)
求解结果 solutions 超复杂,很长!!!!可用pretty()函数美化一下输出进行查看。
有三个解,每个解里都有一个很长的共同的字表达式。
[r, sigma] = subexpr(solutions)
sigma =
((d/(2*a) + b^3/(27*a^3) - (b*c)/(6*a^2))^2 + (- b^2/(9*a^2) +
c/(3*a))^3)^(1/2) - b^3/(27*a^3) - d/(2*a) + (b*c)/(6*a^2)
14. subs()
替换符号量
R = subs(S,old,new)
用 new 替换 S 中的 old
R = subs(S, new)
用 new 替换 S 中的缺省变量
R = subs(S)
用workspace中的量替换 S 中的相应变量
注意:替换的结果被赋值给R,S本身并没有改变
clear all
syms a x b
f = a*x+b
c = subs(f,a,5)
c=
b + 5*x
d = subs(f,5)
d=
5*a + b
b = 2
e = subs(f)
e=
a*x + 2
若要同时替换多个变量,写成向量形式即可
syms a b
subs(cos(a) + sin(b), [a, b], [sym('alpha'), 2])
ans =
sin(2) + cos(alpha)
也可以用元胞数组:
subs(cos(a) + sin(b), {a, b}, {sym('alpha'), 2})
ans =
sin(2) + cos(alpha)
当用数组替换一个符号标量时,会将数组中的每个元素进行替
换,返回一个同维数组
syms a t
subs(exp(a*t) + 1, a, -magic(3))
ans =
[ exp(-8*t) + 1, exp(-t) + 1, exp(-6*t) + 1]
[ exp(-3*t) + 1, exp(-5*t) + 1, exp(-7*t) + 1]
[ exp(-4*t) + 1, exp(-9*t) + 1, exp(-2*t) + 1]
syms x y
subs(x*y, {x, y}, {[0 1; -1 0],
[1 -1; -2 1]})
ans =
[ 0, -1]
[ 2, 0]
将数组A的第一个元素替换为数组B,输出数组会扩充为4×4数组
A = sym('A', [2,2])
B = sym('B', [2,2])
A44 = subs(A, A(1,1), B)
A=
[ A1_1, A1_2]
[ A2_1, A2_2]
B=
[ B1_1, B1_2]
[ B2_1, B2_2]
A44 =
[ B1_1, B1_2, A1_2, A1_2]
[ B2_1, B2_2, A1_2, A1_2]
[ A2_1, A2_1, A2_2, A2_2]
[ A2_1, A2_1, A2_2, A2_2]
(二)求解符号方程(组)
1. equationsToMatrix
2. linsolve
3. solve
4. vpasolve
5. poles
6. finverse
7. dsolve
将线性方程组转为矩阵形式
解矩阵形式的线性方程组
解符号方程(组)
求符号方程(组)数值解
求表达式或函数的奇点
求反函数
解微分方程(组)
1. equationsToMatrix 将线性方程组转为矩阵形式
[A,b] = equationsToMatrix(eqns,vars)
[A,b] = equationsToMatrix(eqns)
A = equationsToMatrix(eqns,vars)
A = equationsToMatrix(eqns)
syms x y z
[A, b] = equationsToMatrix([x + y - 2*z == 0, x + y + z
== 1, 2*y - z + 5 == 0], [x, y, z])
A=
[ 1, 1, -2]
[ 1, 1, 1]
[ 0, 2, -1]
b=
0
1
-5
X=
3
-7/3
1/3
X = A\b
syms s t
[A,b]=equationsToMatrix([s - 2*t + 1 == 0, 3*s - t == 10])
A=
[ 1, -2]
[ 3, -1]
b=
-1
10
X = A\b
X=
21/5
13/5
2. linsolve 解矩阵形式的线性方程组
X = linsolve(A,B)
求解 AX = B
[X,R] = linsolve(A,B)
求解 AX = B,若A为方阵,返回A
的条件数(1范数)的倒数R,否则返回A的秩R
syms
eqn1
eqn2
eqn3
x
=
=
=
y z
2*x+y+z == 2;
-x+y-z == 3;
x+2*y+3*z == -10;
2x + y + z = 2
−x + y − z = 3
x + 2y + 3z = −10
[A,B] = equationsToMatrix([eqn1, eqn2, eqn3], [x, y, z])
A=
[ 2, 1, 1]
[ -1, 1, -1]
[ 1, 2, 3]
B=
2
3
-10
X = linsolve(A,B)
X=
3
1
-5
也可以直接创建符号数组A和b,求线性方程组的符号解:
x  y  z  a

3 x  y  6 z  7
 y  3z  4

clc;
clear all;
A = sym('[1,1,1;3,-1,6;0,1,3]')
b = sym('[a;7;4]')
s = A\b
s=
-14/15 + 3/5*a
-3/5 + 3/5*a
23/15 - 1/5*a
3. solve
解方程(组)
注意:建议用syms定义符号变量,然后创建方程和方程组。未
来版本的solve将不支持使用字符串来表示方程和方程组。
S = solve(eqn,var) 解方程eqn,变量为var
Y = solve(eqns,vars)解方程组eqns,解Y为结构数组
[y1,...,yN] = solve(eqns,vars) 解方程组,解被赋值为独立变量
syms x
solve(x^4 + 1 == 2*x^2 - 1)
ans =
(1 + i)^(1/2)
(1 - i)^(1/2)
-(1 + i)^(1/2)
-(1 - i)^(1/2)
syms x y
S = solve(x + y == 1, x - 11*y == 5)
S=
x: [1x1 sym]
y: [1x1 sym]
S = [S.x S.y]
S=
[ 4/3, -1/3]
syms a b
[b, a] = solve(a + b == 1, 2*a - b == 4, b, a)
b=
-2/3
a=
5/3
solve 在得不到解析解时,会给出数值解
syms x
solve(sin(x) == x^2 - 1)
ans =
-0.63673265080528201088799090383828
例:求范德瓦尔斯气体的三个临界参量
a 

P

V  b   R T

2 
V 

 2P 
 P 
,
0
  0, 
2 
 V T
 V T
Pc , V c , T c  ?
clc;
clear all;
syms a b R T V
P(V,T) = R*T/(V-b) - a/V^2
PdV = diff(P,V)
PdV2 = diff(PdV,V)
S = solve(PdV==0,PdV2==0,V,T)
Vc = S.V
Tc = S.T
Pc = P(Vc,Tc)
4. vpasolve
求符号方程(组)的数值解
(R2012b以后版本中的新函数)
S = vpasolve(eqn,var)
解方程eqn,变量为var
Y = vpasolve(eqns,vars)解方程组eqns,解Y为结构数组
[y1,...,yN] = vpasolve(eqns,vars) 解方程组,解被赋值为独立变量
syms x
vpasolve(sin(x^2) == 1/2, x)
ans =
-226.94447241941511682716953887638
syms x y
S = vpasolve([x^3 + 2*x == y, y^2 == x], [x, y])
S=
x: [6x1 sym]
y: [6x1 sym]
查看具体数值: S.y
S.y
5. poles
求表达式或函数的奇点
poles(f,var)
求f的奇点(变量为var)
P = poles(f,var)
求f的奇点,结果赋值给向量P
[P,N] =poles(f,var)
求f的奇点、极点的阶,结果赋值给P、N
[P,N,R]= poles(f,var) 求f的奇点、阶、留数,结果赋值给P、N、R
poles(f,var,a,b)
求f在区间(a,b)内的奇点
P = poles(f,var,a,b)
[P,N] =poles(f,var,a,b)
[P,N,R]= poles(f,var,a,b)
syms x a
[Poles,Orders,Residues] = poles(a/x^2/(x - 1),x)
Poles =
1
0
Orders =
1
2
Residues =
a
-a
6. finverse
求反函数
g = finverse(f)
g = finverse(f,var)
返回f的反函数,自变量为默认变量
返回f的反函数,自变量为var
syms x
f(x) = 1/tan(x);
g = finverse(f)
g(x) =
atan(1/x)
syms u v
finverse(exp(u - 2*v), u)
ans =
2*v + log(u)
7. dsolve
解常微分方程(组)
S = dsolve(eqn)
解常微分方程eqn
S = dsolve(eqn,cond) cond为初值(边界)条件
Y = dsolve(eqns)
解方程组eqns,结果赋值给结构Y
Y = dsolve(eqns,conds)
[y1,...,yN] = dsolve(eqns)
解方程组eqns,结果赋值给变量y1,y2,...
[y1,...,yN] = dsolve(eqns,conds)
syms a x(t)
dsolve(diff(x) == -a*x)
ans =
C2*exp(-a*t)
syms a b y(t)
dsolve(diff(y) == a*y, y(0) == b)
ans =
b*exp(a*t)
syms f(t) g(t)
[f(t), g(t)] = dsolve(diff(f) == 3*f + 4*g,...
diff(g) == -4*f + 3*g, f(0) == 0, g(0) == 1)
f(t) =
sin(4*t)*exp(3*t)
g(t) =
cos(4*t)*exp(3*t)
(三)符号计算在微积分中的应用
•
•
•
•
•
•
•
•
1. limit
2. diff
3. int
4. symsum
5. symprod
6. taylor
7. taylortool
8. 积分变换
求极限
求微分
求积分
级数求和
级数求积
函数的泰勒级数展开
单变量泰勒级数展开图形工具
1.limit
求极限
limit(expr,x,a)
limit(expr,a)
limit(expr)
limit(expr,x,a,‘left’)
limit(expr,x,a,‘right’)
求x→a时expr的极限
求默认变量→a时expr的极限
求默认变量→0时expr的极限
求x→a-时expr的左极限
求x→a+时expr的右极限
syms x h
limit(sin(x)/x)
limit((sin(x + h) - sin(x))/h, h, 0)
ans =
1
ans =
cos(x)
syms x
limit(1/x, x, 0, 'right')
limit(1/x, x, 0, 'left')
ans =
Inf
ans =
-Inf
2. diff
求导数
diff(F)
求F的一阶导数,变量为默认变量
diff(F,var)
求F的一阶导数,变量为var.
diff(F,n)
求F的n阶导数,变量为默认变量
diff(F,var,n) 求F的n阶导数,变量为var,也可写成diff(F,n,var)
diff(F,var1,...,varN) 求F对变量var1,...,varN的混合偏导数
syms x
f(x) = sin(x^2);
df = diff(f)
syms x y
diff(x*cos(x*y), y, 2)
syms x y
diff(x*sin(x*y), x, y)
df(x) =
2*x*cos(x^2)
ans =
-x^3*cos(x*y)
ans =
2*x*cos(x*y) - x^2*y*sin(x*y)
3. int
求积分
int(expr,var)
int(expr,var,a,b)
求expr的不定积分,变量为var
求expr的定积分
若不指定var,则有symvar确定默认变量。a:积分下限,b:积分上限。
当 a 或 b 取 inf(或 -inf)时,计算的就是广义积分。
syms x z
int(x/(1 + z^2), x)
int(x/(1 + z^2), z)
ans =
x^2/(2*(z^2 + 1))
ans =
x*atan(z)
syms x
int(x*log(1 + x), 0, 1)
ans =
1/4
syms x
int(1/(x^2-2*x+3),-inf,inf)
ans =
(pi*2^(1/2))/2
练习:求下面的积分,给出50位精度的数值
2
x
2
  
1
x
2
x y
( x  y  z ) dzdydx
2
xy
先对 z 变量的积分
然后对 y 变量的积分
最后对 x 变量的积分
2
2
4、级数求和
symsum(expr,var)
求exp指定项的和,变量为var (0,1,2,…,var-1)
symsum(expr,var,a,b) 求exp指定项的和,变量为var (a,a+1,a+2,…,b)
syms k
symsum(k)
ans =
k^2/2 - k/2
syms k
symsum(k^2, 0, 10)
symsum(1/k^2,1,Inf)
syms k x
symsum(x^k/sym('k!'), k, 0, Inf)
ans =
385
ans =
pi^2/6
ans =
exp(x)
级数求积
5. symprod
symprod(expr,var)求exp指定项的积,变量为var (1,2,…,var)
symprod(expr,var,a,b)求exp指定项的积,变量为var (a,a+1,…,b)
syms k
symprod(k)
symprod((2*k - 1)/k^2)
ans =
factorial(k)
ans =
(1/2^(2*k)*2^(k + 1)*factorial(2*k))/(2*factorial(k)^3)
syms k
symprod(1 - 1/k^2, k, 2, Inf)
symprod(k^2/(k^2 - 1), k, 2, Inf)
ans =
1/2
ans =
2
6、taylor
函数的泰勒级数展开
在 x  a处 展 开 式 : f  x   f  a  
在 x  0处 展 开 : f  x   f  0  
f ' a 

f ''  a 
1!
f '0 
1!



2!

f ''  0 
 x  a
n0


n

2!

n0
x
n
f
n
f
n
a
n!
0
n!
taylor(f)
求 f 在默认变量 = 0 处的 5 阶泰勒展开式
taylor(f,v)
求 f 在变量 v = 0 处的 5 阶泰勒展开式
taylor(f,v,a) 求 f 在变量 v = a 处 5 阶泰勒展开式
taylor(_, ‘Order‘,n) 求 n-1 阶泰勒展开式(不指定Order时,n默认为6)
syms x
f = sin(x)/x;
t6 = taylor(f)
t6 =
x^4/120 - x^2/6 + 1
t10 = taylor(f, 'Order', 10)
t10 =
x^8/362880 - x^6/5040 + x^4/120 - x^2/6 + 1
7. taylortool
单变量泰勒展开图形工具
8. 积分变换
(1) fourier()
(2) ifourier()
(3) laplace()
(4) ilaplace()
(5) ztrans()
(6) iztrans()
Fourier积分变换
逆Fourier积分变换
Laplace变换
逆Laplace变换
Z-变换
逆Z-变换
这些函数的具体使用,请参考matlab的帮助文件。
五、绘制符号函数的图形
1. ezplot
绘制直角坐标系二维曲线
2. ezpolar 绘制极坐标二维曲线
3. ezplot3 绘制三维参数曲线
4. ezmesh 绘制三维参数曲线
5. ezmeshc 绘制带等高线的三维网格图
6. ezsurf
绘制三维表面图
7. ezsurfc 绘制带等高线的三维表面图
8. ezcontour
绘制等高线
9. ezcontourf
绘制填充等高线
1. ezplot
绘制直角坐标系二维曲线
ezplot(f)
绘制f的图形,单(或双)变量区间为 [–2π 2π]
ezplot(f,[min,max])
单(或双)变量区间为 [min,max]
ezplot(f,[xmin,xmax,ymin,ymax]) 双变量变量区间为 [xmin,xmax,ymin,ymax]
ezplot(f,fign)
在第fign个图形窗口绘制f的图形
ezplot(x,y)
绘制参数方程x = x(t), y = y(t) 的图形, 0 <= t <= 2π
ezplot(x,y,[tmin,tmax]) 绘制 x = x(t) , y = y(t)的图形, tmin <= t <= tmax.
ezplot(f,figure_handle) 在句柄figure_handle确定的窗口中绘制f的图形
syms x y
f(x, y) = sin(x + y)*sin(x*y);
ezplot(f)
syms t
x = t*sin(5*t);
y = t*cos(5*t);
ezplot(x, y)
2. ezpolar 绘制极坐标二维曲线
ezpolar(f) plots the polarcurve r = f(θ)over the default domain 0 < θ < 2π.
ezpolar(f, [a, b]) plots f for a < θ < b.
syms t
ezpolar(1 + cos(t))
1


练习: 
2

sin
7
t

cos
30
t
 
  ,t 
8 
2

1  

100   t 

2



100
3 
 1
 2  , 2  


syms t
f = 100/(100+(t-1/2*pi)^8)*(2-sin(7*t)-1/2*cos(30*t))
ezpolar(f,[-pi/2,3*pi/2])
3. ezplot3 绘制三维参数曲线
ezplot3(x,y,z)
plots the spatialcurve x = x(t), y = y(t),and z = z(t) over the default
domain 0 < t < 2π.
ezplot3(x,y,z,[tmin,tmax])
plotsthe curve x = x(t), y = y(t),and z = z(t) over the domain
tmin < t < tmax.
ezplot3(...,'animate')
producesan animated trace of the spatial curve.
syms t
ezplot3(sin(t), cos(t), t,[0,6*pi])
绘制三维网格图
4. ezmesh
ezmesh(f) creates a graph of f(x,y), –2π < x < 2π, –2π < y < 2π.
ezmesh(f, domain) plots f overthe specified domain. domain canbe either a 4-by-1 vector
[xmin, xmax, ymin, ymax]or a 2-by-1 vector [min, max] (where, min < x < max, min < y < max).
ezmesh(x,y,z) plots the parametric surface x = x(s,t), y = y(s,t),and z = z(s,t)over the square
–2π < s < 2π, –2π < t < 2π.
ezmesh(x,y,z,[smin,smax,tmin,tmax]) or ezmesh(x,y,z,[min,max]) plotsthe
parametric surface using the specified domain.
ezmesh(...,n) plots f overthe default domain using an n-by-n grid.The default value for n is
60.
ezmesh(...,'circ') plots f overa disk centered on the domain.
f
 x, y  
2
xe
x y
2
syms x y
ezmesh(x*exp(-x^2-y^2),[-2.5,2.5],40)
colormap([0 0 1])
绘制带等高线的三维网格图
5. ezmeshc
ezmeshc(f) creates a graph of f(x,y), –2π < x < 2π, –2π < y < 2π.
ezmeshc(f,domain) plots f overthe specified domain. domain canbe either a 4-by-1
vector [xmin, xmax, ymin, ymax]or a 2-by-1 vector [min, max] (where, min < x < max, min < y
< max).
ezmeshc(x,y,z) plots the parametric surface x = x(s,t), y = y(s,t),and z = z(s,t)over the
square –2π < s < 2π, –2π < t < 2π.
ezmeshc(x,y,z,[smin,smax,tmin,tmax]) or ezmeshc(x,y,z,[min,max])
plotsthe parametric surface using the specified domain.
ezmeshc(...,n) plots f overthe default domain using an n-by-n grid.The default value for n
is 60.
ezmeshc(...,'circ') plots f overa disk centered on the domain.
f
 x, y  
y
1 x  y
2
2
syms x y
ezmeshc(y/(1 + x^2 + y^2),[-5,5,-2*pi,2*pi])
6. ezsurf
绘制三维表面图
ezsurf(f) plots over the default domain –2π < x < 2π, –2π < y < 2π.
ezsurf(f,domain) plots f overthe specified domain. domain canbe either a 4-by-1 vector
[xmin, xmax, ymin, ymax]or a 2-by-1 vector [min, max] (where, min < x < max, min < y <
max).
ezsurf(x,y,z) plots the parametricsurface x = x(s,t), y = y(s,t),and z = z(s,t)over the square
–2π < s < 2π, –2π < t < 2π.
ezsurf(x,y,z,[smin,smax,tmin,tmax]) or ezsurf(x,y,z,[min,max]) plots the
parametricsurface using the specified domain.
ezsurf(...,n) plots f overthe default domain using an n-by-n grid.The default value for n is
60.
ezsurf(...,'circ') plots f overa disk centered on the domain.
syms x y
ezsurf(real(atan(x+i*y)))
7. ezsurfc 绘制带等高线的三维表面图
ezsurfc(f) creates a graphof f(x,y), –2π < x < 2π, –2π < y < 2π
ezsurfc(f,domain) plots f overthe specified domain. domain canbe either a 4-by-1
vector [xmin, xmax, ymin, ymax] ora 2-by-1 vector [min, max] (where, min < x < max, min < y
< max).
ezsurfc(x,y,z) plots the parametricsurface x = x(s,t), y = y(s,t),and z = z(s,t)over the square
–2π < s < 2π, –2π < t < 2π.
ezsurfc(x,y,z,[smin,smax,tmin,tmax]) or ezsurfc(x,y,z,[min,max]) plotsthe
parametric surface using the specified domain.
ezsurfc(...,n) plots f overthe default domain using an n-by-n grid.The default value for n is
60.
ezsurfc(...,'circ') plots f overa disk centered on the domain.
f
 x, y  
y
1 x  y
2
2
syms x y
ezsurfc(y/(1 + x^2 + y^2),[-5,5,-2*pi,2*pi],35)
8. ezcontour
绘制等高线
ezcontour(f) plots the contourlines of f(x,y),–2π < x < 2π, –2π < y < 2π.
ezcontour(f,domain) plots f(x,y) over the specified domain. domain canbe either a 4by-1 vector [xmin, xmax, ymin, ymax] or a 2-by-1 vector [min, max] (where, min < x < max,
min < y < max).
ezcontour(...,n) plots f overthe default domain using an n-by-n grid.The default value for
n is 60.
ezcontour automatically adds a title andaxis labels.
2
2
2
2
1   x 1 2  y 2
2
x
 x   y 1
3
5  x y
f  x , y   3 1  x  e
 10   x  y  e
 e
3
5

syms x y
f = 3*(1-x)^2*exp(-(x^2)-(y+1)^2)...
- 10*(x/5 - x^3 - y^5)*exp(-x^2-y^2)...
- 1/3*exp(-(x+1)^2 - y^2);
ezcontour(f,[-3,3],49)
9. ezcontourf
绘制填充等高线
ezcontourf(f) plots the contourlines of f(x,y), –2π < x < 2π, –2π < y < 2π.
ezcontourf(f,domain) plots f(x,y)over the specified domain. domain canbe either a 4by-1 vector [xmin, xmax, ymin, ymax]or a 2-by-1 vector [min, max] (where, min < x < max, min
< y < max).
ezcontourf(...,n) plots f overthe default domain using an n-by-n grid.The default value for
n is 60.
ezcontourf automatically adds a title andaxis labels.
syms x y
f = 3*(1-x)^2*exp(-(x^2)-(y+1)^2)...
- 10*(x/5 - x^3 - y^5)*exp(-x^2-y^2)...
- 1/3*exp(-(x+1)^2 - y^2);
ezcontourf(f,[-3,3],49)

similar documents