二重积分的数值方法

《数值分析课程设计》报告专业：学号：学生姓名：指导教师：一、题目数值积分中二重积分探究。二、理论数值积分就是用数值方法近似计算定积分()bafxdx。其原理很简单，就是将积分核()fx用插值多项式()nPx替代，用多项式()bnaPxdx的结果近似定积分()bafxdx的值。一般常用的方法是，将积分区间,ab等分为n个子区间，即取步长()/hban，子区间端点为kxakh（k=0,1,…，n），在每个子区间上套用插值积分公式，再将个区间的结果累加起来。比较常用的有梯形公式，是在每个子区间上用1阶多项式（即直线段）近似()fx并积分的结果：11100()()()2()()22nnnkkkkkhhTfxfxfafxfb另外一种在实际应用中很受欢迎的方法是，在每个子区间上用2阶多项式（即抛物线）近似()fx并积分，得到著名的辛普森（Simpson）公式：1111/211/2000[()4()()][()4()2()()]66nnnnkkkkkkkkhhSfxfxfxfafxfxfb其中1/2(1/2)kxakh。三、方法、算法与程序设计Ⅰ.辛普森公式求二重积分考虑二重积分(,)RfxydA，它是曲面(,)zfxy与平面区域R围成的体积，对于矩形区域{(,)|,}Rxyaxbcyd，可将它写成累次积分(,)((,))bdacRfxydxfxydydx。若用复合辛普森公式，可分别将[,]ab,[,]cd分成N,M等份，步长bahN,dckM，先对积分(,)dcfxydy，应用复合辛普森公式，令iycik，1/21()2iycik，则1101/201(,)[(,)4(,)2(,)(,)],6MMdiiMciikfxydyfxyfxyfxyfxy从而得1101/201(,)(,)4(,)2(,)(,)6MMbdbbbbiiMacaaaaiikfxydydxfxydxfxydxfxydxfxydx。对每个积分再分别用复合辛普森公式即可求得积分值。MATLAB程序见附录1，MATLAB中自带自适应辛普森公式dblquad()，对于变量区域同样适用。对于变量区域{(,)|,()()}Rxyaxbcxydx，写成累次积分的形式：()()(,){(,)}bdxRacxIfxydxdyfxydydx进行数值计算的表达式为：11(,,(),())(,)NnMmnmmnmIabcxdxwvfxy上面的表达式中mw、nv表示权重，取决于一维积分方法。我们常用复合辛普森公式，先对内积分进行计算，在计算外积分，与矩形区域情况基本一致。Ⅱ高斯求积公式求二重积分在高斯求积公式中，若取权函数()1x，区间为[1,1]，则得公式110()()nkkkfxdxAfx。勒让德多项式是区间[1,1]上的正交多项式，因此，勒让德多多项式1()nPx的零点就是求积公式的高斯点。若取1()Pxx的零点00x做节点构造求积公式11()2(0)fxdxf；若取221()(31)2Pxx的零点13构造求积公式1111()()()33fxdxff；当4n时，求积公式为11()0.2369269(0.9061798)0.4786287(0.5384693)0.5688889(0)0.4786287(0.5384693)0.2369269(0.9061798)fxdxfffff同样先用高斯求积公式求内积分，再求外积分，可得二重积分值。四、算例、应用实例算例：计算二重积分xyDedxdy。（1）若区域{01,01}Dxy，试分别用复合辛普森公式（取n=4）及高斯求积公式（取n=4）求积分。（2）若区域22{1:0,0}Dxyxy用复合辛普森公式（取n=4）求积分。解：（1）xyDedxdy=1357113****111***0*1*8888424000[4444222]6yyyyyyyxyyyhedxdyeeeeeeeeedy对各个积分应用复合辛普森公式。也可应用MATLAB中的函数进行计算，程序见附录2。先将区域{01,01}Dxy变换为区域{(,)|1,1}Duvuv，其中21,21uxvy，等价于11(1),(1)22xuyv，有11111(1)(1)40011uvxyIedxdyedudv。对于,uv取2n时的高斯求积公式节点及系数，即000.9061798uv，110.5384693uv，220uv，330.5384693uv，440.9061798uv，040.2369269AA，130.4786287AA，20.5688889A用4n的高斯积分公式计算积分I，114411(1)(1)(1)(1)441100uvuvijijIedudvAAe（2）21100yxyxyDedxdyedydx，21[0,1]y等分为4等份，对应值为1130,,,,1424的(0,1,2,3,4)kyk值，用1,kkyy节点应用辛普森公式对内积分求积，再用复合辛普森公式对外积分求积，也可用MATLAB中的函数实现，结果和程序如下（附录3）。五、参考文献【1】数值分析李庆扬，王朝能，易大义清华大学出版社【2】数值分析课程设计陈越，童若锋浙江大学出版社【3】MATLAB教程张志涌北京航空航天大学出版社六、附录附录1：functionq=DblSimpson(f,a,A,b,B,m,n)if(m==1&&n==1)%辛普森公式q=((B-b)*(A-a)/9)*(subs(sym(f),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(f),findsym(sym(f)),{(A-a)/2,(B-b)/2}));else%复合辛普森公式q=0;fori=0:n-1forj=0:m-1x=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(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('sqrt(1-y.^2)','y');Q=quad2dggen(fun,xlower,xupper,0,1,1e-4);