1、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
2、(S,v,a,b) 对通项 S 求和,指定其中 v 为变量且 v 从 a 变到 b.例：求 ,键入10kik=sym(k) % k 是一个符号变量;symsum(k) 得 ans = 1/2*k2-1/2*k 例：求 ,键入:102ksymsum(k2,0,10)得 ans = 385 例：求 键入0!kxsymsum(xk/sym(k!),k,0,inf),得 ans = exp(x) 这最后的一个例子是无穷项求和.求极限limit(P) 表达式 P 中自变量趋于零时的极限limit(P,a) 表达式 P 中自变量趋于 a 时的极限limit(P,x,a,left) 表达式 P 中自变量 x
3、 趋于 a 时的左极限limit(P,x,a,right) 表达式 P 中自变量 x 趋于 a 时的右极限例：求 ,键入xsinlm0P=sym(sin(x)/x);limit(P) 得 ans = 1 例：求 键入x1lim0P=sym(1/x);limit(P,x,0,right) 得 ans = inf 例：求 ,键入:hxhsin)si(l0P=sym(sin(x+h)-sin(x)/h);h=sym(h);limit(P,h,0) 得 ans = cos(x) 例：求 , 键入)lim ,)1(li-xxx eav=sym(1+a/x)x,exp(-x);limit(v,x,inf,
4、left) 得 ans = exp(a), 0 求导数diff(S,v) 求表达式 S 对变量 v 的一阶导数.diff(S,v,n) 求表达式 S 对变量 v 的 n 阶导数.例如：设 A= ,求 键入命令:21cosxebadAsyms a b x;A= 1/(1+a),(b+x)/cos(x);1,exp(x2);diff(A,x) 得 ans = 0, 1/cos(x)+(b+x)/cos(x)2*sin(x)0, 2*x*exp(x2) 例：求 y=sinx+ex 的三阶导数,键入命令:diff(sin(x)+x*exp(x),3) 得 ans = -cos(x)+3*exp(x)+
5、x*exp(x) 例：设 ,求 A 的先对 x 再对 y 的混合偏导数.可键入命令：xyineA1siS=sym(x*sin(y),xn+y;1/x/y,exp(i*x*y);dsdxdy=diff(diff(S,x),y) 得: dsdxdy = cos(y), 0 1/x2/y2, 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(x2);p1=diff(p,x)得：p
6、1 = f(x2)+2*x2*D(f)(x2)例：求 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 阶 Macla
7、urin 展开.例：求 sinxe-x 的 7 阶 Maclaurin 展开.可键入f=sym(sin(x)*exp(-x);F=taylor(f,8) 得F = x-x2+1/3*x3-1/30*x5+1/90*x6-1/630*x7 例：求 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*co
8、s(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 的格式为m
9、taylor(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 z0maple(readlib(mtaylor);maple(mtaylor(sin(x*y*z),x=x0,y=y0,z=z0,2) 得：ans = sin(x0*y0
10、*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 进行定积分.例：求 ,可键入dx2)1(int(-2*x/(1+x2)2)得 ans = 1/(1+x2) 例：求 ,可键入dzx)1(2键入 int(x/(1+z2),z) 得 ans = atan(z)*x 例：求 ,可键入10)ln(dx
11、int(x*log(1+x),0,1) 得ans = 1/4 例：求 可键入：txdlnsi2int(2*x,sin(t),log(t) 得：ans = log(t)2-sin(t)2 对（符号）矩阵积分例：求 ,输入dtetaint(exp(t),exp(a*t),得： ans = exp(t), 1/a*exp(a*t) 求符号方程的解线性方程组的求解线性方程组的形式为 A*X=B；其中 A 至少行满秩.X=linsolve(A,B) 输出方程的特解 X.例：解方程组 .键入1cosintA=sym(cos(t),sin(t);sin(t),cos(t);B=sym(1;1);c=lins
12、olve(A,B) c = 1/(sin(t)+cos(t) 1/(sin(t)+cos(t) 代数方程的求解solve(P,v)对方程 P 中的指定变量 v 求解.v 可省略.solve(p1,P2,Pn,v1,v2,vn) 对方程 P1,P2,Pn 中的指定变量 v1, v2vn 求解.例：解 ,可输入rxpsinsolve(p+sin(x)=r) 得：ans =-asin(p-r) 例：解 ,可输入：0342xyP1=x2+x*y+y=3;P2=x2-4*x+3=0; x,y=solve(P1,P2) 得：x = 1 3y = 1 -3/2 解 ,可输入：102vuaP1=a+u2+v2
13、=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 只给出一个数值解.这一点可以从表示解的数字不被方括号括住而确定.例：解 键入：20)sin(yxexx,y=solve(sin(x+y)-exp(x)*y=0,x2-y=2) 得：x = -6.0173272500593065641097297117905y = 34.2082272
14、34306296508646214438330 由于这两个数字没有被 括住,所以它们是数值解.另外,可利用 solve 来解线性方程组的通解.例：解 键入2467149532XP1=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 symx2: 1x1 symx3: 1x1 symx4: 1x1 sym 可以看到：屏幕提示“有 3 个方程 4 个变量”,
15、意为解不唯一（有时会提示解不唯一）.且输出的是解的结构形式.为进一步得到解,可输入：u.x1,u.x2,u.x3,u.x4, 得：ans = x1ans = -5*x1-4*x4ans = 11*x1+9*x4+2ans = x4 这样就得到了原方程组的通解. 解符号微分方程解符号微分方程的命令格式为： dsolve(eq1,eq2, ).其中 eq 表示相互独立的常微分方程、初始条件或指定的自变量.默认的自变量为 t.如果输入的初始条件少于方程的个数,则在输出结果中出现常数 c1,c2 等字符.关于微分方程的表达式有如下的约定：字母 y 表式函数,Dy 表示 y 对 t 的一阶导数；Dny
16、表示 y 对 t 的 n阶导数.例如：求 xdty的解可键入：x,y=dsolve(Dx=y,Dy=-x) 得x =cos(t)*C1+sin(t)*C2y =-sin(t)*C1+cos(t)*C2 dsolve 中的输入宗量最多只能有 12 个,但这并不妨碍解具有多个方程的方程组,因为可以把多个方程或初始条件定义为一个符号变量进行输入.例如求 , , f(0)=0 , g(0)=1 的解.可输入指令：gfdt43gfdt3P=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
17、(3*t)*cos(4*t) 注意：微分方程表达式中字母 D 必须大写.例如求解微分方程 0()y,(0)1,y()3dx可输入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+y2=1,y(0)=0,x) ans = sin(x) -sin(x) 对于无法求出解析解的非线性微分方程,屏幕将提示出错信息.微分方程的数值解及其它问题的数值解 常微分方程的数值解MATLAB 提供了求微分方
18、程数值解的指令：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 是初始值,
19、以向量的形式输入.tol 是用来控制精度的参数,可缺省.缺省时 ode23 默认 tol=1.e-3;ode45 默认 tol=1.e-6.trace 用来控制是否显示中间结果,可缺省.缺省时,默认 trace=0,不显示.输出结果 t 和 x 分别是时间向量和相应的状态向量.虽然 ode45 比 ode23 的精度高,但它的运算速度更快.例：求著名的 Van der pol 方程 ,并绘出其解的图形 .xyy)1(2第一步：在编辑器中编写名为 fname 的 M 文件.function X=fname(t,x)X=zeros(2,1);X(1)=(1-x(2)2)*x(1)-x(2);X(2
20、)=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,不显示图形.例如：求 dx302xe第一步：在编辑器中建立被积函数的 M 文件.取名为 fname 即在编辑器中输入：function y=fname(x)y=exp(-x2);第二步：将此文件存放于自己的文件夹中.第三步：在 MATLAB 环境下调用 fname.即输入s=quad8(fname,0,3)就可以得到结果：s =8862
