MATLAB-ch013(数值计算-微积分.ppt

第第1313讲讲 数值计算数值计算 微积分微积分 张建瓴 13.1 数值积分 在工程教学和应用中，除了进行数据逼近外，还要求逼近 曲线下面的面积，这就是积分问题。 一、数值积分方法 典型的数值积分方法有：用常数(0阶多项式)近似函 数矩形法；用直线(一阶多项式)近似函数曲线的梯形 法；用抛物线(二阶多项式)近似函数曲线的Simpson 法，以及用一般多项式近似函数的Romberg法等。 y=sin(x3)*sqrt(x) x求y, 表13-1列出了函数数值积分的一些命令。 表13-1 函数数值积分的命令 常见的一元数值积分命令 MATLAT提供了在有限区间内，数值计算某函数下的面积（ 积分）的三种函数：trapz，quad和quad8。 二、一（元）维数值积分 1、trapz函数 函数trapz通过计算若干梯 形面积的和来近似某函数的 积分，这些梯形如图13-1所 示，是通过使用函数humps 的数据点形成。 图13-2 较好的梯形逼近曲线下的面积示意图 从图中可明显地看出，单个梯形的面积在某一段欠估计了 函数真正的面积，而在其它段又过估计了函数的真正面积 。如同线性插值，当梯形数目越多时，函数的近似面积越 准确。例如，在图13-1中，如果我们大致增加一倍数目的 梯形，我们得到如下（如图13-2）所示的更好的近似结果 。 trapz 函数 对如上所示的两个曲线，用trapz在区间-1x=-1: 0.17: 2; % rough approximation y=humps(x); area=trapz(x, y) % call trapz just like the plot command area = 25.9174 trapz 函数的应用 x=-1: 0.07: 2; % better approximation y=humps(x); area=trapz(x , y) area = 26.6243 上述两个结果不同是基于对图形的观察，粗略近似可能低 估了实际面积。除非特别精确，没有准则说明哪种近似效 果更好。 trapz 函数的应用 MATLAB提供的求积函数命令quad和quad8在使用时， 其递推的层次限制在十层以内，达到这个限制则会 提示警告信息，并且这两个函数命令都不能解决可 积的奇异值问题，例如，求解 。 quad函数和quad8函数 函数quad和quad8完整的调用格式为： （1）q=quad(fun,a,b,tol,trace,pl,p2,) 采用 Simpson法计算积分； （2）q=quad8(fun,a,b,tol,trace,p1,p2,) 采用 八样条Newton-Cotes公式求数值积分。 其中：fun是被积函数，可以是表达式字符串、内联函数 、M函数文件名，被积函数的自变量，一般采用字母x； a、b分别是积分的上、下限，都是确定的值； quad和quad8函数调用格式 tol是一个二元向量，它的第一个元素用来控制相对误差 ，第二个元素用来控制绝对误差，缺省时积分的相对精度 为0.001； trace如果取非零值时，将以动态图形的形式展现积分的 整个过程，若取零值，则不画图，其缺省值是0； pl和p2是向被积函数传递的参数。 在上面的调用格式中，前三个输入参数是调用时必须的， 而后面的输入参数可缺省。 quad和quad8的参数 MATLAB的函数quad和quad8是基于数学上的正方形概念来 计算函数的面积。为获得更准确的结果，两个函数在所需 的区间都要计算被积函数。 与简单的梯形比较，这两个函数进行更高阶的近似，而且 quad8比quad更精确。这两个函数的调用方法与fzero相同 ，即 area=quad(humps,-1,2) % find area between -1 and 2 area = 26.3450 quad和quad8函数的调用 area=quad8(humps,-1,2) area = 26.3450 注意： 这两个函数返回完全相同的估计面积，而且这个估计值在 两个trapz面积的估计值之间。 quad和quad8函数的调用 求函数的数值积分 （1）建立函数funq function y=funq(x) y=x.3+x.2+2; （2）对被积函数funq进行数值积分 q=quad(funq,-1,1,le-4) 使用quad命令求数值积分 q= 4.6667 例13-1 example13_1.m q8=quad8(funq,-1,1,le-4,1) 用quad8命令求数值 积分 q8= 4.6667 程序的运行结果显示出积分的过程如图13-3所示。 例13-1 一元函数积分中存在的问题，同样存在于多重积分中。 1、积分限为常数的二重积分 多重积分使用函数dblquad，其完整的调用格式为： result=dblquad(fun,inmin, inmax, outmin, outmax, tol, method) 其中：输入参数fun是被积函数，可以直接用字符串内联 函数或M函数文件表达，但不论什么形式，该被积函数应 有两个变量，即内变量和外变量。内变量接受向量输入， 外变量接受标量输入。被积函数的输出是与内变量同长的 向量。 三、多重数值积分 输入参数inmin，inmax是内变量的下限和上限； outmin、outmax是外变量的下限和上限； tol的含义与命令quad中的情况相同； method是积分方法选项，如“quad”和“quad8”等。 注意： 该命令不适用于内积分区间上、下限为函数的情况。 dblquad函数的参数 求积分上下限都为常数的二重积分，被积函数为 y*sin(x)+s*cos(y)，其中x的取值范围是到2，y的 取值范围是0到。 （1）建立名为integrnd的M文件 fimction out=integrnd(x,y) out=y*sin(x)+x*cos(y) （2）用函数dblquad命令来求integrnd的二重积分 result=dblquad(integrnd,pi,2*pi,0,pi) result= -9.8698 例13-6 example13_6.m 对于内积分上下限是外积分变量的函数的积分问题 ，求解过程较为麻烦。 一般方法都是先求出 ，然 后再求 。 2、内积分上下限为函数的二重积分 【例13-10】计算。 （1）编写内积分区间上下限的M函数文件 x_low.m function f=x_low(y) f=sqrt(y); （2）编写被积函数 被积函数函数采用内联函数 ff=inline(x.2+y.2,x,y)表示。 （3）被积函数用内联函数表达时，运行以下指 令，即得结果 ff=inline(x.2+y.2,x,y);%被 积函数的内联函数表达 SS=double_int(ff,x_low,2,1,4) SS = 9.5810 （4）本题用符号计算很容易获得高精度解 Ssym=vpa(int(int(x2+y2,x,sqrt(y),2),y,1,4) Ssym = 9.580952380952381 例13-10 example13_10.m 1、“完整”离散序列的数值卷积 （1）求和运算上下界的确定 （2）卷积C(n)“非平凡”区间的确定 （3）“截尾”离散序列的数值卷积 （4）多项式乘法与离散卷积的算法同构 （5）连续时间函数的数值卷积 （6）“零阶”近似算法 （7）“积分-插值”混合算法 （8）SIMULINK卷积法 【例13-12】example12_12.m 有序列和。 （A）求这两个完整序列的卷积，并图示。 （B）假设A(7)及其后的4个非零值未知，而成为 截尾序列，求卷积并图示。 %完整序列卷积 a=ones(1,10);n1=3;n2=12; %完整A(n)序列的非平凡 值和区间端点 b=ones(1,8);n3=2;n4=9; %完整B(n)序列的非平凡 值和区间端点 c=conv(a,b);nc1=n1+n3;nc2=n2+n4; %计算卷积和确 定卷积非平凡区间端点 kc=nc1:nc2; %构成非平凡区间的序号自变量 % 截尾序列卷积 aa=a(1:6);nn1=3;nn2=8;%截 尾A(n)序列的非平凡值和区间端点 cc=conv(aa,b);ncc1=nn1+n3; 图13-7 “完整”序列卷积和“截尾”序列卷积 nx=nn2+n4; %“非平凡”区间右端点 ncc2=min(nn1+n4,nn2+n3);%截尾序列卷积被正确计 算区间的右端点 kx=(ncc2+1):nx;kcc=ncc1:ncc2;N=length(kcc); stem(kcc,cc(1:N),r,filled) axis(nc1-2,nc2+2,0,10),grid,hold on stem(kc,c,b),stem(kx,cc(N+1:end),g),hold off 【例13-8】求函数u(t)=e-tU(t)和h(t)=te- t/2U(t)的卷积。本例展示： （A）符号Laplace变换求卷积的理论表示； （B）SIMULINK卷积法的执行过程和它的快速精 确性。 （C）从理论符号解产生相应的理论数值序列。 （1）符号卷积法（得解析结果） syms tao;t=sym(t,positive); % 把t定义为“取正”符号变量 US1=laplace(exp(-t); %u(t)的L氏变换 HS1=laplace(t*exp(-t/2) %h(t)的L氏变换 yt1=simple(ilaplace(US1*HS1) %L氏反变换得卷 积的理论解 图13-8 卷积解算模型exm5835_2_2.mdl HS1 = 1/(1/2+s)2 yt1 = 4*exp(-t)+(2*t-4)*exp(-1/2*t) （2）SIMULINK卷积法 根据h(t)的传递函数，构作SIMULINK模型 exm5835_2.mdl，并运行，可从示波器上看到卷积结 果。 （3）比较两种卷积方法的结果，并用图形表示 t=yt2(:,1); %exm5835_2.mdl示波器所保存的仿真采样时 间序列 yyt1=eval(vectorize(char(yt1);%理论解的 数值序列 dy,kd=max(abs(yyt1-yt2(:,2); dy12=dy/abs(yyt1(kd) dy12 = 2.8978e-006 【例13-9】用“零阶”近似法求u(t)=e-tU(t)和 h(t)=te-t/2U(t)的卷积。本例演示： （A）连续函数的有限长度采样。 （B）卷积数值计算三个误差（“截尾”误差、“零阶” 近似误差、计算机字长误差）的影响。 （C）卷积“无截尾误差”区间、“非平凡”区间端点的 确定。 （D）离散卷积和连续卷积之间的关系。 （E）指令conv的使用。 （F）绘图分格线的运用。 （1）离散卷积运算前的准备 由于数值计算只能对有限长度序列进行，所以 必须对被卷函数做截尾处理。假如把5%作为函数的截 断阈值，那么u(t)在t=3处截尾，h(t)在t=11处截尾 ，即取t2=3，t4=11。 取采样周期T=0.001。 （2）实施计算 %“零阶”法计算卷积 t2=3; t4=11;T=0.01; tu=0:T:t2;N2=length(tu); th=0:T:t4;N4=length(th); u=exp(-tu);h=th.*exp(-th/2); tx=0:T:(t2+t4);Nx=length(tx); yt3=T*conv(u,h); % 把计算结果与理论值比较 t=tx;yyt1=eval(vectorize(char(yt1); % 例13-8第（1）步的计算结果 dy,kd=max(abs(yyt1(1:N2)-yt3(1:N2); dy13(1)=dy/abs(yyt1(kd); 图13-9 “零阶”近似法卷积在不同区间上的精度图示 dy,kd=max(abs(yyt1(N2+1:N4)-yt3(N2+1:N4); dy13(2)=dy/abs(yyt1(N2+kd); dy,kd=max(abs(yyt1(N4+1:Nx)-yt3(N4+1:Nx); dy13(3)=dy/abs(yyt1(N4+kd); （3）不同区间段的误差及图示 disp(与理论结果的相对误差) disp(blanks(4),0,3段 3,11段 11,14段), disp(dy13) plot(t,yyt1,:b,t,yt3,r) set(gca,Xtick,0,3,11,14), grid 与理论结果的相对误差 0,3段 3,11段 11,14段 0.0068 0.0810 0.6974 。 三、卷积 13.1 数值微分 在工程试验或工程应用中，有时要根据已知的数据 点，求某一点的一阶或高阶导数，这时就要用到数值 微分。与积分相反，数值微分非常困难。积分描述了 一个函数的整体或宏观性质，而微分则描述了一个函 数在一点处的斜率，即函数的微观性质。因此积分对 函数的形状在小范围内的改变不敏感。而微分却很敏 感。一个函数小的变化，容易产生相邻点的斜率的大 的改变。 由于微分这个固有的困难，所以尽可能避免数值 微分，特别是对实验获得的数据进行微分。在这种情 况下，最好用最小二乘曲线拟合这种数据，然后对所 得到的多项式进行微分。或用另一种方法，对该数据 进行三次样条拟合。 数值微分的基本思路是先用逼近或拟合等方法将 已知数据在一定范围内的近似函数求出，再用特定的 方法对此近似函数进行微分。 己知函数某些节点的值，只要将用曲线拟合得到的多 项式微分，再对微分后的多项式求值，即可求出在拟 合范围以内任意一点的任意阶微分。这种方法一般只 用在低阶数值微分。 例13-20example13_20.m 用5阶多项式拟合 函数cos(x)，并利用多项式的求导来求处的一阶和 二阶导数。 x=0:0.3:4; y=cos(x); p=polyfit(x,y,5); 生成拟合多项式 pp=polyder(p) 对拟合多项式求一 阶微分 pp= -0.0330 0.257 0.0027 -1.0257 0.0069 polyvat(pp,pi) 求pi处的一阶导数 ans= 0.0025 ppp=polyder(pp) 求多项式的二阶微 分 pp= -0.1321 0.6200 0.0053 -1.0257 polyval(ppp,pi) 求pi处的二阶导数 ans= 1.0150 例13-21example13_21.m x=0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 y=-.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.30 11.2; % data n=2; % order of fit p=polyfit(x,y,n) % find polynomial coefficients p = -9.8108 20.1293 -0.0317 xi=linspace(0, 1, 100); z=polyval(p , xi); % evaluate polynomial plot(x,y,o,x,y,xi,z,:) xlabel(x), ylabel(y=f(x), title(Second Order Curve Fitting ) 图13-9 二次曲线拟合 在这种情况下，运用多项式微分函数polyder求 得微分。 pd=polyder(p) pd= -19.6217 20.1293 y=-9.8108x2+20.1293x-0.0317的微分是dy/dx=- 19.6217x+20.1293。由于一个多项式的微分是另一个 低一阶的多项式，所以还可以计算并画出该函数的微 分。 图13-11 曲线拟合多项式微分 z=polyval(pd, xi); % evaluate derivative plot(xi, z) xlabel(x), ylabel(dy/dx), title(Derivative of a curve Fit Polynimial) %(微分曲线如图13-11所 示) 在这种情况下，拟合的多项式为二阶，使其微分 为一阶多项式。这样，微分为一条直线，它意味该微 分与x成线性变化。 一、多项式求导法求数值微分 二、用diff函数计算差分求数值微分 给定一些描述某函数的数据，MATLAB提供了一个 计算其非常粗略的微分的函数。这个函数命名为diff ，它计算数组中元素间的差分。 【调用格式】： D=diff(X) 因为微分定义为： 则y=f(x)的微分可近似为： 它是y的有限差分除以x的有限差分。因为diff计 算数组元素间的差分，所以在MATLAB中，可近似求得 函数的微分。继续前一个例子： dy=diff(y)./diff(x); % compute differences and use array division xd=x(1:length(x)-1); % create new x axis since dy is shorter than y plot(xd,dy); title(Approximate Derivative Using DIFF) 图13-12 用diff得到的近似微分 ylabel(dy/dx), xlabel(x) 例13-25example13_25.m 用diff近似求函 数的微分。 x=0:0.3:4; y=cos(x); dy=diff(y)./diff(x) 应用数组近似计算 微分 dy= Columns 1 through 7 -0.1489 -0.4333 -0.6791 -0.8642 - 0.9721 -0.9931 -0.9255 Columns 8 through l3 -0.7752 -0.5556 -0.2864 0.0084 0.3024 0.5694 因函数diff计算的是数组元素间的差分，故所得 输出比原数组少了一个元素。这样，在绘制微分曲线 时，必须舍弃x数组中的一个元素。当舍弃x的第一个 元素时，上述过程给出向后差分近似。而舍弃x的量 后一个元素，则给出向前差分近似。如果比较这两条 曲线，就很容易发现，用有限差分近似微分会导致很 差的结果，特别是在处理被噪音污染了的数据方面。 一、多项式求导法求数值 微分 13.3 符号微积 分 MATLAB的数学工具箱提供了微积分运算的基本函数：包 括微分、积分、求和、极限和泰勒展开等。 1、符号自变量的确定 在MATLAB中进行数学运算时，自变量的选取一般 根据上下文得到。例如，在表达式f=x3和y=sin(at+b) 中，如果对其进行求导而没有指定自变量，则MATLAB 将根据数学约定，分别得到求导式：f=3x2和 y=acos(at+b)，即分别假设这两个表达式中的自变量 为x和t，其他符号a和b都看作常数或参数。 根据数学约定，自变量一般都是小写字母，并且 在拉丁字母表的后面(如x、y和z)。在MATLAB中，可 以用函数findsym找出自变量的符号。 2、findsym函数 函数findsym用来找出自变量的符号。 例13-30example12_30.m 如何确定自变量 。 syms a b c x y f=a*x3+b*x-c; 定义一个符号表达 式 g=sin(x+3*y); 定义另一个符号 表达式 diff(f) 对第一个表达式求导，但并没 有说明对哪一个变量运算 ans= 3*a*x2+b findsym(f,1) 寻找f中的第一个符号 变量 ans= x findsym(g,1) 寻找g中的第一个符号 变量 ans= x findsym(f) 给出表达式f中的所有符号 变量 ans= a,b,c,x 函数findsym的运算规则是在符号表达式中的缺 省变量为靠近小写字母x的优先，在x后面的字母优先 。 F=(1:6).3 建立一个向量 F= 1 8 27 64 125 216 diff(F) 计算元素之间的差分 ans= 7 19 37 61 91 4、多元向量函数的雅可比式 MATLAB中可使用jacobian函数计算变换的雅可比 矩阵J。其调用格式为： J=jacobian(f,v) 例13-42example13_42.m 求函数的 jacobian矩阵 syms r m n x=r*sin(m)*sin(n); y=r*cos(m)*sin(n); z=r*cos(m); J=jacobian(x;y;z,r,m,n) 计算雅可比矩阵 J= sin(m)*sin(n) r*cos(m)*sin(n) r*sin(m)*cos(n) cos(m)*sin(n) -r*sin(m)*sin(n) r*cos(m)*cos(n) cos(m) -r*sin(m) 0 【注意】：jacobiaa函数的第一个输入参数必须 是列向量，第二个输入参数必须是行向量。 因jacobian行列式是非常复杂的三角函数表达式 ，所以用simple命令可对其进行三角变换和简化。 三、符号积分 积分有不定积分、定积分、广义积分和重积分等 。一般说来，积分比微分更难求。符号积分与数值积 分相比，符号积分指令简单，适应性强。但可能占用 的机器时间较长。有时可能得不出“封闭”解，但都给 警告提示和积分原式。 在MATLAB中，函数int(f)命令用来进行符号积分 。该命令试图找到一个符号表达式，使得diff(F)=f ，也就是说int(f)返回f的不定积分，或者是f的微分 的逆运算。与微分类似，int(f,a)使用符号对象a作 为积分自变量。积分命令的调用格式为： （1）int(s) 求符号表达式s关于findsym确定 的默认符号变量进行不定积分； （2）int(s,v) 求符号表达式s关于变量v的不 定积分； （3）int(s,a,b) 求符号表达式s关于默认变 量从a到b的定积分，其中a和b是数值； （4）int(s,v,a,b) 求符号表达式s关于变量v 的从a到b的定积分。 【说明】： （1）正如函数diff一样，积分命令int对符号数 组或矩阵的每一个元素进行运算； （2）a、b分别是积分的上，下限，允许它们取 任何值或符号表达式； （3）当积分限由一具体值变为正负无穷时，定 积分便转化为广义积分，此时也只需将积分限变为无 穷，即可得到相应函数的广义积分值。 例13-45example13_45.m 符号积分运算 syms x u t; A=cos(x*t),sin(x*t);-sin(x*t),cos(x*t); %定义一个符号矩阵 S=sin(u+3*x); %定义一个符号 表达式 int(s) %以默认变量，进行积分 ans= -l/3*cos(u+3*x) int(s,u) %以u为变量进行积分 ans= -cos(u+3*x) int(s,pi/2,pi) %以变量x从pi/2到 pi进行积分 ans= 1/3*cos(u)+1/3*sin(u) int(u,pi/2,pi) %以变量u从pi/2到pi 进行积分 ans= 4*cos(x)3-3*cos(x)- 4*sin(x)*cos(x)2+sin(x) int(s,x,2,sin(t) 以变量x从2到 sin(t)进行积分 ans= - 4/3*cos(u)*cos(sin(t)3+cos(u)*cos(sin(t)+4/ 3*sin(u)*sin(sin(t)*cos(sin(t)2 -1/3*sin(u)*sin(sin(t)+l/3*cos(u+6) int(A,t) 以变量t对矩阵A进行积分 ans= 1/x*sin(x*t) -cos(x*t)/x cos(x*t)/x l/x*sin(x*t) syms x k f=exp(-(k*x)2); int(f,x,-inf,inf) 四、符号求和 当符号变量的和存在时，可用函数symsum进行符 号求和。其调用格式为： （1）symsum(s) 关于默认变量对通项s求不定 和； （2）symsum(s,v) 关于指定变量v对通项s求不 定和； （3）symsum(s,a,b)或symsum(s,v,a,b) 关于 默认变量或指定变量v在a，b之间取值时，对通项s求 和。 例13-47example12_47.m 求下面级数的和 。 syms k n simple(symsum(k2,0,n) %求和 0+12+22+n2 ans= 1/6*n*(n+1)*(2*n+1) symsum(1/k2,1,inf) %求和 1+1/22+1/32+ ans= 1/6*pi2 五、泰勒级数 在MATLAB中，提供有函数taylor，用于对变量进 行泰勒级数的展开。taylor函数的调用命令为： （1）taylor(f) 求f对默认变量的泰勒展开， 六阶小量以余项式给出； （2）taylor(f,n) 求f对默认变量的前(n-1) 项非零泰勒展开多项式； （3）taylor(f,a) 求f对默认变量在点a处的 泰勒展开多项式； （4）taylor(f,x) 用指定变量x代替findsym 所确定的变量。 例13-50example12_50.m 求下面的泰勒级 数 taylor(exp(-x) 求对默认变量的前6个 非零项 ans= 1-x+l/2*x2-1/6*x3+1/24*x4-1/120*x5 taylor(log(x),8,1) 求log(x)的在x=l点 处的泰勒展开式的前面8个非零项 ans= x-1-l/2*(x-1)2+1/3*(x-1)3-l/4*(x- 1)4+1/5*(x-1)5-l/6*(x-1)6+1/7(x-1)7 taylor(xt,3,t) 求xt对指定变量t的泰勒 展开式的前面3个非零项 ans= 1+log(x)*t+l/2*log(x)2*t2 六、符号积分示例 【例13-51】example12_51.m 求。演示：积分 指令对符号函数矩阵的作用。 syms a b x; f=a*x,b*x2;1/x,sin(x); disp(The integral of f is); pretty(int(f) The integral of f is 2 3 1/2 a x 1/3 b x log(x) -cos(x) 【例13-52】example12_52.m 求。演示如何使用 mfun指令获取一组积分值。 （1）求一般积分结果 F1=int(1/log(t),t,0,x) F1 = -Ei(1,-log(x) （2）利用mfun指令求x=0.5, 0.6, 0.7, 0.8, 0.9时的定积分 x=0.5:0.1:0.9 F115=-mfun(Ei,1,-log(x) x = 0.5000 0.6000 0.7000 0.8000 0.9000 F115 = -0.3787 -0.5469 -0.7809 -1.1340 -1.7758 【例13-53】example12_53.m 求积分。注意： 内积分上下限都是函数。 syms x y z F2=int(int(int(x2+y2+z2,z,sqrt(x*y),x2*y),y,sqrt( x),x2),x,1,2) VF2=vpa(F2) %积分结果用32位数字表示 图13-13 交互式近似积分 F2 =1610027357/6563700- 6072064/348075*2(1/2)+14912/4641*2(1/4)+64/225* 2(3/4) VF2 = 224.92153573331143159790710032805 【例13-54】example13_54.m 利用rsums求积分 。（与例12-11结果比较） syms x positive; px=0.5/log(0.5*x); rsums(px) 小 结 本章从工程教学的角度，详细介绍了MATLAB在数 值微积分等方面的应用。通过本章的学习，应该重点 掌握如下内容： （1）熟练掌握一维和多重数值积分的命令，以 及用多项式求导法求数值微分和用diff计算插分法求 数值微分。 （2）熟练掌握符号微积分应用的内容，包括符 号自变量的确定、求函数的极限、对符号表达式(符 号数组和多元向量函数)求导数和积分、符号积分， 符号求和等，同时熟练掌握通过调用函数taylor求函 数的泰勒级数展开式。 一、符号自变量的确 定 1、符号自变量的确定 在MATLAB中进行数学运算时，自变量的选取一般 根据上下文得到。例如，在表达式f=x3和y=sin(at+b) 中，如果对其进行求导而没有指定自变量，则MATLAB 将根据数学约定，分别得到求导式：f=3x2和 y=acos(at+b)，即分别假设这两个表达式中的自变量 为x和t，其他符号a和b都看作常数或参数。 根据数学约定，自变量一般都是小写字母，并且 在拉丁字母表的后面(如x、y和z)。在MATLAB中，可 以用函数findsym找出自变量的符号。 2、findsym函数 函数findsym用来找出自变量的符号。 例13-30example12_30.m 如何确定自变量 。 syms a b c x y f=a*x3+b*x-c; 定义一个符号表达 式 g=sin(x+3*y); 定义另一个符号 表达式 diff(f) 对第一个表达式求导，但并没 有说明对哪一个变量运算 ans= 3*a*x2+b findsym(f,1) 寻找f中的第一个符号 变量 ans= x findsym(g,1) 寻找g中的第一个符号 变量 ans= x findsym(f) 给出表达式f中的所有符号 变量 ans= a,b,c,x 函数findsym的运算规则是在符号表达式中的缺 省变量为靠近小写字母x的优先，在x后面的字母优先 。 二、极限 在实际工作中，极限的求法有很多技巧，因而往 往比较复杂。MATLAB采用函数limit直接计算函数的 极限，其调用格式为： （1）limit(f,x,a) 求符号表达式f当xa时 的极限； （2）limit(f,a) 对系统默认变量且该默认变 量a时表达式f的极限； （3）limt(f) 对系统默认变量且该默认变量 a=0时表达式f的极限； （4）limit(f,x,a,right)或 limit(f,x,a,left) 求x从右侧或从左侧趋近a时 表达式f的极限。 例13-32example12_32.m 求极限 syms x h limit(sin(x+h)-sin(x)/h,h,0) 求 当h0时的极限 ans= cos(x) limit(x-2)/(x2-4),2) 求对默认变量 2时的极限 ans= 1/4 limit(sin(x)/x) 求对默认变量0时 的极限 ans= 1 limit(1/x,x,0) ans= NaN limit(1/x,x,0,right) % x从右侧趋近0时 的极限 ans= inf limit(1/x,x,0, left) %x从左侧趋近0 时的极限 ans= -inf 【说明】： （1）系统缺省情况下limit仍等同于 limit(f,x,0)。 （2）对于极限，不论从左侧趋近a还是从右侧趋 近a，其结果都相同，则f(x)在a的极限存在；如果两 侧的值不同，则说明f(x)在a点的极限不存在。表达 式f(x)在其奇异点处没有极限。 三、导数和微分 由高等数学可知，函数f(x,y,z,)在某一点 (x0,y0,z0,)的增长率就是此函数在该点的导数。对 函数的严格定义为： 可以用前面的limit命令来求各种函数的导数。 而在MATLAB中，求函数的导数或微分由专门的函数 diff来完成。 1、符号表达式的导数和微分 diff函数有四种调用格式，分别为： （1）diff(f) 对默认的自变量求微分，并且 自变量由函数findsym确认； （2）diff(f,a) 对自变量a求微分； （3）diff(f,n) 对默认的自变量求n阶微分， 并且自变量由函数findsym确认； （4）diff(f,a,n) 对变量a求n阶微分。 当a缺省时，自变量自动由findsym确认；n缺省 时，默认n=1。 例13-35example13_35.m 符号表达式的求 导与微分 syms a b c x f; f=a*x4+x3-b*x+c; 定义一个符号 表达式 diff(f) 对默认的变量求微分，等价 于diff(f,x,1) ans= 4*a*x3+3*x2-b diff(f,a) 对变量a求微分 ans= x4 diff(f,2) 对默认的变量求二次微分，等 价于diff(f,x,2) ans= 12*a*x2+6*x diff(f,a,3) 对变量a求三次微分 ans= 0 2、符号数组和符号矩阵的导数和微分 函数diff也可对数组进行运算。如果f是符号数 组或矩阵，则求导对数组或矩阵中的每个元素进行。 例13-37example13_37.m 符号矩阵的微分 syms a b x; F=cos(a*x) sin(a*x);-sin(b*x) cos(a2*x); diff(F) %对F求微分 ans= -sin(a*x)*a cos(a*x)*a -cos(b*x)*b -sin(a2*x)*a2 diff(F,a,2) %对变量a求二阶微分 ans= -cos(a*x)*x2 -sin(a*x)*x2 0 -4*cos(a2*x)*a2*x2- 2*sin(a2*x)*x diff(diff(f,a),x) 求二阶混合导数 ans= 4*x3 3、计算数值向量的数值差分 函数diff也可以计算数值向量或矩阵的数值差分 。对于一个数值向量或矩阵F，diff(F)计算F(2:m,:) -F(1:m-1,:)的数值差分。其调用格式为： Y=diff(F,n,dim) 其中：F是向量或数组；n是差分阶数；dim是指定沿 着数组的那一维进行差分。而n和dim可以省略。 【注意】： 向量或数组的有限差分是按列进行的。 例13-39example13_39.m 数值向量的差分 F=(1:6).3 建立一个向量 F= 1 8 27 64 125 216 diff(F) 计算元素之间的差分 ans= 7 19 37 61 91 4、多元向量函数的雅可比式 MATLAB中可使用jacobian函数计算变换的雅可比 矩阵J。其调用格式为： J=jacobian(f,v) 例13-42example13_42.m 求函数的 jacobian矩阵 syms r m n x=r*sin(m)*sin(n); y=r*cos(m)*sin(n); z=r*cos(m); J=jacobian(x;y;z,r,m,n) 计算雅可比矩阵 J= sin(m)*sin(n) r*cos(m)*sin(n) r*sin(m)*cos(n) cos(m)*sin(n) -r*sin(m)*sin(n) r*cos(m)*cos(n) cos(m) -r*sin(m) 0 【注意】：jacobiaa函数的第一个输入参数必须 是列向量，第二个输入参数必须是行向量。 因jacobian行列式是非常复杂的三角函数表达式 ，所以用simple命令可对其进行三角变换和简化。 三、符号积分 积分有不定积分、定积分、广义积分和重积分等 。一般说来，积分比微分更难求。符号积分与数值积 分相比，符号积分指令简单，适应性强。但可能占用 的机器时间较长。有时可能得不出“封闭”解，但都给 警告提示和积分原式。 在MATLAB中，函数int(f)命令用来进行符号积分 。该命令试图找到一个符号表达式，使得diff(F)=f ，也就是说int(f)返回f的不定积分，或者是f的微分 的逆运算。与微分类似，int(f,a)使用符号对象a作 为积分自变量。积分命令的调用格式为： （1）int(s) 求符号表达式s关于findsym确定 的默认符号变量进行不定积分； （2）int(s,v) 求符号表达式s关于变量v的不 定积分； （3）int(s,a,b) 求符号表达式s关于默认变 量从a到b的定积分，其中a和b是数值； （4）int(s,v,a,b) 求符号表达式s关于变量v 的从a到b的定积分。 【说明】： （1）正如函数diff一样，积分命令int对符号数 组或矩阵的每一个元素进行运算； （2）a、b分别是积分的上，下限，允许它们取 任何值或符号表达式； （3）当积分限由一具体值变为正负无穷时，定 积分便转化为广义积分，此时也只需将积分限变为无 穷，即可得到相应函数的广义积分值。 例13-45example13_45.m 符号积分运算 syms x u t; A=cos(x*t),sin(x*t);-sin(x*t),cos(x*t); %定义一个符号矩阵 S=sin(u+3*x); %定义一个符号 表达式 int(s) %以默认变量，进行积分 ans= -l/3*cos(u+3*x) int(s,u) %以u为变量进行积分 ans= -cos(u+3*x) int(s,pi/2,pi) %以变量x从pi/2到 pi进行积分 ans= 1/3*cos(u)+1/3*sin(u) int(u,pi/2,pi) %以变量u从pi/2到pi 进行积分 ans= 4*cos(x)3-3*cos(x)- 4*sin(x)*cos(x)2+sin(x) int(s,x,2,sin(t) 以变量x从2到 sin(t)进行积分 ans= - 4/3*cos(u)*cos(sin(t)3+cos(u)*cos(sin(t)+4/ 3*sin(u)*sin(sin(t)*cos(sin(t)2 -1/3*sin(u)*sin(sin(t)+l/3*cos(u+6) int(A,t) 以变量t对矩阵A进行积分 ans= 1/x*sin(x*t) -cos(x*t)/x cos(x*t)/x l/x*sin(x*t) syms x k f=exp(-(k*x)2); int(f,x,-inf,inf) 四、符号求和 当符号变量的和存在时，可用函数symsum进行符 号求和。其调用格式为： （1）symsum(s) 关于默认变量对通项s求不定 和； （2）symsum(s,v) 关于指定变量v对通项s求不 定和； （3）symsum(s,a,b)或symsum(s,v,a,b) 关于 默认变量或指定变量v在a，b之间取值时，对通项s求 和。 例13-47example12_47.m 求下面级数的和 。 syms k n simple(symsum(k2,0,n) %求和 0+12+22+n2 ans= 1/6*n*(n+1)*(2*n+1) symsum(1/k2,1,