二重积分的数值方法

《二重积分的数值方法》由会员分享，可在线阅读，更多相关《二重积分的数值方法（6页珍藏版）》请在装配图网上搜索。

1、数值分析课程设计报告专业：-学号：_ 学生姓名:指导教师一、题目 数值积分中二重积分探究。二、理论数值积分就是用数值方法近似计算定积分Jbf (x)dx。其原理很简单，就是a将积分核f (x)用插值多项式Pn(x)替代，用多项式I bPn(x)dx的结果近似定积分aIb f(x)dx 的值。a一般常用的方法是，将积分区间a,b等分为n个子区间，即取步长h = (b - a) / n，子区间端点为xk = a + kh (k=0丄，n),在每个子区间上套用插值积分公式，再 将个区间的结果累加起来。比较常用的有梯形公式，是在每个子区间上用1阶多项式(即直线段)近似f(x) 并积分的结果：=2 艺/

2、(f(电+1)=2k=0f (a) + 2 艺 f (xk) + f (b)k=0另外一种在实际应用中很受欢迎的方法是，在每个子区间上用 2 阶多项式(即抛物线)近似f (x)并积分，得到著名的辛普森(Simpson)公式:S = 6 艺f (Xk) + 4 f (Xk +1/2)+ f (Xk +1) = 6 f (a) + 4 艺 f (Xk +1/2)+ 2 艺 f (Xk) + f(b) k =0k =0k =0其中 Xk +1/2 = a + (k +1/ 2)h。三、方法、算法与程序设计I .辛普森公式求二重积分 考虑二重积分I f (x,y)dA，它是曲面z = f (x,y)与

3、平面区域R围成的体积,R对于矩形区域R = (x, y) I a x b,c y d，可将它写成累次积分f ( x, y)dy)dx 。II f (x, y)dx = Ib (Idac R若用复合辛普森公式，可分别将a,b,c,d分成N,M等份，步长丫 b - a d - ch =, k =NM先对积分J df (x,y)dy，应用复合辛普森公式，令yi = c + ik , yi +1/2 = c + (i +丄)k , c2则Jd f (x, y )dy = f (x, y 0) + 4 S f (x, c6i=0y,+1/2)+2 S f (x，yJ+f (x，yM 几i=1从而得JbJ

4、df (x,y)dydx = 6 Jbf(x,y0)dx + 4SJbf(x,y, + 口2)dx + 2SJbf(x,y,)dx + Jbf(x,yM)dxa cU aaaai =0i =1。对每个积分再分别用复合辛普森公式即可求得积分值。MATLAB程序见附录1, MATLAB中自带自适应辛普森公式dblquad(),对于变量区域同样适用。对于变量区域R = (x, y) I a x b,c(x) y d(x)，写成累次积分的形式:I = JJ f (x, y)dxdy = JbJd(x) f (x, y)dydxRa c( x)进行数值计算的表达式为:I(a,b,c(x),d(x)=Sw

5、mSN vnf (xm, ymn)m =1 n=1上面的表达式中 wm、vn 表示权重，取决于一维积分方法。我们常用复合辛普森 mn公式，先对内积分进行计算，在计算外积分，与矩形区域情况基本一致。II高斯求积公式求二重积分在高斯求积公式中，若取权函数p (x) = 1，区间为-1,1，则得公式J1 f (x)dx 沁-1SAkf(xk)k=0勒让德多项式是区间-1,1上的正交多项式，因此，勒让德多多项式Pn +1( x) 的零点就是求积公式的高斯点。若取P1( x) = x的零点x 0 = 0做节点构造求积公式J1 f (x)dx 沁 2f (0)；-1若取 P2(x)旳2 -1)的零点土构造

6、求积公式J1 f (x)dx 沁 f (-占 + f(*)；当n二4时，求积公式为n f (x)dx u 0.2369269 f (0.9061798) + 0.4786287 f (0.5384693) + 0.5688889 f (0) + 0.4786287 f (0.5384693) + 0.2369269 f (0.9061798) 同样先用高斯求积公式求内积分，再求外积分，可得二重积分值。四、算例、应用实例算例：计算二重积分JJe-xydxdy。D（1）若区域D = 0 x 1,0 y 1，试分别用复合辛普森公式（取n=4）及高 斯求积公式（取n=4）求积分。（2）若区域D = x

7、2 + y2 0, y 0用复合辛普森公式（取n=4）求积分。解：（1）Df1 f1 e-xydxdy =00601357113e -0* y + 4e8*y + 4e8*y + 4e8*y + 4e8*y + 2e4*y + 2e2*y + 2e4*y+ e1*y dy对各个积分应用复合辛普森公式。也可应用MATLAB中的函数进行计算，程序见附录2。先将区域D = 0 x 1,0 y 1变换为区域D = (u, v)| 1 u, v 1，其中u = 2 x 1,v = 2 y 1，等价于 x = (u +1), y = (v +1)，有 22I = f1 f1 e-xydxdy = f1 f

8、1 e-4(u+1)(v+1)dudv。0 0 1 1对于u, v取n = 2时的高斯求积公式节点及系数，即u0 =v0 = 0.9061798 ， u1= v1= 0.5384693 ， u2 = v2 = 0， u3 = v3 = 0.5384693u4 = v4 = 0.9061798 ， A0= A4=0.2369269， A1= A3 = 0.4786287 ，A2 = 0.5688889用 n = 4的高斯积分公式计算积分 I，I=f1f1 e4(u+1)(v+1) dudv u 工区 AAe 4(u+1)(v+1)1 1=0 =0 i ji=0 j =000 D113点-y2 u

9、 0,1等分为4等份，对应值为 , 的yk(k = 0,1,2,3,4)值，用*4 2 4 kyk,yk +1节点应用辛普森公式对内积分求积，再用复合辛普森公式对外积分求积, 也可用MATLAB中的函数实现，结果和程序如下(附录3)。五、参考文献【1】 数值分析 李庆扬，王朝能，易大义 清华大学出版社【2】 数值分析课程设计 陈越，童若锋 浙江大学出版社【3】MATLAB教程 张志涌北京航空航天大学出版社六、附录附录1：function q=DblSimpson(f,a,A,b,B,m,n) if(m=1 & n=1) %辛普森公式 q=(B-b)*(A-a)/9)*(subs(sym(f),

10、findsym(sym(f),a,b)+. subs(sym(f),findsym(sym(f),a,B)+. subs(sym(f),findsym(sym(f),A,b)+. subs(sym(f),findsym(sym(f),A,B)+.4*subs(sym(f),findsym(sym(f),(A-a)/2,b)+. 4*subs(sym(f),findsym(sym(f),(A-a)/2,B)+.4*subs(sym(f),findsym(sym(f),a,(B-b)/2)+. 4*subs(sym(f),findsym(sym(f),A,(B-b)/2)+.16*subs(sym

11、(f),findsym(sym(f),(A-a)/2,(B-b)/2);else %复合辛普森公式q=0;for i=0:n-1for j=0:m-1 x=a+2*i*(A-a)/2/n; y=b+2*j*(B-b)/2/m;x1=a+(2*i+1)*(A-a)/2/n; y1=b+(2*j+1)*(B-b)/2/m;x2=a+2*(i+1)*(A-a)/2/n; y2=b+2*(j+1)*(B-b)/2/m;q=q+subs(sym(f),findsym(sym(f),x,y)+. subs(sym(f),findsym(sym(f),x,y2)+. subs(sym(f),findsym(

12、sym(f),x2,y)+.subs(sym(f),findsym(sym(f),x2,y2)+. 4*subs(sym(f),findsym(sym(f),x,y1)+.4*subs(sym(f),findsym(sym(f),x2,y1)+. 4*subs(sym(f),findsym(sym(f),x1,y)+.4*subs(sym(f),findsym(sym(f),x1,y2)+. 16*subs(sym(f),findsym(sym(f),x1,y1);endendendq=(B-b)*(A-a)/36/m/n)*q; 附录2：fun=inline(exp(x*y); d=dblquad(fun,0,1,0,1);附录3：f = inline(exp(-x*y),x,y); xlower = inline(0,y);xupper = inline(sqr t(l-y.2),y);Q=quad2dggen(fun,xlower,xupper,0,1,1e-4);