应用第4版微积分问题的计算机求解

1、第4章 微积分问题的计算机求解4.1极限问题的解析解单变量函数的极限已知函数 f ( x ) ，则极限问题的一般描述为其中， x0可以是一个确定的值，也可以是无穷大。对某些函数来说，还可以定义单边极限（或称左右极限），或前者表示x从左侧趋近于x0点，称为左极限，后者相应地称为右极限。极限问题在MATLAB符号运算工具箱中可以使用1imit函数直接求出，该函数的调用格式为1limit(f,x,a)：求2limit(f,a)：求f(x)的极限值。符号函数f(x)的变量为函数findsym(f)确定的默认自变量a，即xa3limit(f)：求f(x)的极限值。符号函数f(x)的变量为函数findsy

2、m(f)确定的默认自变量。没有指定变量的目标值，系统默认自变量趋近于0，即a0。4limit(f,x,a,right)：求。right表示变量x从右边趋近于a5limit(f,x,a,left)：求。left表示变量x从左边趋近于a例1 求解极限问题。syms x a b;f = x*(1+a/x) x*sin(b/x);L=limit(f,x,inf)例2求syms a m xf = (x (1/m) - a (1/m) / (x-a)limit (f,x,a)例3求syms xf= x * (sqrt(x 2 + 1) - x)limit(f,x,inf,left) %为什么不right？

3、例4求syms a m xf = (sqrt(x) - sqrt(a) + sqrt(x - a) / sqrt(x * x - a * a)limit(f,x,a,right)例5 求解单边极限问题syms x;limit(exp(x 3) - 1)/(1 - cos(sqrt(x - sin(x),x, 0,right)例6 分别求出tanu函数关于 / 2点处的左右极限。syms u;f = tan(u);L1 = limit(f, u,pi/2, left)L2 = limit(f, u,pi/2, right )4. 1. 2多变量函数的极限多元函数的极限可以同样用MATLAB中的l

4、imit函数直接求解。假设有二元函数f (x, y)，若要求出二元函数的极限则可以嵌套使用limit函数。例如，L1 = limit (limit(f,x,x0),y,y0)或L2limit(limit(f,y,y0), x,x0)如果x0或y0不是确定的值，而是另一个变量的函数，例如 x g(y)，则上述的极限求取顺序不能交换。例7 求出二元函数极限值syms x y a;f = exp( -1/(y 2+ x 2)*sin(x) 2 / x 2*(1+1/y 2 ) (x+ a 2 *y 2);L = limit(limit(f,x,1/sqrt(y),y,inf)4. 2 函数导数的解析

5、解4. 2. 1 函数的导数和高阶导数如果函数和自变量都已知，且均为符号变量，则可以用diff函数解出给定函数的各阶导数。diff 函数的调用格式为diff(s)：没有指定变量和导数阶数，则系统按findsym函数指示的默认变量对符号表达式s求一阶导数。diff(s,v)：以v为自变量，对符号表达式s求一阶导数。diff(s,n)：按findsym函数指示的默认变量对符号表达式s求n阶导数，n为正整数。diff(s,v,n)：以v为自变量，对符号表达式s求n阶导数。例1，求syms xy=sqrt(1+exp(x)diff(y) %求1。未指定求导变量和阶数，按默认规则处理例2，求y、 y。s

6、yms xy=x*cos(x)diff(y,x,2) %求2。求y对x的2阶导数diff(y,x,3) %求2。求y对x的3阶导数例3给定函数，试求出 syms x;f=sin(x)/(x 2+4*x+3);f4 = diff(f,x,4)4. 2. 2隐函数的偏导数已知隐函数的数学表达式为，则可以通过隐函数对它们的偏导数求出自变量之间的偏导数。具体可以用下面的公式求出：例4由定义，求，syms a x y zf=x2+y2+z2-a2zx=-diff(f,x)/ diff(f,z) %求z对x的偏导数zy=-diff(f,y)/ diff(f,z) %求z对y的偏导数例5二元函数，求syms

7、 x y;f=(x 2-2*x)*exp(- x 2 - y 2 - x*y);- simple(diff(f,x)/diff(f,y)4 . 2 . 3 参数方程的导数若已知参数方程，则可以由递推公式求出对于简单的一阶和二阶导数，可以直接用下面的公式，例3，求、syms a b t xy1=a*cos(t); y2=b*sin(t);F1 = diff(y2) / diff(y1) %按参数方程求导公式求y对x的导数F2 = (diff(y1) *diff(y2, 2) -diff(y1, 2)*diff(y2) /(diff(y1) 3 %求y对x的2阶导数对于更高阶参数方程的导数，可以编

8、写出如下的递归调用函数来求解。function result=paradiff(y,x,t,n)if mod(n, 1) = 0, error(n should positive integer, please correct)elseif n = 1, result = diff(y,t)/diff(x,t);else, result=diff(paradiff(y,x,t,n -1),t)/diff(x,t);endend例4 已知参数方程，求syms t;y=sin(t)/(t + 1) 3;x=cos(t)/(t + 1) 3;f=paradiff(y,x,t,3); %调用parad

9、iff函数文件n,d=numden(f); %分子和分母分别存放在n与d中F=simple(n)/simple(d)4. 3 积分问题的解析解4.3.1 不定积分MATLAB 符号运算工具箱中提供了一个int函数，可以直接用来求符号函数的不定积分。该函数的调用格式为 F = int ( fun , x ）如果被积函数 fun 中只有一个变量，则调用语句中的 x 可以省略。另外，该函数得出的结果F(x）是积分原函数，实际的不定积分应该是F(x) + C 构成的函数族，其中，C是任意常数。对于可积的函数，MATLAB符号运算工具箱提供的int函数可以用计算机代替繁重的手工推导，立即得出原始问题的解

10、而对于不可积的函数来说，MATLAB也是无能为力的。例1求x=sym(x)f=(3-x2)3int(f)例2求syms alpha tf=exp(alpha*t) int(f)例3求syms x tf=5*x*t/(1+x2)int(f,t)例4 考虑两个不可积问题（1） （2） 解（1）syms x;int(exp(-x2/2)该解中含有erf函数，它的定义为这样似乎可以写出积分的解析表达式。但事实上，这样的结果在工程中是不能用的，必须得出相应的数值解。解（2）syms a x;int(x*sin(a*x4)*exp(x2/2)运行后，将出现如下的错误信息：Warning: Explicit

11、 integral could not be found说明积分不成功。定积分与无穷积分计算在MATLAB中，仍然可以使用int函数来求解定积分或无穷积分问题，该函数的具体调用格式为I = int ( f, x, a, b)其中，x为自变量，(a, b）为定积分的积分区间，求解无穷积分时，允许将a , b 设置成 - Inf 或Inf。当a,b中有一个符号表达式时，函数返回一个符号函数。例1求x=sym(x)int(abs(1-x),1,2)例2求x=sym(x)f=1/(1+x2)int(f,-inf,inf)例3求x=sym(x)f=x3/(x-1)10I=int(f,2,3)double

12、(I) %将符号结果转换为数值例4求syms x tint(4*x/t,t,2,sin(x)多重积分对于多重积分，需要根据实际情况先选择积分顺序，可积的部分作为内积分，然后再处理外积分。每步积分均采用int函数处理。如果交换积分顺序后仍然不能积出解析解，则说明原积分没有解析解，需要采用数值方法求解。例1求二重积分syms t x c1 c2, %c1和c2是积分常数y = int(int(t 2 * exp(- 3 * t * x), t) + c1, x) + c2例2 计算积分。D是由X轴、Y轴和抛物线在第一象限内所围成的区域，如下图所示。解：先把二重积分化成二次积分：采用先对变量y积分的

13、二次积分：syms x y,int(int(3 * x 2 * y 2, y, 0, 1 - x 2) , x, 0, 1)4.4函数的级数展开与级数求和问题求解 Taylor幂级数展开如果在x = a点附近进行Taylor幂级数展开，则得出其中各个系数认可以如下求出。Taylor幂级数展开可以用符号运算工具箱的taylor函数给出，其调用格式为taylor (f, x, k) %按 x = 0 进行taylor幕级数展开taylor (f, x, k, a) %按 x = a 进行taylor幕级数展开其中， f为函数的符号表达式，x为自变量，若函数只有一个自变量，则x可以省略。k为需要展开

14、的项数，默认值为 6 项。例1求的5阶泰勒级数展开式。x=sym(x)f1=sqrt(1-2*x+x3)-(1-3*x+x2)(1/3)taylor(f1,x,5)例2将在x=1处按5次多项式展开。x=sym(x)f1=(1+x+x2)/ (1-x+x2)taylor(f1,6,1)%展开到x-1的5次幂时应选择n=6 Fourier级数展开给定函数f (x)，其中，x- L,L，且周期为 T 2L，可以人为地对该函数在其他区间上进行周期延拓，使得f (x) = f(kT + x)，k为任意整数。这样可以把函数展开成无穷三角函数和形式，即把函数展开为下面级数的形式：其中，该级数称为Fourie

15、r级数，an，bn称为Fourier系数。若x (c, d)，则可以计算出 L = (d - c) / 2。这时可以引入新变量，使得，则可以将映射成- L,L区间上的函数。这样可以对其进Fourier级数展开。然后再将转换成x的函数即可。MATLAB没有直接提供求解Fourier系数与级数的现成函数。可以由上述公式很容易地编写求解Fourier级数的函数文件。function A,B,F=fseries(f,x,n,a,b)if nargin=3, a=-pi; b=pi; end % nargin是实际函数输入参数个数L=(b-a)/2; if a+b, f=subs(f,x,x+L+a);

16、 end%如果x区间对于y轴非对称，就将x置换为x+L+aA=int(f,x,-L,L)/L; B=; F=A/2; %求a0/2 %以下循环是求an和bnfor i=1:nan=int(f*cos(i*pi*x/L),x,-L,L)/L; bn=int(f*sin(i*pi*x/L),x,-L,L)/L; A=A, an; B=B,bn; F=F+an*cos(i*pi*x/L)+bn*sin(i*pi*x/L);endif a+b, F=subs(F,x,x-L-a); end %如果x区间对于y轴非对称，再将x置换回x-L-asubs(S, x, x) 表示在符号表达式S中用x代替x。其

17、中x是符号变量或表示变量名的字符串，x是符号或数值变量或表达式。该函数的调用格式为A,B,F=fseries(f,x,p,a,b)，其中，f为给定函数，x为自变量， p为展开项数，a , b为x的区间，默认值为, - 。A , B为Fourier 系数，F为展开式。例1 求给定函数的Fourier级数前12项的展开。syms x; f=x*(x-pi)*(x-2*pi);A,B,F=fseries(f,x,12,0,2*pi);F下面的语句给出了12项 Fourier级数展开对原函数的拟合情况ezplot(f, 0, 2 * pi) %画原函数曲线hold onezplot(F, 0, 2 *

18、 pi) %画级数展开的曲线函数的拟合效果是很理想的。如果想比较更大区间内的拟合效果，比如x ( - , 3)，则可以用下面语句ezplot(f, - pi, 3 * pi) %画原函数曲线hold onezplot(F, - pi, 3 * pi) %画级数展开的曲线这时的拟合效果在（0，2）区间内仍然很理想，但在其他区间内，差别非常大。这是因为Fouricr 级数是定义在周期延拓基础上的，所以在其他区间和原函数完全不同。例2 考虑( - , )区间的方波信号，假设 x 0时，y = 1，否则，y - l。试对该方波信号进行Fourier级数拟合，并观察用多少项能有较好的拟合效果。解：给定的

19、函数可以由表示，由下面语句可以生成 x 轴数据点，求出理论的方波数值，syms x; f=abs(x)/x; % 定义方波信号xx=-pi:pi/200:pi; xx=xx(xx=0); xx=sort(xx,-eps,eps); % 剔除零点yy=subs(f,x,xx); plot(xx,yy), hold on % 绘制出理论值并保持坐标系这里，sort函数的作用是将元素按顺序排列下面是通过不同阶次的Fourier 级数展开去拟合原来的方波函数。for n=1:20; a,b,f1=fseries(f,x,n); y1=subs(f1,x,xx); plot(xx,y1);end最后的“

20、xx,y1”是n=20时的结果，现用有颜色的曲线凸出出来：plot(xx,y1,r);从得出的结果看，当阶次等于10左右就能得出较好的拟合，再增加阶次也不会有显著的改善。如果比较区间扩展到( -2 , 2 )，可以用下面语句来比较拟合情况。取n = 14的情况。syms x; f=abs(x)/x; % 定义方波信号xx=-pi:pi/200:pi; xx=xx(xx=0); xx=sort(xx,-eps,eps); % 剔除零点yy=subs(f,x,xx); plot(xx,yy), hold on % 绘制出理论值并保持坐标系a,b,f1=fseries(f,x,14); f1syms

21、 x; f=abs(x)/x; % 定义方波信号xx=-2*pi:pi/200:2*pi; xx=xx(xx=0); xx=sort(xx,-eps,eps); % 剔除零点yy=subs(f,x,xx); figure, plot(xx,yy), hold on % 绘制出理论值并保持坐标系y1=subs(f1,x,xx); plot(xx,y1)可以看出，在指定区间以外的周期延拓区间， Fourier 级数与原函数无关。4.4.3级数求和的计算级数的和可以表达为下面的形式：其中，fk为级数的通项， k为级数自变量，k0和kn为级数求和的起始项与终止项。MTATLAB符号运算工具箱中提供了s

22、ymsum函数可以用于已知通项的有穷或无穷级数的和。其调用格式为：symsum(fk,k,k0,kn)可以将起始项或终止项设置成无穷量inf 。如果fk中只含有一个变量，则在调用symsum函数时可以省略k量。例 求下列级数之和。1n=sym(n)s1= symsum(1/n2,n,1,inf)2n=sym(n)s2= symsum(-1)(n+1)/n,1,inf) %未指定求和变量，默认为n3syms n xS3= symsum(n*xn,n,1,inf) %此处的求和变量n不能省略4n=sym(n)S4= symsum(n 2,1,100) %计算有限级数的和4. 5曲线积分与曲面积分的

23、计算MATLAB并未直接提供曲线积分和曲面积分的现成函数。曲线积分和曲面积分可以转换成一般积分问题。这样就可以利用 MATLAB 语言的符号运算工具箱来求解曲线积分和曲面积分的解析解。4.5.1 曲线积分及MATLAB求解一、第一类曲线积分第一类曲线积分问题起源于对不均匀分布的空间曲线总质量的求取。假设在空间曲线L上的密度函数为f (x, y, z) ，则其总质量可以由下面的式子直接求出该式即为第一类曲线积分。其中，s为曲线上某点的弧长，所以这类曲线积分又称为对弧长的曲线积分。若 x , y , z 均由参数方程x = x (t)，y = y (t)，z = z (t)给出，则可以将这些量直接

24、代入f (x, y, z)函数，而弧长可以表示成，简记作则可以将这类曲线积分变换成对参数t的普通定积分问题若被积函数为二元函数f (x, y)，也可以用相应的转换方法将其转换成普通积分问题。故可以用 MATLAB求出第一类曲线积分的值。例1 求，其中L为螺线，(, )。syms t;symsapositive; x=a*cos(t); y=a*sin(t); z=a*t;I=int(z2/(x2+y2)*sqrt(diff(x,t)2+diff(y,t)2+diff(z,t)2),t,0,2*pi)（附：定义特殊类型的符号变量语法： syms a 类型; 或者 a=sym(a,类型); 两种语

25、句效果相同, 注意的是他们的区别在于sym中的类型一定要加单引号！ 这里的类型可以是“real”， “unreal”，“positive”。 这样定义的好处是， 如果定义a为positive类型， 那么在之后的计算中， a都只会被赋予正的值。 例如， 如果要解一个方程： a*a=1,那么给出的解就只有a=1， 而自动将 a 0。解：按上面的顺序将积分曲面S的4个积分平面记为S1、S2、S3、S4，原积分可以简写为下面4个曲面积分之和在S1、S2、S3平面，由于被积函数的值为0，故这些积分也为0，所以只需求S4的曲面积分。S4平面的数学表示为：，则，syms x y;syms a positiv

26、e;z=a-x-y;I=int(int(x*y*z*sqrt(1+diff(z,x)2+diff(z,y)2),y,0,a-x),x,0,a)结果应为：若曲面由参数方程，给出，则曲面积分可以由下面的公式求出式中，是曲面S向u-平面上投影，例6 求曲面积分，其中S为螺旋曲面，的(,)部分。syms u v;syms a positive;x=u*cos(v); y=u*sin(v); z=v;E=simple(diff(x,u)2+diff(y,u)2+diff(z,u)2);F=diff(x,u)*diff(x,v)+diff(y,u)*diff(y,v)+diff(z,u)*diff(z,v

27、);G=simple(diff(x,v)2+diff(y,v)2+diff(z,v)2);I=int(int(x2*y+z*y2)*sqrt(E*G-F2),u,0,a),v,0,2*pi)二、第二类曲面积分第二类曲面积分又称为对坐标的曲面积分。其数学定义为其中，正向曲面S+由给出，被积函数为行向量，而为列向量。这类曲面积分问题可以转换成第一类曲面积分其中，z由代替，且，而这样，各个cos（）的分母上的和积分式中的相抵消，整个曲面积分可以写成若曲面由参数方程，给出，则可以求出，其中，。经过整理这时整个曲面积分可以简化成式中，是曲面S向u-平面上投影例7 求曲面积分，其中S是椭球面的上半部，且沿

28、椭球面的上面。解：可以引入参数方程，,且(,)。这样，原始曲面积分问题可以转换为一般双重积分问题，其中，syms u v; syms a b c positive;x = a*sin(u)*cos(v); y=b*sin(u)*sin(v); z=c*cos(u);C=diff(x,u)*diff(y,v) - diff(y,u)*diff(x,v); R=x*y+zI=int(int(R*C,u,0,pi/2) ,v,0,2*pi)结果应为：4.6 数值微分前面介绍了已知原函数，可以通过diff函数求取各阶导数解析解的方法。应该指出的是，这种解析解方法的前提是原型函数为已知的。如果函数表达式

29、未知，只有实验数据，这样的求导问题就不能用前面的方法获得问题的解析解。要求解这样的问题，需要用数值算法求解。由于在MATLAB中没有现成的数值微分函数，数值微分求解是借助差分运算实现的。 数值微分算法假设y是自变量x的函数，。如果自变量由x变化到x + x，相应地函数由变化到，则在x点处y对x的导数定义为一般来说，函数的导数仍然是函数。设f(x)的导数。高等数学关心的是g(x)的形式及性质，而数值分析则关心怎样计算g(x)在一组离散点X=(x1, x2, ,xn)的近似值G(g1, g2, gn)，以及所计算的近似值有多大误差。如果x等间隔地测量，得到X0, x1, xi, xn，测出了一组y

30、数据为y0, y1, yi, yn，x的间隔为x1 x1x0,x2 =x2x1 , xi-1=xixi-1,xi =xixi-1,xn=xnxn-1且x1 x2 = =xi-1 =xi =xn=h相应地，y的间隔为y1 y2y1,y2 =y3 -y2 , yi-1 =yiyi-1,yi =yiyi-1,yn=ynyn-1那么函数f(x)在xi点的导数可以通过下面的任一公式来实现引进记号 函数在x点处以h(h0)为步长的向前差分 函数在x点处以h(h0)为步长的向后差分函数在x点处以h(h0)为步长的中心差分 函数在x点处以h(h0)为步长的向前差分商 函数在x点处以h(h0)为步长的向后差分商

31、函数在x点处以h(h0)为步长的中心差分商在h0的极限过程，函数在x的微分df接近在该点的任意差分、或，而函数在x的导数f接近在该点的任意差分商，在MATLAB中，没有直接提供求数值导数的函数，只有计算向前差分的函数diff，其调用格式为：DX=diff(X)：计算向量X的向前差分，DX(i)=X(i+1)-X(i)，i=1,2,n-1。DX=diff(X,n)：计算X的n阶向前差分。例如，diff(X,2)=diff(diff(X)。DX=diff(A,n,dim)：计算矩阵A的n阶差分，dim=1时(缺省状态)，按列计算差分；dim=2，按行计算差分。例1生成以向量V=1,2,3,4,5,

32、6为基础的范得蒙矩阵，按列进行差分运算。V=vander(1:6) %生成范得蒙矩阵VDV=diff(V) %计算V的一阶差分从上例看到，每差分一阶，矩阵的行或列的长度少1。例2设，用符号方法和数值方法求函数f(x)的数值导数，并在同一个坐标系中做出f(x)的图像。方法1 用subs函数实现符号变量向数值变量的转变syms x1f1 = sqrt(x1.3+2*x1.2-x1+12)+(x1+5).(1/6)+5*x1+2 %符号原函数g1 = diff(f1,x1, 1) %符号导函数x=-3:0.01:3; %生成采样点x2=x,3.01; %生成采样点f=subs(f1,x1,x2);

33、%实现符号变量向数值变量的转变g=subs(g1,x1,x); %由符号导数转变的精确导数值。为什么f的自变量x2的长度比g的自变量x的长度多一个?dx=diff(f)/0.01; %直接对f(x)求数值导数。plot(x,dx,b.,x,g,r-); legend(数值导数, 精确导数) %作图。为什么这里dx的自变量用x而不用x2?方法2 用inline函数实现符号变量向数值变量的转变第1步 用符号方法求函数f(x)的导数：syms xf = sqrt(x.3+2*x.2-x+12)+(x+5).(1/6)+5*x+2 %原函数g = diff(f,x, 1) %导函数得到用符号方法求函数

34、f(x)的导数表达式为：g =1/2/(x3+2*x2-x+12)(1/2)*(3*x2+4*x-1)+1/6/(x+5)(5/6)+5第2步 求函数f(x)的数值导数，并将两种方法得到的结果作图比较：f=inline(sqrt(x.3+2*x.2-x+12)+(x+5).(1/6)+5*x+2); %将符号表达式f变成内联函数，供数值运算用。Inline的用法见下面说明g=inline(3*x.2+4*x-1)./sqrt(x.3+2*x.2-x+12)/2+1/6./(x+5).(5/6)+5); %将符号表达式g变成内联函数，为作图做准备x=-3:0.01:3; %生成采样点dx=dif

35、f(f(x,3.01)/0.01; %直接对f(x)求数值导数，为什么加一个3.01?gx=g(x); %求函数f的导函数g在采样点的导数(精确值)plot(x,dx,b.,x,gx,r-); legend(数值导数, 精确导数) %作图(附：内联函数inline是MATLAB提供的一个对象(Object)。它的性状表现和函数文件一样，但内联函数的创建比较容易。inline(CE)CE是字符串，CE为不包含赋值符号“=”的表达式。上式把串表达式转化为输入宗量自动生成的内联函数。上述调用格式将自动地对CE进行辫识，把CE中由字母数字组成的连续字符认做变量。即除“预定义变量名（如 i , j ,

36、pi ) ”和“常用函数名（如 sin )”以外的由字母数字组成的连续字符将被认做变量。但注意：若连续字符后紧接“左圆括号”，那么将不被当作输入宗量。如 x ( 1 ) ，就不会认做输入宗量处理。inline(CE,arg1,arg2,)上述调用格式把串表达式转化为arg1,arg2等指定输入宗量的内联函数；这种调用格式是创建内联函数的最稳妥、可靠途径。输入宗量字符可表达得更自如。)上述所介绍的基于向前差分和向后差分的数值微分算法在求高阶微分时的精度一般都是很低的。这时应该采用中心差分的算法。下面介绍两种中心差分的算法。中心差分计算函数的一阶导数用上面介绍公式，即：如果再在xi-1,xi,xi

37、+1这些点中间插入一点变成, , , , , 。则，xi =xi+1xi由于i是通项，取i和取2i不应该影响结果，所以用表示由Taylor级数将上式展开其中， 由此，表示，含的项非常小的意思。可以看出，这种算法的精度为。用这种中心差分算法的高阶微分公式为另一种具有精度级的中心差分算法为针对上述具有精度级的中心差分算法，可以编写下述函数文件： function dy,dx=diff_ctr(y, Dt, n) yx1=y 0 0 0 0 0; yx2=0 y 0 0 0 0; yx3=0 0 y 0 0 0; yx4=0 0 0 y 0 0; yx5=0 0 0 0 y 0; yx6=0 0 0

38、 0 0 y; switch n case 1 dy = (-diff(yx1)+7*diff(yx2)+7*diff(yx3)-diff(yx4)/(12*Dt); L0=3; case 2%-(y2-y1)+7(y1-y0)+7(y0-y-1)-(y-1-y-2) dy=(-diff(yx1)+15*diff(yx2)-15*diff(yx3)+diff(yx4)/(12*Dt2);L0=3; case 3 dy=(-diff(yx1)+7*diff(yx2)-6*diff(yx3)-6*diff(yx4)+. 7*diff(yx5)-diff(yx6)/(8*Dt3); L0=5; ca

39、se 4 dy = (-diff(yx1)+11*diff(yx2)-28*diff(yx3)+28*diff(yx4)-. 11*diff(yx5)+diff(yx6)/(6*Dt4);L0=5; end dy=dy(L0+1:end-L0); dx=(1:length(dy)+L0-2-(n2)*Dt; %冒号两边先运算这个M函数的调用格式dy, dx = diff_ctr(y, t, n)其中，y为给定的等间距的实测数据构成的向量，t为自变量的间距；n为所需的导数阶次。向量dy为得出的导数向量，而dx为相应的自变量向量。注意这两个向量的长度比y短。例3 给定函数，用数值微分求原函数的14

40、阶导数，并和解析解比较。h=0.05; x=0:h:pi;%生成一个横坐标点组成的向量xsyms x1; y=sin(x1)/(x12+4*x1+3);yy1=diff(y); f1=subs(yy1,x1,x); % 求各阶导数的解析解与对照数据yy2=diff(yy1); f2=subs(yy2,x1,x); yy3=diff(yy2); f3=subs(yy3,x1,x);yy4=diff(yy3); f4=subs(yy4,x1,x);下面的语句是通过数值的方法获得已知数据点处的各阶导数，同时绘制出曲线，将数值导数和由解析解计算出来的导数在相同的坐标系下绘制出来，供比较y=sin(x)

41、./(x.2+4*x+3); % 生成已知数据点y1,dx1=diff_ctr(y,h,1); subplot(221),plot(x,f1,dx1,y1,:);y2,dx2=diff_ctr(y,h,2); subplot(222),plot(x,f2,dx2,y2,:)y3,dx3=diff_ctr(y,h,3); subplot(223),plot(x,f3,dx3,y3,:);y4,dx4=diff_ctr(y,h,4); subplot(224),plot(x,f4,dx4,y4,:)可以看出，由中心差分算法获得的导数是很精确的，其误差从图上是看不出来的。4.7 数值积分一、数值积分

42、基本原理在数学上，定积分定义为其中，a, b为积分区间，在实际计算中，不可能取，。因此，数值积分的目的是通过在有限个采样点上计算的值来逼近在区间a, b内的定积分。定义：设，如果：具有下面的性质：则表达式称为数值积分或面积公式。称为积分的截断误差，称为面积节点，称为权。根据应用的需要，节点的选择有很多种方法。常用的梯形公式、辛普生公式、以及布尔公式，都选择等距的节点；而高斯一勒让德公式中的节点选择某些勒让德多项式的O点。对于任何应用，都需要了解一些关于数值解的精度问题。定义：如果一个正整数n使得对所有次数的多项式满足，而对某些次数为n+1的多项式有，那么这个正整数n称为面积公式的精度。通过研究

43、当为多项式时的情形可以预测的形式，考虑任意i次多项式：若，则对所有x，有，且。因此截断误差的一般形式为：其中，K是一合理选择的常数，n为精度。这一结果的证明可在数值积分的高级教程中找到。二、闭型牛顿一柯蒂斯公式经过分析发现，如果存在一个唯一通过M+1个等距点、且次数小于等于M的多项式。当用该多项式来近似区间a, b内的时，的积分就近似等于的积分，这一结果的公式称为牛顿一柯蒂斯公式。如果使采样点和，这时的牛顿一柯蒂斯公式称为闭型牛顿一柯蒂斯公式。下图为多项式次数为M =1,2,3,4时利用闭型牛顿一柯蒂斯公式的情况。( a ) x0, x1 = 0.0, 0.5 上用的梯形公式求积分( b )

44、x0, x2 = 0.0, 1.0 上用的辛普生公式求积分( c ) x0, x3 = 0.0, 1.5 上用的辛普生3/8公式求积分( d ) x0, x4 = 0.0, 2.0 上用的布尔公式求积分下面给出了上述几个闭型牛顿一柯蒂斯面积公式。设为等距节点，k = 0, 1， ， M，相邻节点距离为h，且。前4个闭型牛顿一柯蒂斯面积公式为：（梯形公式）（辛普生公式）（辛普生3/8公式）（布尔公式）牛顿一柯蒂斯公式精度：设充分可微，使牛顿一柯蒂斯面积公式的包含更高阶导数项。梯形公式的精度为n =1 。若，则：其中，表示函数自身和它的前n阶导数在区间连续的所有函数的集合。辛普生公式的精度为n =

45、3。若，则：辛普生3/8公式的精度为n = 3。若，则：布尔公式的精度为n = 5。若，则：n阶牛顿柯特斯公式为其中Ci称为柯特斯系数当n=1和2时，分别对应梯形公式和辛普生公式当n=4时，特别称为牛顿柯特斯公式求内曲线下的面积的直观方法，是用区间上的一系列梯形的面积来近似。这就是组合梯形的概念。组合梯形公式：设区间被等距节点, k = 0, 1， ， M，分为宽度为的 M 个子区间。那么M个子区间的组合梯形公式可以表示为下面三种等价方式中的任何一种： (a) (b) (c)这是区间内积分的一种逼近，写为：组合梯形方法组合梯形公式误差：即：组合梯形程序：function s = trapr1

46、(f, a, b, M) %Input - f is the integrand input as a sring f% - a and b are upper and lower limits of integration% - M is the number of subintervals%Output - s is the trapezoidal rules sumh =(b-a)/M; s = 0; for k=1:(M-1);x = a + h * k; s = s + feval(f, x); %求a、b点以外等间距点的f值ends = h* (feval(f, a) + feva

47、l(f, b) / 2 + h*s; %按组合梯形公式计算(feval函数功能：y1,y2, = feval (FN,arg1, arg2, ) 表示用参量arg1, arg2等执行FN函数指定的计算。FN只能是函数名)例用组合梯形方法求定积分，取15个子区间(1) 组合梯形程序文件 traprl.m (见上面)(2) 被积函数文件fx.m。function f=fx(x)f=x.*sin(x)./(1+cos(x).*cos(x);(2) 调用函数traprl求定积分。I=trapr1(fx,0,pi,15)组合辛普生公式：设区间由, k = 0, 1， ， 2M，分为2M个宽度为的等距子区

48、间。那么2M个子区间的组合辛普生公式可表示为下面三种等价方式之一： (a) (b) (c)这是区间内积分的一种逼近，写为：a辛普生方法 b 组合辛普方法组合辛普生公式误差：即：function s=simpr1(f, a, b, M) %Input -f is the integrand inpnt as a string f% - a and b are upper and lower limits of integration% - M is the number of subintervals% Output - s is the simpson rule sum h=(b-a)/(2*

49、M);s1 = 0; s2 = 0; for k=1:M;x=a+h*(2*k-1); s1=s1+feval (f, x); end for k=1:(M-1);x=a+h*2*k;s2=s2+feval(f,x); ends=h*(feval(f,a)+feval(f,b)+4*s1+2*s2)/3;%按组合辛普生公式计算例用组合辛普生方法求定积分，取15个子区间(1)组合辛普生程序文件 simpr1.m (见上面)(2) 被积函数文件fx.m。function f=fx(x)f=x.*sin(x)./(1+cos(x).*cos(x);(2) 调用函数simpr1求定积分。I=simpr1(fx,0,pi,15)如果将区间分成的子区间