《数值分析课程设计》
报 告
专业： 学号： 学生姓名： 指导教师：
一、题目
数值积分中二重积分探究。
二、理论
数值积分就是用数值方法近似计算定积分?f(x)dx。其原理很简单，就是
ab将积分核f(x)用插值多项式Pn(x)替代，用多项式?Pn(x)dx的结果近似定积分
ab?baf(x)dx的值。
一般常用的方法是，将积分区间?a,b?等分为n个子区间，即取步长
h?(b?a)/n，
子区间端点为xk?a?kh（k=0,1,…，n），在每个子区间上套用插值积分公式，再将个区间的结果累加起来。
比较常用的有梯形公式，是在每个子区间上用1阶多项式（即直线段）近似
f(x)并积分的结果：
n?1hn?1h?? Tn???f(xk)?f(xk?1)???f(a)?2?f(xk)?f(b)?
2k?02?k?0?另外一种在实际应用中很受欢迎的方法是，在每个子区间上用2阶多项式（即抛物线）近似f(x)并积分，得到著名的辛普森（Simpson）公式：
hn?1Sn??[f(xk)?f4xk(?16k?0/n?1h2?)fxk?(1?)]fa[?(?)fx4k?6k?01/2(??f)xk?0n?1k2?fb()()]
其中xk?1/2?a?(k?1/2)h。
三、方法、算法与程序设计 Ⅰ.辛普森公式求二重积分
考虑二重积分?f(x,y)dA，它是曲面z?f(x,y)与平面区域R围成的体积，
R对于矩形区域R?{(x,y)|a?x?b,c?y?d}，可将它写成累次积分
)???(??f(x,ydxRacbdfx(y,dy)dx。)
若用复合辛普森公式，可分别将[a,b],[c,d]分成N,M等份，步长
- 2 -
h?b?ad?c,k?， NMd1先对积分?f(x,y)dy，应用复合辛普森公式，令yi?c?ik，yi?1/2?c?(i?)k，
c2则
?从而得
dcM?1M?1kf(x,y)dy?[f(x,y0)?4?f(x,yi?1/2)?2?f(x,yi)?f(x,yM)],
6i?0i?1??abdcM?1bM?1bbk?b?f(x,y)dydx???f(x,y0)dx?4??f(x,yi?1/2)dx?2??f(x,yi)dx??f(x,yM)dx?aaa6?ai?0i?1?。对每个积分再分别用复合辛普森公式即可求得积分值。MATLAB程序见附录
1，MATLAB中自带自适应辛普森公式dblquad()，对于变量区域同样适用。 对于变量区域R?{(x,y)|a?x?b,c(x)?y?d(x)}，写成累次积分的形式：
I???f(x,y)dxdy??{?Rabd(x)c(x)f(x,y)dy}dx
进行数值计算的表达式为：
I(a,b,c(x),d(?x?))m?1Mm?wnn?1Nvmf(xm n,y)上面的表达式中wm、vn表示权重，取决于一维积分方法。我们常用复合辛普森公式，先对内积分进行计算，在计算外积分，与矩形区域情况基本一致。 Ⅱ高斯求积公式求二重积分
在高斯求积公式中，若取权函数?(x)?1，区间为[?1,1]，则得公式
?1?1f(x)dx??Akfx(k。)
k?0n勒让德多项式是区间[?1,1]上的正交多项式，因此，勒让德多多项式Pn?1(x)的零点就是求积公式的高斯点。
若取P1(x)?x的零点x0?0做节点构造求积公式
?1?1；f(x)d?x2f(0)
11若取P2(x)?(3x2?1)的零点?构造求积公式
23
?1?1f(x)dx?f(?11)?f()； 33 - 3 -
当n?4时，求积公式为
?1?1f(x)dx?0.2369269f(?0.9061798)?0.4786287f(?0.5384693)?0.5688889f(0)?0.4786287f(0.5384693)?0.2369269f(0.9061798) 同样先用高斯求积公式求内积分，再求外积分，可得二重积分值。
四、算例、应用实例
算例：
计算二重积分??e?xydxdy。
D（1）若区域D?{0?x?1,0?y?1}，试分别用复合辛普森公式（取n=4）及高斯求积公式（取n=4）求积分。
（2）若区域D?{x2?y2?1:x?0,y?0}用复合辛普森公式（取n=4）求积分。 解： （1）
1357113?*y?*y?y*?y*?*y?*y?*yh1?0*y?1*y?0?0edxdy?6?0[e?4e8?4e8?4e8?4e8?2e4?2e2?2e4?e]dy11?xy??eD?xydxdy=
对各个积分应用复合辛普森公式。
也可应用MATLAB中的函数进行计算，程序见附录2。
先将区域D?{0?x?1,0?y?1}变换为区域D?{(u,v)|?1?u,v?1}，其中 u?2x?1,v?2y?1，等价于x? I??111(u?1),y?(v?1)，有 2211?(u?1v)?(4?1?110?10e?xydxdy???edudv。
1)对于u,v取n?2时的高斯求积公式节点及系数，即
u0?v0??0.9061798，u1?v1??0.5384693，u2?v2?0，u3?v3?0.5384693， u4?v4?0.9061798，A0?A4?0.2369269，A1?A3?0.4786287， A2?0.5688889
用n?4的高斯积分公式计算积分I，
I??（2）
1?1?1?1e1?(u?1)(v?1)4dudv???AiAeji?0j?0441?(u?1)(v?1)4
??edxdy??D?xy100?1?y2e?xydydx，
- 4 -
1131?y2?[0,1]等分为4等份，对应值为0,,,,1的yk(k?0,1,2,3,4)值，用
424yk,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),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;
for i=0:n-1 for 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(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}); end end end
- 5 -
q=((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);
- 6 -
