MATLAB的符号运算.doc

MATLAB的符号运算
前面介绍的内容基本上是MATLAB的数值计算功能,参与运算过程的变量都是被赋了值的数值变量.在MATLAB环境下,符号运算是指参与运算的变量都是符号变量,即使是数字也认为是符号变量. 数值变量和符号变量是不同的.
1 符号微积分
下面着重介绍一些与微积分有关的指令,这些指令都需要符号表达式作为输入宗量.
求和
symsum(S) 对通项S求和,其中k为变量且从0变到k-1.
symsum(S,v) 对通项S求和,指定其中v为变量且v从0变到v-1.
symsum(S,a,b) 对通项S求和其中k为变量且从a变到b.
symsum(S,v,a,b) 对通项S求和,指定其中v为变量且v从a变到b.
例：求,键入
k=sym('k') % k是一个符号变量;
symsum(k)
得 ans = 1/2*k^2-1/2*k
例：求,键入:
symsum(k^2,0,10)
得 ans = 385
例：求键入
symsum('x'^k/sym('k!'),k,0,inf),
得 ans = exp(x)
这最后的一个例子是无穷项求和.
求极限
limit(P) 表达式P中自变量趋于零时的极限
limit(P,a) 表达式P中自变量趋于a时的极限
limit(P,x,a,'left') 表达式P中自变量x趋于a时的左极限
limit(P,x,a,'right') 表达式P中自变量x趋于a时的右极限
例：求,键入
P=sym('sin(x)/x');
limit(P)
得 ans = 1
例：求 键入
P=sym('1/x');
limit(P,'x',0,'right')
得 ans = inf
例：求,键入:
P=sym('(sin(x+h)-sin(x))/h');h=sym('h');
limit(P,h,0)
得ans = cos(x)
例：求, 键入
v=sym('[(1+a/x)^x,exp(-x)]');
limit(v,'x',inf,'left')
得 ans = [ exp(a), 0]
求导数
diff(S,v) 求表达式S对变量v的一阶导数.
diff(S,v,n) 求表达式S对变量v的n阶导数.
例如：设A=,求键入命令:
syms a b x;
A= [1/(1+a),(b+x)/cos(x);1,exp(x^2)];
diff(A,'x')
得 ans = [0 1/cos(x)+(b+x)/cos(x)^2*sin(x)]
[0, 2*x*exp(x^2)]
例：求y=sinx+ex的三阶导数,键入命令:
diff('sin(x)+x*exp(x)',3)
得 ans = -cos(x)+3*exp(x)+x*exp(x)
例：设,求A的先对x再对y的混合偏导数.可键入命令：
S=sym('[x*sin(y),x^n+y;1/x/y,exp(i*x*y)]');
dsdxdy=diff(diff(S,'x'),'y')
得: dsdxdy = [ cos(y), 0]
[ 1/x^2/y^2, i*exp(i*x*y)-y*x*exp(i*x*y)]
例：求y=(lnx)x的导数.可键入命令：
p='(log(x))^x';
p1=diff(p,'x')
得：p1 = log(x)^x*(log(log(x))+1/log(x))
例：求y=xf(x2)的导数.可键入命令：
p='x*f(x^2)';p1=diff(p,'x')
得：p1 = f(x^2)+2*x^2*D(f)(x^2)
例：求xy=ex+y的导数.可键入命令：
p='x*y(x)-exp(x+y(x))';p1=diff(p,'x')
得：p1 = y(x)+x*diff(y(x),x)-(1+diff(y(x),x))*exp(x+y(x))
再键入
p2='y+x*dy-(1+dy)*exp(x+y)=0';
dy=solve(p2,'dy')%把dy作为变量解方程
得 dy= -(y-exp(x+y))/(x-exp(x+y))
求Taylor展开式
taylor(f,v) f对v的五阶Maclaurin展开.
taylor(f,v,n) f对v的n-1阶Maclaurin展开.
例：求sinxe-x 的7阶Maclaurin展开.可键入
f=sym('sin(x)*exp(-x)');F=taylor(f,8) 得
F = x-x^2+1/3*x^3-1/30*x^5+1/90*x^6-1/630*x^7
例：求sinxe-x 在x=1 处的7阶Taylor展开.可键入
f=sym('sin(x)*exp(-x)');F=taylor(f,8,1) 得
F = sin(1)*exp(-1)+(-sin(1)*exp(-1)+cos(1)*exp(-1))*(x-1)
-cos(1)*exp(-1)x-1)^2
+(1/3*sin(1)*exp(-1)+1/3*cos(1)*exp(-1))*(x-1)^3
-1/6*sin(1)*exp(-1)*(x-1)^4
+(1/30*sin(1)*exp(-1)-1/30*cos(1)*exp(-1))*(x-1)^5
+1/90*cos(1)*exp(-1)*(x-1)^6
+(-1/630*cos(1)*exp(-1)-1/630*sin(1)*exp(-1))*(x-1)^7
多元函数的Taylor展开
MATLAB不能直接进行多元函数的Taylor展开.必须先调用MAPLE函数库中的mtaylor命令.方法为：
在MATLAB的工作窗口中键入
maple('readlib(mtaylor)')
mtaylor的格式为
mtaylor(f,v,n)
f为欲展开的函数式
v为变量名.写成向量的形式：[var1=p1,var2=p2,…,varn=pn],展开式将在（p1,p2,…,pn）处进行.如只有变量名,将在0点处展开.n为展开式的阶数（n-1阶）.要完成Taylor展开,只需键入maple('mtaylor（f,v,n）')即可.
例：在（x0,y0,z0）处将F=sinxyz进行2阶Taylor展开.键入
syms x0 y0 z0
maple('readlib(mtaylor)');
maple('mtaylor(sin(x*y*z),[x=x0,y=y0,z=z0],2)') 得：
ans = sin(x0*y0*z0)+cos(x0*y0*z0)*y0*z0*(x-x0)+cos(x0*y0*z0)*x0*z0*(y-y0)
+cos(x0*y0*z0)*x0*y0*(z-z0)
求积分
int(P) 对表达式P进行不定积分.
int(P,v) 以v为积分变量对P进行不定积分.
int(P,v,a,b) 以v为积分变量,以a为下限,b为上限对P进行定积分.
例：求,可键入
int('-2*x/(1+x^2)^2')
得 ans = 1/(1+x^2)
例：求,可键入
键入int('x/(1+z^2)','z')
得 ans = atan(z)*x
例：求,可键入
int('x*log(1+x)',0,1) 得
ans = 1/4
例：求可键入：
int('2*x','sin(t)','log(t)') 得：
ans = log(t)^2-sin(t)^2
对（符号）矩阵积分
例：求,输入
int('[exp(t),exp(a*t)]'),得：
ans = [ exp(t), 1/a*exp(a*t)]
求符号方程的解
ⅰ线性方程组的求解
线性方程组的形式为A*X=B；其中A至少行满秩.
X=linsolve(A,B) 输出方程的特解X.
例：解方程组.键入
A=sym('[cos(t),sin(t);sin(t),cos(t)]');B=sym('[1;1]');
c=linsolve(A,B)
c =[ 1/(sin(t)+cos(t))]
[ 1/(sin(t)+cos(t))]
ⅱ 代数方程的求解
solve(P,v)
对方程P中的指定变量v求解.v可省略.
solvep1,P2,…,Pn,v1,v2,…,vn)
对方程P1,P2,…Pn中的指定变量v1, v2…vn求解.
例：解,可输入
solve('p+sin(x)=r') 得：
ans =-asin(p-r)
例：解,可输入：
P1='x^2+x*y+y=3';P2='x^2-4*x+3=0';
[x,y]=solve(P1,P2) 得：
x = [ 1]
[ 3]
y = [ 1]
[ -3/2]
解,可输入：
P1='a+u^2+v^2=0';P2='u-v=1';[u,v]=solve(P1,P2,'u','v') 得：

u = [ 1/2+1/2*(-1-2*a)^(1/2)]
[ 1/2-1/2*(-1-2*a)^(1/2)]
v = [ -1/2+1/2*(-1-2*a)^(1/2)]
[ -1/2-1/2*(-1-2*a)^(1/2)]
对于有些无法求出解析解的非线性方程组,MATLAB只给出一个数值解.这一点可以从表示解的数字不被方括号括住而确定.
例：解键入：
[x,y]=solve('sin(x+y)-exp(x)*y=0','x^2-y=2') 得：
x = -6.0173272500593065641097297117905
y = 34.208227234306296508646214438330
由于这两个数字没有被[ ]括住,所以它们是数值解.
另外,可利用solve来解线性方程组的通解.
例：解键入
P1='2*x1+7*x2+3*x3+x4=6'; P2='3*x1+5*x2+2*x3+2*x4=4';
P3='9*x1+4*x2+x3+7*x4=2';
u=solve(P1,P2,P3,'x1','x2','x3','x4')
Warning: 3 equations in 4 variables.
u = x1: [1x1 sym]
x2: [1x1 sym]
x3 [1x1 sym]
x4: [1x1 sym]
可以看到：屏幕提示“有3个方程4个变量”,意为解不唯一（有时会提示解不唯一）.且输出的是解的结构形式.为进一步得到解,可输入：
u.x1,u.x2,u.x3,u.x4, 得：
ans = x1
ans = -5*x1-4*x4
ans = 11*x1+9*x4+2
ans = x4
这样就得到了原方程组的通解.
⑷ 解符号微分方程
解符号微分方程的命令格式为： dsolve('eq1','eq2',…).
其中eq表示相互独立的常微分方程、初始条件或指定的自变量.默认的自变量为t.如果输入的初始条件少于方程的个数,则在输出结果中出现常数c1,c2等字符.关于微分方程的表达式有如下的约定：字母y表式函数,Dy表示y对t的一阶导数；Dny表示y对t的n阶导数.
例如：求
的解可键入：[x,y]=dsolve('Dx=y','Dy=-x') 得
x =cos(t)*C1+sin(t)*C2
y =-sin(t)*C1+cos(t)*C2
dsolve中的输入宗量最多只能有12个,但这并不妨碍解具有多个方程的方程组,因为可以把多个方程或初始条件定义为一个符号变量进行输入.
例如求 , , f(0)=0 , g(0)=1 的解.可输入指令：
P='Df=3*f+4*g,Dg=-4*f+3*g';v='f(0)=0,g(0)=1';
[f,g]=dsolve(P,v)
f = exp(3*t)*sin(4*t)
g = exp(3*t)*cos(4*t)
注意：微分方程表达式中字母D必须大写.
例如求解微分方程
可输入
y=dsolve('D3y=-y','y(0)=1,Dy(0)=0,D2y(0)=0','x') 得:
y = (1/3+2/3*exp(1/2*x)*cos(1/2*3^(1/2)*x)*exp(x))/exp(x)
最后看一个解非线性微分方程的例子：
dsolve('(Dy)^2+y^2=1','y(0)=0','x')
ans = [ sin(x)]
[ -sin(x)]
对于无法求出解析解的非线性微分方程,屏幕将提示出错信息.
微分方程的数值解及其它问题的数值解
ⅰ 常微分方程的数值解
MATLAB提供了求微分方程数值解的指令：
[t,x]=ode23('fname',[t0,tf],x0,tol,trace)
[t,x]=ode45('fname',[t0,tf],x0,tol,trace)
这两个格式中的输入参数意义完全一样.下面介绍这两个格式的有关内容及各参数的意义.这两个格式都采用Runge--Kutta法求解微分方程的数值解.它们是针对一阶微分方程组设计的.因此,如果待解的是高阶微分方程,那么首先要化成形式为x'=f(t,x)的一阶微分方程组.称为“状态方程”.
‘fname’是f(t,x)的函数名.该函数以x'为输出,以t,x为输入变量,注意次序不能颠倒.
t0和tf分别是积分的起始值和终止值.x0是初始值,以向量的形式输入.tol是用来控制精度的参数,可缺省.缺省时ode23默认tol=1.e-3;ode45默认tol=1.e-6.trace用来控制是否显示中间结果,可缺省.缺省时,默认trace=0,不显示.输出结果t和x分别是时间向量和相应的状态向量.虽然ode45比ode23的精度高,但它的运算速度更快.
例：求著名的Van der pol方程,并绘出其解的图形.
第一步：在编辑器中编写名为fname的M文件.
function X=fname(t,x)
X=zeros(2,1);
X(1)=(1-x(2)^2)*x(1)-x(2);
X(2)=x(1);
第二步：将此文件存放于自己的文件夹中听候调用.
第三步：在MATLAB的命令窗口调用这个函数,即键入如下命令：
[t,x]=ode45('fname',[0,20],[0,0.5]);
plot(t,x)
ⅱ 数值积分
quad('fname',a,b,tol,trace) Simpson法求数值积分.
quad8('fname',a,b,tol,trace) Newton-Cotes法求数值积分.
fname是被积函数文件名b,a分别是积分上下限用tol来控制积分精度.可缺省.缺省时默认tol=0.001.用trace来控制是否用图形显示积分过程.可缺省.缺省时默认trace=0,不显示图形.
例如：求 dx
第一步：在编辑器中建立被积函数的M文件.取名为fname即在编辑器中输入：
function y=fname(x)
y=exp(-x^2);
第二步：将此文件存放于自己的文件夹中.
第三步：在MATLAB环境下调用fname.即输入
s=quad8('fname',0,3)
就可以得到结果：
s =8862