<%@LANGUAGE="VBSCRIPT" CODEPAGE="936"%>西南科技大学理学院-数学在线首页
   
 
数学软件

 

 

教学辅导

第8讲 符号计算


符号数学工具箱是操作和解决符号表达式的符号数学工具箱(函数)集合,有复合、简化、微分、积分以及求解代数方程和微分方程的工具。另外还有一些用于线性代数的工具,求解逆、行列式、正则型式的精确结果,找出符号矩阵的特征值而无由数值计算引入的误差。工具箱还支持可变精度运算,即支持符号计算并能以指定的精度返回结果。
符号数学工具箱中的工具是建立在功能强大的称作Maple软件的基础上。它最初是由加拿大的滑铁卢(Waterloo)大学开发的。当要求MATLAB进行符号运算时,它就请求Maple去计算并将结果返回到MATLAB命令窗口。因此,在MATLAB中的符号运算是MATLAB处理数字的自然扩展。
8.1 符号表达式
符号表达式是代表数字、函数、算子和变量的MATLAB字符串,或字符串数组。不要求变量有预先确定的值,符号方程式是含有等号的符号表达式。符号算术是使用已知的规则和给定符号恒等式求解这些符号方程的实践,它与代数和微积分所学到的求解方法完全一样。符号矩阵是数组,其元素是符号表达式。
MATLAB在内部把符号表达式表示成字符串,以与数字变量或运算相区别;否则,这些符号表达式几乎完全象基本的MATLAB命令。
下表列有几则符号表达式例子以及MATLAB等效表达式。
符号表达式 MATLAB表达式
'1/(2*x^n)'
y='1/sqrt(2*x)'
'cos(x^2)-sin(2*x)'
M=sym('[a,b;c,d]')
f=int('x^3/sqrt(1-x)','a','b')
MATLAB符号函数使我们能用多种方法来操作符号表达式,比如,
>>diff('cos(x)') %differentiate cos(x) with respect to x
ans=-sin(x)
>>M=sym('[a,b;c,d]') %create a symbolic matrix M
M=
[a,b]
[c,d]
>>determ(M) %find the determinant of the symbolic matrix M
ans=a*d-b*c
要注意的是,以上第一例的符号表达式是用单引号以隐含方式定义的。它告诉MATLAB 'cos(x)'是一个字符串并说明diff('cosx')是一个符号表达式而不是数字表达式;然而在第二个例子中,用函数sym显式地告诉MATLAB M=sym('[a,b;c,d]')是一符号表达式。在MATLAB可以自己确定变量类型的场合下,通常不要求显式函数sym。
MATLAB中函数function argument形式是与function( 'argument ')等价的。其中,function是一个函数,argument是一字符串。
例如,diff cos(x)和diff( 'cos(x) ')两者都意味diff (sym 'cos(x) ')。但第一种形式显然更便于输入。然而,很多时候sym是必要的。在上述的第二个例子中,
>>M=[a,b;c,d] % M is a numeric matrix using value of a through d
???Uundefine function or variable a.
>>M= '[a,b;c,d] ' % M is a character string, but not a symbolic matrix
M=[a,b;c,d]
>>M=sym('[a,b;c,d]') %M is a symbolic matrix
M= [a,b]
[c,d]
M以三种方式定义: 数字型(如果a、b、c、d已预先确定)、字符串型或符号矩阵型。许多符号函数非常巧妙能够自动将字符转变为符号表达式。但在某些情况下,尤其是建立符号数组时,必须用函数sym,特别地将字符串变为符号表达式。隐含形式,例如diff cos(x),对于那些不需要参考先前结果的简单任务,最有用。但是最简单形式(无引号)要求一个参量,它是一个单字符的字符串、不包含插入的空格。
>>diff x^2+3*x+5 %the argument is equivalent to 'x^2+3*x+5 '
ans=2*x+3
>> diff x^2 + 3*x + 5 % spaces break the argument into separate strings
???Error using==>diff
Too manyinput arguments
无变量的符号表达式称作符号常量。符号常量常常与整数很难区别,例如
>>f=symop('(3*4-2)/5+1') %reduce a symbolic constant to its simplest form
f=3
>>isstr(f) %is f a string? (1=yes,0=no)
ans=1
在这个例子中,f代表符号常数'3';而不是数字3。MATLAB是以字符ASCII码形式来存储字符串的。所以,如果对字符串进行数字运算,则在运算中,采用各字符串的ASCII码值。因为数字51是字符 ' 3 ' 的ASCII表示,所以f加1在数值上不能得到期望的结果
>>f+1
ans=52
当字符表达式中含有多于一个的变量时,只有一个变量是独立变量。如果不告诉MATLAB哪一个变量是独立变量,MATLAB将基于以下规则选择一个:
在符号表达式中缺省的独立变量是唯一的,除去i和j的小写字母,不是单词的一部分。如果没有这种字母,就选择x作为独立变量。如字符不是唯一的,就选择在字母顺序中最接近x的字母。如果有相连的字母,就选择在字母表中较后的那一个。
缺省的独立变量,有时称作自由变量,在表达式'1/(5+cos(x))'中是' x';在'3*y+z'中是'y';在'a+sin(t)'是't'。在表式' sin(pi/4)-cos(3/5) '中自由符号变量是'x ',因为此式是一个符号常数无符号变量。可利用函数symvar询问MATLAB在符号表达式中哪一个变量它认为是独立变量。
>>symvar('a*x+y*) %find the default symbolic variable
ans=x
>>symvar('a*t+s/(u+3)') %u is the closest to 'x'
ans=u
>>symvar('sin(omega)') % 'omega' is not a singlee character。
ans=x
>>symvar('3*i+4*j') %i and j are equel to sqrt(-1)
ans=x
>>symvar('y+3*s','t') % find the variable closest to t rather than x
ans=s
如果利用规则symvar不能找到一个缺省独立变量,它便假定无独立变量并返回x。这一结论对含有由多个字母组成的变量,如:alpha或s2的表达式,或不含变量的符号常数均成立。如果需要,绝大多数命令都使用用户选项以指定独立变量。
>>diff('x^n') %differentiate with respect to the default variable 'x'
ans=x^n*n/x
>> diff('x^n','n') % differentiate x^n with respect to 'n'
ans=x^n*log(x)
>>diff('sin(omega)') %differentiate using the default variables (x)
ans=0
>>diff( 'sin(omega)','omega') %specify the independent variable
ans=cos(omega)
8.2 符号表达式运算
一旦创建了一个符号表达式,或许想以某些方式改变它;也许希望提取表达式的一部分,合并两个表达式或求得表达的数值。有许多符号工具可以帮助完成这些任务。
提取分子和分母:
如果表达式是一个有理分式(两个多项式之比),或是可以展开为有理分式(包括哪些分母为1的分式),可利用numden来提取分子或分母。例如,给定如下的表达式:

在必要时,numden将表达式合并、有理化并返回所得的分子和分母。进行这项运算的MATLAB语句是:
>>m='x^2' %create a simple expression
m=x^2
>>[n,d]=numden(m) %extract the numerator and denominator
n=x^2
d=1
>>f='a*x^2/(b-x)' %create a rational expression
f=a*x^2/(b-x)
>>[n,d]=numden(f) %extract the numerator and denominator
n=a*x^2
d=b-x
前二个表达式得到期望结果。
>>g='3/2*x^2+2/3*x-3/5' %rationalize and extract the parts
g=3/2*x^2+2/3*x-3/5
>>[n,d]=numden(g)
n=45*x^2+20*x-18
d=30
>>h='(x^2+3)/(2*x-1)+3*x/(x-1)' %the sum of rational polynomials
h=(x^2+3)/(2*x-1)+3*x/(x-1)
>>[n,d]=numden(h) %rationalize and extract
n=x^3+5*x^2-3
d=(2*x-1)*(x-1)
在提取各部分之前,这二个表达式g和h被有理化,并变换成具有分子和分母的一个简单表达式。
>>k=sym('[3/2,(2*x+1)/3;4/x^2,3*x+4]') % try a symbolic array
k=[ 3/2,(2*x+1)/3]
[4/x^2, 3*x+4]
>>[n,d]=numden(k)
n=[3, 2*x+1]
[4, 3*x+4]
d=[ 2,3]
[x^2,1]
这个表达式k是符号数组,numden返回两个新数组n和d,其中n是分子数组,d是分母数组。如果采用s=numden(f)形式,numden仅把分子返回到变量s中。
8.3 标准代数运算
很多标准的代数运算可以在符号表达式上执行,函数symadd、symsub、symlnul和symdiv为加、减、乘、除两个表达式,sympow将一个表达式上升为另一个表达式的幂次。例如: 给定两个函数

>>f='2*x^2+3*x-5' %define the symbolic expression
f=2*x^2+3*x-5
>>g='x^2-x+7'
g=x^2-x+7
>>symadd(f,g) % find an expression for f+g
ans=3*x^2+2*x+2
>>symsub(f,g) % find an expression for f-g
ans=x^2+4*x-12
>>symmul(f,g) %find an expression for f*g
ans=(2*x^2+3*x-5)*(x^2-x+7)
>>symdiv(f,g) %find an expression for f/g
ans=(2*x^2+3*x-5)/(x^2-x+7)
>>sympow(f, ' 3*x ' ) % find an expression for
ans=(2*x^2+3*x-5)^3**
另一个通用函数可让用户用其它的符号变量、表达式和算子创建新的表达式。symop取由逗号隔开的、多至16个参量。各个参量可为符号表达式、数值或算子(' + '、' - '、'*'、' / '、' ^ '、' ( '或' ) '),然后symop可将参量联接起来,返回最后所得的表达式.
>>f='cos(x)' % create an expression
f=cos(x)
>>g='sin(2*x)' % create another expression
g=sin(2*x)
>>symop(f,'/ ',g,'+',3) %combine them
ans=cos(x)/sin(2*x)+3
所有这些运算也同样可数组参量进行。
8.4 高级运算
MATLAB具有对符号表达式执行更高级运算的功能。函数compose把f(x)和g(x)复合成f(g(x))。函数finverse求表达式的函数逆,而函数symsum求表达式的符号和。
给定表达式
>>f='1/(1+x^2'; % create the four expression
>>g='sin(x)';
>>h='1/(1+u^2)';
>>k='sin(v)';
>>compose(f,g) % find an expression for f(g(x))
ans=1/(1+sin(x)^2)
>>compose(g,f) % find an expression for g(f(x))
ans=sin(1/(1+x^2))
compose也可用于含有不同独立变量的函数表达式。
>>compose(h,k,'u','v') % given h(u),k(v),find(k(v))
ans=1/(1+sin(v)^2)
表达式譬如f(x)的函数逆g(x),满足g(f(x))=x。例如, 的函数逆是ln(x),因为ln( )=x。sin(x)的函数逆是arcsin(x),函数 的函数逆是arcsin 。函数fincerse返回表达式的函数逆。如果解不是唯一就给出警告。
>>finverse('1/x') %the inverse of 1/x is 1/x since '1/(1/x)=x '
ans=1/x
>>finverse('x^2') % g(x^2)=x has more than one solution
Warning: finverse(x^2) is not unique
ans=x^(1/2)
>>finverse( ' a*x+b ' ) % find the solution to ' g(f(x))=x '
ans=-(b-x)/a
>>finverse('a*b+c*d-a*z'),'a') %find the solution to 'g(f(a))=a '
ans=-(c*d-a)/(b-z)
symsun函数求表达式的符号和有四种形式:symsun(f)返回 ;symsum(f,'s ')返回 ,symsun(f,a,b)返回 ;最普通的形式symsun(f,'s',a,b)返回 。
试一试 ,它应返回: 。
>>symsum('x^2')
ans=1/3*x^3-1/2* x^2+1/6*x
又怎么样呢?它应返回 。
>>sym('(2*n-1)^2',1,'n')
ans=11/3*n+8/3-4*(n+1)^2+4/3*(n+1)^3
>>factor(ans) % change the form ( we will revisit 'factor' later on)
ans=1/3*n*(2*n-1)*(2*n+1)
最后让我们试一试 ,其返回应是 。
>>symsum('1/(2*n-1)^2 ',1,inf)
ans=1/8*pi^2
8.5 变换函数
本节提出许多工具,将符号表达式变换成数值或反之。有极少数的符号函数可返回数值。然而请注意,某些符号函数能自动地将一个数字变换成它的符号表达式,如果该数字是函数许多参量中的一个。
函数sym可获取一个数字参量并将其转换为符号表达式。函数numneric的功能正好相反,它把一个符号常数(无变量符号表达式)变换为一个数值。
>>phi='(1+sqrt(5))/2' % the ' golden ' ratio
phi=(1+sqrt(5))/2 % convert to a numeric value
>>numeric(phi)
ans=1.6180
函数eval将字符串传给MATLAB以便计算。所以eval是另一个可用于把符号常数变换为数字或计算表达式的函数。
>>eval(phi) % execute the string ' (1+sqrt(5))/2
ans=1.6180
正如所期望那样,numeric和eval返回相同数值。
符号函数sym2poly将符号多项式变换成它的MATLAB等价系数向量。函数poly2syrn功能正好相反,并让用户指定用于所得结果表达式中的变量。
>>f=' 2*x^2+x^3-3*x+5 ' % f is the symbolic polynomials
f=2*x^2+x^3-3*x+5
>> n=sym2poly(f) % extract eht numeric coefficient vector
n=1 2 -3 5
>>poly2sym(n) % recreate the polynomials in x (the default)
ans=2*x^2+x^3-3*x+5
>> poly2sym(n,' s ') % recreate the polynomials in s
ans=s^3+2*s^2-3*s+5
8.6 变量替换
假设有一个以x为变量的符号表达式,并希望将变量转换为y。MATLAB提供一个工具称作subs,以便在符号表达式中进行变量替换。其格式为subs(f,new,old),其中f是符号表达式,new和old是字符、字符串或其它符号表达式。'新'字符串将代替表达式f中各个'旧'字符串。以下有几个例子:
>> f= ' a*x^2+b*x+c ' % create a function f(x)
f=a*x^2+b*x+c
>> subs(f,' s ',' x ') % substitute ' s ' for ' x ' in the expression f
ans=a*s ^2+b*s+c
>>subs(f,' alpha ',' a ') % substitute ' alpha ' for ' a ' in f
ans=alpha*x^2+b*x+c
>>g=' 3*x^2+5*x-4 ' % create another function
g=3*x^2+5*x-4
>>h=subs(g,' 2 ',' x ') % substitute ' 2 ' for ' x ' in g
h=18
>>isstr(h) % show that the result is a symbolic expression
ans=1
最后一个例子表明subs如何进行替换,并力图简化表达式。因为替换结果是一个符号常数,MATLLAB可以将其简化为一个符号值。注意,因为subs是一个符号函数,所以它返回一个符号表达式。尽管看似数字,实质上是一个符号常数。为了得到数字,我们需要使用函数numeric或eval来转换字符串。
>>numeric(h) % convert a symbolic expression to a number
ans=18
>>isstr(ans) % show that the result is a numeric value
ans=0
8.7 微分和积分
微分和积分是微积分学研究和应用的核心,并广泛地用在许多工程学科。MATLAB符号工具能帮助解决许多这类问题。
微分
符号表达式的微分以四种形式利用函数diff:
>>f= ' a*x^3+x^2-b*x-c ' % define a symbolic expression
f=a*x^3+x^2-b*x-c
>>diff(f) %differentiate with respect to the default variable x
ans=3*a*x^2+2*x-b
>>diff(f,'a ') % differentiate with respect to a
ans=x^3
>>diff(f,2) % differentiate twice with respect to x
ans=6*a*x+2
>>diff(f,' a ',2) % differentiate twice with respect to a
ans=0
函数diff也可对数组进行运算。如果F是符号向量或数组,diff(F)对数组内的各个元素进行微分。
>> F=sym('[a*x, b*x^2;c*x^3,d*s]') % create a symbolic array
F=[ a*x, b*x^2]
[c*x^3, d*s]
>>diff(F) % differentiate the element with respect to x
ans=[ a,2*b*x]
[3*c*x^2, 0]
注意函数diff也用在MATLAB计算数值向量或矩阵的数值差分。对于一个数值向量或矩阵M,diff(M)计算M(2:m,:)-M(1:m-1,:)的数值差分,如下所示:
>>m=[(1:8).^2)] % create a vector
M=1 4 9 16 25 36 49 64
>>diff(M) % find the differences between elements
ans=3 5 7 9 11 13 15
如果diff的表达式或可变参量是数值,MATLAB就非常巧妙地计算其数值差分;如果参量是符号字符串或变量,MATLAB就对其表达式进行微分。
积分
积分函数int(f),其中f是一符号表达式,它力图求出另一符号表达式F使diff(F)=f。正如从研究微分学所了解的,积分比微分复杂得多。积分或逆求导不一定是以封闭形式存在,或许存在但软件也许找不到,或者软件可明显地求解,但超过内存或时间限制。当MATLAB不能找到逆导数时,它将返回未经计算的命令。
>>int('log(x)/exp(x^2)') % attempt to integrate
ans=log(x)/exp(x^2)
同微分一样,积分函数有多种形式。形式int(f)相对于缺省的独立变量求逆导数;形式(f,'s')相对于符号变量s积分;形式int(f,a,b)和int(f,'s',a,b),a,b是数值,求解符号表达式从a到b的定积分;形式int(f,'m ' ,'n')和形式int(f,'s','m','n'),其中m,n是符号变量,求解符号表达式从m到n的定积分。
>>f='sin(s+2*x)') % crate a symbolic function
f=sin(s+2*x)
>>int(f) % integrate with respect to x
ans=-1/2*cos(s+2*x)
>>int(f,' s ') % integrate with respect to s
ans=-cos(s+2*x)
>>int(f,pi/2,pi) %integrate with respect to x from /2 to
ans=-cos(x)
>>int(f,' s ',pi/2,pi) %integrate with respect to s from /2 to
ans=cos(2*x)-sin(2*x)
>>int(f,' m ','n') % integrate with respect to x from m to n
ans=-1/2*cos(s+2*n)+1/2*cos(s+2*m)
正如函数diff一样,积分函数int对符号数组的每一个元素进行运算。
>> F=sym( '[a*x,b*x^2;c*x^3,d*s] ' ) % create a symbolic array
F=[ a*x,b*x^2]
[c*x^3, d*s]
>>diff(F) % ubtegrate the array elements with respect to x
ans=[1/2*a*x^2,1/3*b*x^3]
[1/4*c*x^4, d*s*x]

8.8 符号表达式简化和格式化
有时MATLAB返回的符号表达式难以理解,有许多工具可以使表达式变得更易读懂。第一个就是函数pretty,该命令以类似于数学课本上的形式来显示符号表达式。
>>f=taylor('log(x+1)/(x-5)') % 6 terms is the default
f=-1/5*x+3/50*x^2-41/750*x^3+293/7500*x^4-1207/37500*x^5+O(x^6)
>>pretty(f)
符号表达式可用许多等价形式来提供。在不同的场合,某种形式可能胜于另一种。MATLAB用许多命令来简化或改变符号表达式。
>>f=sym('(x^2-1)*(x-2)*(x-3)') % create a function
f=(x^2-1)*(x-2)*(x-3)
>>collect(f) % collect all like terms
ans=x^4-5*x^3+5*x^2+5*x-6
>>hornor(ans) % change to Horner or nested representation
ans=-6+(5+(5+(-5+x)*x)*x)*x
>>factor(ans) % express as a product of polynomials
ans=(x-1)*(x-2)*(x-3)*(x+1)
>>expand(f) % distribute products over sums
ans=x^4-5*x^3+5*x^2+5*x-6
Simplify是功能强大、通用的工具。它利用各种类型代数恒等式,包括求和、积分和分数幂、三角、指数和log函数、Bessel函数、超几何函数和 函数等来简化表达式。
下面例子说明函数的乘幂。
>>simplify('log(2*x/y)')
ans=log(2)+log(x)-log(y)
>>simplify('sin(x)^2+3*x+cos(x)^2-5')
ans=3*x-4
>>simplify('(-a^2+1)/(1-a)')
ans=a+1
其中要讨论的最后一个函数是最有用的,但也是最不正统的。函数simple试用了几种不同的简化工具,然后选择在结果表达式中含有最少字符的那种形式。如立方根:
>>f='(1/x^3+6/x^2+12/x+8)^(1/3)' %create the expression
f=(1/x^3+6/x^2+12/x+8)^(1/3)
>>simple(f) % simplify it
simplify: (2*x+1)/x
ans=(2*x+1)/x
>>simple(ans) % try it once again - another method may help
combine(trig)
2+1/x
ans=2+1/x
正如所见,simple试用了几种可简化表达式的简化方式,并让看到每一个尝试的结果。有时,它帮助多次使用函数simple并对第一次的结果作不同的简化操作,如上所作。simple对于含有三角函数的表达式尤为有用。让我们试一下
>> simple('cos(x)+sqrt(-sin(x)^2)') % simplify a trig expression
simplify:cos(x)+(cos(x)^2-1)^(1/2)
radsimp:cos(x)+i*sin(x)
combine(trig):cos(x)+(-1/2+1/2*cos(2*x))^(1/2)
factor:cos(x)+(-sin(x)^2)^(1/2)
expand:cos(x)+(-sin(x)^2)^(1/2)
convert(exp):
1/2*exp(i*x)+1/2/exp(i*x)+1/4*4^(1/2)*(exp(i*x)-1/exp(i*x))
convert(tan):
(1-tan(1/2*x)^2)/(1+tan(1/2*x)^2)+(-4*tan(1/2*x)^2/(1+tan(1/2*)^2)^2)^(1/2)
ans=cos(x)+i*sin(x)
>> simple(ans) % one more time
convert(exp):exp(i*x)
convert(tan):
(1-tan(1/2*x)^2)/(1+tan(1/2*x)^2)+2*i*tan(1/2*x)/(1+tan(1/2*x)^2)
ans =exp(i*x)
8.9 方程求解
用MATLAB所具有的符号工具可以求解符号方程。
求解单个代数方程。如果表达式不是一个方程式(不含等号),则在求解之前函数solve将表达式置成等于0。
>>solve('a*x^2+b*x+c') % solve for the roots of the quadratic eqution
ans=[1/2/a*(-b+(b^2-4*a*c)^1/2)]
[1/2/a*(-b-(b^2-4*a*c)^1/2)]
结果是符号向量,其元素是方程的2个解。如果想对非缺省x变量求解,solve必须指定变量。
>>solve( 'a*x^2+b*x+c','b') % solve for b
ans=-(a*x^2+c)/x
带有等号的符号方程也可以求解。
>>f=solve('cos(x)=sin(x)') % solve for x
f=1/4*pi
>>t=solve('tan(2*x)=sin(x)')
t=[ 0]
[acos(1/2+1/2*3^(1/2))]
[acos(1/2=1/2*3^(1/2))]
并得到数值解。
>>numeric(f)
ans=0.7854
>>numeric(t)
ans=0
0 + 0.8314i
1.9455
注意在求解周期函数方程时,有无穷多的解。在这种情况下,solve对解的搜索范围限制在接近于零的有限范围,并返回非唯一的解的子集。
如果不能求得符号解,就计算可变精度解。
>>x=solve('exp(x)=tan(x)')
x=1.306326940423079
代数方程组求解
可以同时求解若干代数方程,语句solve(s1,s2,.....,sn)对缺省变量求解n个方程,语句solve(s1,s2,...,sn,' v1,v2,...,vn ')对n个' v1,v2,...vn '的未知数求解n个方程。
单个微分方程
常微分方程有时很难求解,MATLAB提供了功能强大的工具,可以帮助求解微分方程。函数dsovle计算常微分方程的符号解。因为我们要求解微分方程,就需要用一种方法将微分包含在表达式中。所以,dsovle句法与大多数其它函数有一些不同,用字母D来表示求微分,D2,D3等等表示重复求微分,并以此来设定方程。任何D后所跟的字母为因变量。独立变量可以指定或由symvar规则选定为缺省。例如,一阶方程dy/dx=1+y2的通解为:
>>dsolve('Dy=1+y^2') % find the general solution
ans=-tan(-x+C1)
其中,C1是积分常数。求解初值y(0)=1的同一个方程就可产生:
>>dsolve('Dy=1+y^2 ','y(0)=1') % add an initial condition
y=tan(x+1/4*pi)
独立变量可用如下形式指定:
>> dsolve(' Dy=1+y^2 ',' y(0)=1 ',' v ') % find solution to dy/dv
ans=tan(v+1/4*pi)
让我们举一个二阶微分方程的例子,该方程有两个初始条件:
=cos(2x)-y (0)=0 y(0)=1
>> y=dsolve(' D2y=cos(2*x)-y ',' Dy(0)=0 ',' y(0)=1 ')
y=-2/3*cos(x)^2+1/3+4/3*cos(x)
>> y=simple(y) % y looks like it can be simplified
y=-1/3*cos(2*x)+4/3*cos(x)
通常,要求解的微分方程含有一阶以上的项,如:
>>y=solve('D2y-2Dy-3*y=0 ')
y=C1*exp(-x)+C2*exp(3*x)
加上初始条件:y(0)=0和y(1)=1可得到:
>>y=solve('D2y-2Dy-3*y=0','y(0)=0,y(1)=1')
y=1/(exp(-1)-exp(3))*exp(-x)-1/(exp(-1)-exp(3))*exp(3*x)
>>y=simple(y) % this looks like a candidate for simplification
y=-(exp(-x)-exp(3*x))/(exp(3)-exp(-1))
>>pretty(y) % pretty it up
exp(-x)-exp(3 x)
- ---------------------
exp(3) -exp(-1)
绘制区域[-6 2]内的图象为:
>>ezplot(y,[-6 2])
微分方程组
函数dsolve也可同时处理若干个微分方程式,下面有两个线性一阶方程。如
>>[f,g]=dsolve('Df=3*f+4*g','Dg=-4*f+3*g')
f=C1*exp(3*x)*sin(4*x)+C2*exp(3*x)*cos(4*x)
g=-C2*exp(3*x)*sin(4*x)+C1*exp(3*x)*cos(4*x)
加上初始条件:f(0)=0和g(0)=1,可得到:
>>[f,g]=dsolve('Df=3*f+4*g ','Dg=-4*f+3*g ','f(0)=0,g(0)=1 ' )
f=exp(3*x)*sin(4*x)
g=exp(3*x)*cos(4*x)
8.10 线性代数和矩阵
在本节中,我们将介绍符号矩阵和MATLAB提供的工具,它用线性代数求解问题。
符号矩阵
符号矩阵和向量是数组,其元素为符号表达式,可用函数sym来产生。
>>A=sym('[a,b,c;b,c,a;c,a,b]')
A=
[ a, b, c]
[ b, c, a]
[ c, a, b]
>>G=sym('[cos(t),sin(t);-sin(t),cos(t)]')
G=
[ cos(t), sin(t)]
[-sin(t), cos(t)]
函数sym也可以扩展成定义各元素的公式。注意,只有在这种情况下,i,j分别表示行列的位置;且不影响i,j的缺省值。下面的例子建立了3×3的矩阵,其元素依赖于行和列的位置。
>> S=sym(3,3,'(i+j)+(i-j+s)') % create a matrix using a formula
S=
[ 2/s, 3/(-1+s), 4/(-2+s)]
[1/(1+s), 4/s, s/(-1+s)]
[4/(2+s), 5/(1+s), 6/s]
>>S=sym(3,3,'m ','n ','(m-n)/(m-n-t) ') % use m and n in another formula
S=
[ 0, -1/(-1-t), -2/(-2-t)]
[1/(1-t), 0, -1/(-1-t)]
[2/(2-t), 1/(1-t), 0]
函数sym也可以把数值矩阵转换成符号形式
>> M=[1.1,1.2,1.3;2.1,2.2,2.3;3.1,3.2,3.3] % a numeric matrix
M=
1.1000 1.2000 1.3000
2.1000 2.2000 2.3000
3.1000 3.2000 3.3000
>> S=sym(M) %convert to symbolic form
S=
[11/10, 6/5, 13/10]
[21/10, 11/5, 23/10]
[31/10, 16/5, 33/10]
如果数值矩阵的元素可以指定为小的整数之比,则函数sym将采用有理分式表示。如果元素是无理数,则在符号形式中sym将用符号浮点数表示元素。
>>E=[exp(1) sqrt(2)]
E=2.7183 1.4142
>>sym(E)
ans=[3060513257434036*2^(-50), 3184525836262886*2^(-51)]
用函数symsize可以得到符号矩阵的大小(行,列数)。函数返回数值或向量,而不是符号表达式。symsize的四种形式说明如下:
>>S=sym('[a,b,c;d,e,f]') % create a symbolic matrix
S= [a,b,c]
[d,e,f]
>>d=symsize(S) % return the size of S as the 2-element vector d
d=2 3
>>[m,n]=symsize(S) % return the number of rows in m,and column in n
m=2
n=3
>>m=symsize(S,1) % return the number of rows
m=2
>>n=symsize(S,2) % return the number of columns
n=3
数值数组用N(m,n)形式来访问单个元素,但符号数组元素必须用函数如sym(S,m,n)来获取。若能用相同的句法当然好,但MATLAB符号表达式的表示不方便。在内部,符号数组表示成一个字符串数组;而S(m,n)返回单个字母。所以,符号数组的各个元素,必须由符号函数,如sym,来指定,而不是直接指定。
>>G=sym(' [ab,cd;ef,gh] ') % create a 2-by-2 symbolic matrix
G= [ab, cd]
[ef , gh]
>>G(1,2) % this is the second character of the first row of G
ans=a
>>r=sym(G,1,2) % this is the second expressinon in the first row of G
r=cd
记住,在上例中的符号矩阵G,实际是以2×7字符数组存储在计算机中。第一行为'[ab,cd]',所以第二个元素是'a '。
最后,sym可以用于改变符号数组的一个元素。
>>sym(G,2,2,'pq') % change the (2,2) element in G from 'gh 'to'pq'
ans=
[ab, cd]
[ef, pq]
代数运算
用函数symadd,symsub,symmul和symdiv,对符号矩阵可以执行许多通用的代数运算,用sympow可计算乘幂,用transpose计算符号矩阵的转置。
>>G=sym('[cos(t),sin(t0;-sin(t0,cos(t)]') %create a symbolic matrix
>>symadd(G, 't') %add't'to each element
ans=
[ cos(t)+t, sin(t)+t]
[ -sin(t)+t, cos(t)+t]
>>symmul(G,G) % multiply G by G; Synpow(G,2) does the same thing
ans=
[cos(t)^2-sint(t)^2, 2*cos(t)*sin(t)]
[ -2*cos(t)* sin(t), cos(t)^2-sin(t)^2]
>> simple(G) % try to simplify
ans=
[cos(2*t), sin(2*t)]
[-sin(2*t), cos(2*t)]
8.11 线性代数运算
用函数inverse和determ,可计算符号矩阵的逆阵以及行列式。
>>H=sym(hilb(3)) %the symbolic form of the numeric 3-by-3 Hilbert matrix
H=[1,1/2,1/3]
[1/2,1/3,1/4]
[1/3,1/4,1/5]
>>determ(H) % dind the determinant of H
ans=1/2160
>>J=inverse(H) % find the inverse of H
J=[ 9, -36, 30]
[-36, 192,-180]
[ 30,-180, 180]
>>determ(J) % find the determinant of the inverse
ans=2160
用函数linsolve求解齐次线性方程;这是等价于基本的MATLAB的逆斜杠\算子的符号,linsolve(A,B)对X方阵求解矩阵方程A*X=B。 symop将其参量串接起来,并计算所得到的表达式。
>>f='cos(x)' % create an expression
f=cos(x)
>>symop('atan(',f,'+ ',a,')','^2')
ans=atan(cos(x)+a)^2
在用函数symop时,若将数组和标量混合应慎重。例如:
>>M=sym('[a b;c d]');
>>symop(M,'+','t')
ans=
[a+t, b]
[ c, d+t]
是把t加到M的对角线上。
函数charpoly求解矩阵的特征多项式。
>>G=sym('[1,1/2;1/3,1/4] ') % create a symbolic matrix
G=[ 1, 1/2]
[1/3, 1/4]
>>charpoly(G) % find the characteristic polynomial of G
ans=x ' 2-5/4*x+1/12
用函数eigensys可以求得符号矩阵的特征根和特征向量
>> F=sym(' [1/2,1/4;1/4,1/2] ') % create a symbolic matrix
F=[1/2,1/4]
[1/4,1/2]
>>eigensys(F) % find the eigenvalues of F
ans=[3/4]
[1/4]
>>[V,E]=eigensys(F) %find eigenvalues E and eigenvectors v
V=[-1, 1]
[ 1, 1]
E=[1/4]
[3/4]
矩阵的约当(Jordan)标准型是特征值的对角矩阵;转换矩阵的列是特征向量。对于给定的矩阵A, jordan(A)求出非奇异的矩阵V,使得inv(V)*A*V成为约当标准型,函数jordan有两种形式
>>jordan(F) %find the Jordan form of F,above
ans=[1/4,0]
[0,3/4]
>>[V,J]=jordan(F) % find the Jordan form and eigenvectors
V=[ 1/2,1/2]
[-1/2,1/2]
J=[1/4, 0]
[ 0, 3/4]
标准形和特征向V的列是某些F可能的特征向量。
因为F是非奇异的,F的零空间的基是空矩阵,而列空间的基是单位阵。
>>F=sym('[1/2,1/4;1/4,1/2]') %recreate F
F=[1/2, 1/4]
[1/4, 1/2]
>>nullspace(F) % the nullspace of F is the empty matrix
ans=[]
>>colspace(F) % find the column space of F
ans=[1, 0]
[0, 1]
用函数singvals可求解矩阵奇异值。
>> A=sym(magic(3)) % Generate a 3-by-3 matrix
A= [8, 1, 6]
[3, 5, 7]
[4, 9, 2]
>> singvals(A) % find the singular expressions
ans=[15.00000000000000]
[6.928203230275511]
[3.464101615137752]
函数jacobian(w,v)可相对于v求解w的雅可比(Jacobia)值。其结果的第(i,j)项是df(i)/dv(j)。注意,当f是标量时,f的雅可比值是f的梯度。
>>jacobian( ' u*exp(v) ' ,sym( ' u,v ' ))
ans=[exp(v,u*exp(v))
符号表达式运算小结如下表:
numeric 符号到数值的转换
pretty 显示悦目的符号输出
subs 替代子表达式
sym 建立符号矩阵或表达式
symadd 符号加法
symdiv 符号除法
symmul 符号乘法
symop 符号运算
sympow 符号表达式的幂运算
symrat 有理近似
symsub 符号减法
symvar 求符号变量
collect 合并同类项
expand 展开
factor 因式
simple 求解最简形式
simplify 简化
symsum 和级数
charpoly 特征多项式
horner 嵌套多项式表示
numden 分子或分母的提取
poly2sym 多项式向量到符号的转换
sym2poly 符号到多项式向量的转换
diff 微分
int 积分
jordan 约当标准形
taylor 泰勒级数展开
digits 设置可变精度
vpa 可变精度计算
compose 函数的复合
dsolve 微分方程的求解
finverse 函数逆
linsolve 齐次线性方程组的求解
solve 代数方程的求解
charploy 特征多项式
determ 矩阵行列式的值
eigensys 特征值和特征向量
inverse 矩阵逆
jordan 约当标准形
linsolve 齐次线性方程组的解
transpose 矩阵的转置


copyright© 2005 西南科技大学网络学院版权所有
制作:西南科技大学理学院数学系
Designed by xiandaquan