Matlab课程设计--数值微积分

1、；.数值微积分引言 在微分中，函数的导数是用极限来定义的，如果一个函数是以数值给出的离散形式，那么它的导数就无法用极限运算方法求得，当然也就更无法用求道方法去计算函数在某点处的导数。 一般来说，函数的导数依然是一个函数。设函数f(x)的导数f(x)=g(x),高等数学关心的是g（x）的形式和性质，而数值分析关心的问题是怎样的计算g（x）在一串离散点X=（x1,x2,xn）的近似值G=（g1,g2,.gn）以及所计算的近似值有多大的误差。1.    数值差分与差商任意函数f(x)在x点的导数是通过极限定义的： f(x)=lim f(x+h)-f(x)/h f(x)=

2、lim f(x)-f(x-h)/h f(x)=lim f(x+h/2)-f(x-h/2)/h 上述式子中，均假设h0，如果去掉上述等式右端的h0的极限过程，并引进记号：f(x)=f(x+h)-f(x)f(x)=f(x)-f(x-h)f(x)=f(x+h/2)-f(x-h/2) 称f(x) f(x)、f(x)分别为函数在x点的以h（h>0）为步长的向前差分、向后差分、向中差分。当步长h充分的小时，有  f(x）f(x)/h f(x）f(x)/h f(x）f(x)/h 和差分一样，称f(x)/h, f(x)/h,、f(x)/h分别为函数在x点一h（h>0）

3、为步长的向前差商、向后差商、向中差商。当步长h充分的小时，函数f在x点微分接近于函数在该点的任意种差分，而f在点x的导数接近于函数在该点的任意种差商。 2.    数值微分的实现有两种方式计算任意函数f（x）在给定点x点数值导数。第一种方式是用多项式或样条函数g（x）对f(x)进行逼近，然后用逼近函数g(x)在点x的导数作为f(x)在点x的导数。第二种方式是用f（x）在点x处的某种差商作为其导数。在MATLAB中，没有直接提供数值导数的函数，只有计算向前差分的函数diff，其调用格式为： Diff(X),输入参数x可为矢量或矩阵。若X为矢量，则返回X(2

4、)-X(1),X(3)-X(2),X(n)-X(n-1);若X为矩阵，则返回矩阵每行的差分，即X(2:m，：)-X(1:m-1,:).通常，diff(X)返回沿X的第一个成对维的差值。 Y=diff(X,n):递归调用diff函数n次，生成n次差分。 Y=diff(X,n,dim): 返回X中第dim维的n次差分或导数值 设计1如下： X=3 7 5;0 4 2; Diff(x) ans= -3 -3 -3 设计2如下： x=3 7 5;0 4 2; Diff=(x,2) ans=0       

5、0;         0   设计3如下：x=3 7 5；0 4 2；diff=(x,2,1)ans= empty matrix: 0-by-3  设： 用不同的方法求函数f(x)的数值导数，并在同一个坐标系中画出f(x)的图形 程序如下： f=inline(sqrt(x.3+2*x.2-x+12)+(x+5).(1/6)+5*x+2); g=inline(3*x.2+4*x-1)./sqrt(x.3+2*x.2-x+12)/2+1/6./(x+5).(5

6、/6)+5); x=-3:0.01:3; p=polyfit(x,f(x),5); 用5次多项式p拟合f(x) dp=polyder(p); 对拟合多项式p求导数dpdpx=polyval(dp,x); 求dp在假设点的函数值dx=diff(f(x,3.01)/0.01; 直接对f(x)求数值导数 gx=g(x); 求函数f的导函数g在假设点的导数plot(x,dpx,x,dx,.,x,gx,-); 作图  程序运行后得到图6.3所示的图形。结果表明，用3种方法求得的数值导数比较近似。        &#

7、160;  MATLAB语言提供了计算给定向量差分的函数diff(),其调用的方法是dy=diff(y).假设向量y是yi,i=1,2,3,n构成，则经过diff()函数处理后将得出一个新的向量：yi+1-yi,i=1,2,3,n,显然新得出的向量长度比原来向量的长度y要小1.如：  >>v=vander(1:6) v= 1 1 1 1 1 1 32 16 8 4 2 1 243 81 27 9 3 1 1024 256 64 16 4 1 3125 625 125 25 5 1 7776 1296 216 36 6 1  >>d

8、iff(v) ans= 31 15 7 3 1 0 211 65 19 5 1 0 781 175 37 7 1 0 2101 369 61 9 1 0 4651 671 91 11 1 0 可见，diff（）函数对矩阵的每一列都进行了差分运算，故而结果矩阵的列数是不变的，只是有行数减1.  同时，在MATLAB中，还提供了gradient()函数的调运，这函数可以直接用来求一个矩阵的二维差分，该函数的调用格式为：  dx,dy=gradient(A) >> v=vander(1,6),dx,dy=gradieng(v) v= 1 1 1 1 1 1

9、32 16 8 4 2 1 243 81 27 9 3 1 1024 256 64 16 4 1 3125 625 125 25 5 1 7776 1296 216 36 6 1 dx= 1.0e=003* 0 0 0 0 0 0 -0.0160 -0.0120 -0.0060 -0.0030 -0.0015 -0.0010 -0.1620 -0.1080 -0.0360 -0.0120 -0.0040 -0.0020 -0.7600 -0.4800 -0.1200 -0.0300 -0.0075 -0.0030 -2.5000 -1.5000 -0.3000 -0.0600 -0.0120

10、-0.0040 -6.4800 -3.7800 -0.6300 -0.1050 -0.0175 -0.0050 dy= 31 15 7 3 1 0 121 40 13 4 1 0 496 120 28 6 1 0 1441 272 49 8 1 0 3376 520 76 10 1 0 4651 671 91 11 1 0   6.2.2 数值积分 在科学实验和生产中，经常要求函数f(x)在区间a,b上的定积分： 在高等数学中，计算积分依靠微积分基本定理，只要找到被积函数f(x)的原函数F（x），则可以用牛顿莱布尼茨公式： 来求出定积分。但是，在有些情况

11、下，应用牛顿莱布尼茨公式往往有困难，例如，当被积函数的原函数无法用初等函数表示或被积函数为仅知离散点处函数值的离散函数时。 1.    数值积分的基本原理数值积分研究定积分的数值求解方法。设 I1= 高等数学中知道，当f(x)-p(x) <时，I1-I2<(b-a).这说明。当充分的小时，可用I2近似的代替I1。所以，求任意的函数f（x）在a,b上的定积分时，要是难以使用解析的方法求出f(x)的原函数，则可以寻找一个在a,b上与f(x)逼近，但形式上却简单易求的函数p(x),用p(X)在a,b上的积分值近似的代表f(X)在a,b上的积分值。一

12、般选择被积分函数的选择一次多项式时，称为梯形公式。选择二次多项式时，为辛谱森（simpson）公式。在MATLAB中，我们通常可以使用梯形求积，Smpson求积，Lobotta求积以及二重积分和三重积分运算，另外还有Gass求积和Romberg求积算法。这些方法的基本思想就是把整个的积分空间分割成为若干个子空间，每个空间上的函数积分可以求取，因而在整个的空间函数上积分可以得到，MATLAB基于这种思想采用自适应变步长的方法给出quad()函数来求积分。现在分别举例说明：1.    梯形求积 用trapz函数进行梯形数值积分。 设：   的精确值为2，下面

13、用trapz函数在均匀间隔的网络上求取该积分的数值近似。 X=0:pi/100:pi; Y=sin(X); Z=trapz(X,Y) Z= 1.9998 对于非均匀的间隔情况，例如： X=sort(rand(1.101)*pi); Y=sin(X); Z=trapz(X,Y) Z=1.9981  2.Simpson求积用quad和quad8函数进行自适应递归Simpson求积。要求积分具有下面的形式： Q=以quad函数为例，quad和quad8函数具有下面的语法形式。         

14、;     q=quad(fun,a,b);求函数fun从a到b的积分近似，误差为le-6.fun为M文件函数或匿名函数的句柄，函数x=fun(x)应该接受矢量变量x并返回结果矢量y。              q=quad(fun,a,b):使用绝对误差容限tol代替默认的容限1.0e-6。Tol的值更大时，需要的计算次数少，计算快。       

15、0;      q=quad(fun,a,b,tol,trace):在递归过程中使用非0的trace参数。              q.fcnt=quad(): 返回函数计算次数 用quad函数计算下面积分的数值积分 首先编写函数的M文件myfun.m。 Function y=myfun(x) Y=1./(x.3-2*x-5); 然后在命令窗口键入 Q=quadruple（myfun.0,2） Q= -0.4605

16、 也可以直接使用匿名函数形式，如： F=（x）1.9(x.3-2*x-5)； Q=quad(F,0,2); 用以上两种方法求：  先建立一个函数文件ex.m: Function ex=ex(x) ex=exp(-x.2);然后在MATLAB命令窗口，输入命令： Format long I=quad(ex,0,1) I= 0.74682418072642 I=quadl(ex,0,1) I= 0.74682413398845 也可以不建立关于被积函数文件，而使用语句函数求解，命令如下： q=inline(exp(-x.2); I=quadl(g,0,1) I= 0.7468241339

17、8845 Format short 读者可以用解析方法计算下本例，并可以比较一下。同样。对于下列函数 Humps(x)=   可编写名为humps.m的M文件，内容是 Functions y=humps(x) y=1./(x-.3).2+0.01)+1./(x-.9).2+0.04)-6； 该函数的图形如下： 生成的方法如下： X:-1:0.01:2;Polt(x,humps(x)                现在

18、需要求humps从0到1积分，可使用下面的命令： >>q=quad(humps,0,1) q= 29.8583 注意quad的一般调用格式为 y,n=quad(F,a,b,tol)其中F为描述被积分的字符变量，一般为一个F.m函数文件名，该函数的一般格式为y=F(x);a,b为上下限；tol为变步长积分用的误差限。 MATLAB中常见的一元函数数值积分指令如下 表2  quad自适应辛谱森积分Quad8牛顿8段积分quadl自适应Lobtto积分trapz梯形法求积分sum等宽矩形法求积分fnint样条函数求积分 其中quad()与quad8()这两者的调用格式

19、完全一致，只不过在函数quad8（）中tol的默认值较低。该算法可以更加精确的求出积分的值，但现在在MATLAB一般用新引进的函数quadl（）来代替quad8(). 现在在来介绍下其他的方法3. LOBATTO求积 用quadl函数进行Lobatto求积，该函数的语法格式为：              q=quad(fun,a,b);求函数fun从a到b的积分近似，误差为le-6.fun为M文件函数或匿名函数的句柄，函数x=fun(x)应该接受矢量变量x并返回结果

20、矢量y。              q=quad(fun,a,b):使用绝对误差容限tol代替默认的容限1.0e-6。Tol的值更大时，需要的计算次数少，计算快。              q=quad(fun,a,b,tol,trace):在递归过程中使用非0的trace参数。     

21、         q.fcnt=quad(): 返回函数计算次数 例： 用quadl函数求下面积分的数值积分。    首先编写函数的M文件myfun.m。 Function y=myfun(x) Y=1./(x.3-2*x-5); 然后在命令窗口键入 Q=quadruple（myfun.0,2） Q= -0.4605 也可以直接使用匿名函数形式，如： F=（x）1.9(x.3-2*x-5)； Q=quad(F,0,2)  二重积分用dblquad函数对二重积分进行数值计算，该

22、函数的语法格式为：              q=dblquad(fun,xmin,xmax,ymin,ymax);调用quad函数对fun(x,y)函数进行二重积分运算积分区间为xmianxxmax;yminyymax。Fun参数是一个M文件函数或匿名函数的句柄。Fun(x,y)函数必须接受矢量x与标量y，并返回积分值组成的矢量。           &#

23、160;  q=dblquad(fun,xmin,xmax,ymin,ymax,tol):用容限tol代替默认的1.0e-6.              Q=dblquad(fun,xmin,xmax,ymin,ymax,tol,method):用method参数指定积分方法，代替默认的quad函数。Method还可以指定为quadl或用户指定的函数句柄，该函数必须与quad和quadl具有相同的调用格式。考虑下面的二重积分问题：  使用MATL

24、AB提供的dblquad函数就直接可以求出上述积分的数值解，该函数的调用格式为：I=dblquad(f,a,b,c,d,tol.trace)该函数求f(x,y)在a,b×c,d区域上的二重积分。 例：用dblquad函数求下面二重积分的数值近似。   首先建立函数的M文件integrnd.m,代码为： Function z=integrnd(x,y) z=y*sin(x)+x*cos(y); 然后在命令窗口键入 Q=dblquad(integend,pi,2*pi,0,pi) Q= -9.8696   例： 计算二次定积分 （1）  

25、  建立一个函数文件fxy.mFunction f=fxy(x,y)Global ki;Ki=ki+I;F=exp(-x.2/2).*sin(x.2+y);(2)      调用dblquad函数求解Global ki;ki=0;I=dblquad(fxy,-2,2,-1,1)KiI= 1.57449318974494Ki= 1038如果使用inline函数，则命令如下： f=inline(exp(-x.2/2).*sin(x.2+y),x,y); I=dblquad(f,-2,2m-1,1) I= 1.57449318974494  &