二重积分的数值方法

1、数值分析课程设计报 告 专业： 学号： 学生姓名： 指导教师： 一、题目 数值积分中二重积分探究。二、理论数值积分就是用数值方法近似计算定积分。其原理很简单，就是将积分核用插值多项式替代，用多项式的结果近似定积分的值。一般常用的方法是，将积分区间等分为个子区间，即取步长，子区间端点为（k=0,1,，n），在每个子区间上套用插值积分公式，再将个区间的结果累加起来。比较常用的有梯形公式，是在每个子区间上用1阶多项式（即直线段）近似并积分的结果： 另外一种在实际应用中很受欢迎的方法是，在每个子区间上用2阶多项式（即抛物线）近似并积分，得到著名的辛普森（Simpson）公式： 其中。 三、方法、算法与

2、程序设计.辛普森公式求二重积分 考虑二重积分，它是曲面与平面区域R围成的体积，对于矩形区域，可将它写成累次积分 。若用复合辛普森公式，可分别将,分成N,M等份，步长,，先对积分，应用复合辛普森公式，令，则从而得 。对每个积分再分别用复合辛普森公式即可求得积分值。MATLAB程序见附录1，MATLAB中自带自适应辛普森公式dblquad()，对于变量区域同样适用。对于变量区域，写成累次积分的形式：进行数值计算的表达式为： 上面的表达式中、表示权重，取决于一维积分方法。我们常用复合辛普森公式，先对内积分进行计算，在计算外积分，与矩形区域情况基本一致。高斯求积公式求二重积分在高斯求积公式中，若取权函

3、数，区间为，则得公式 。勒让德多项式是区间上的正交多项式，因此，勒让德多多项式的零点就是求积公式的高斯点。 若取的零点做节点构造求积公式 ；若取的零点构造求积公式 ；当时，求积公式为 同样先用高斯求积公式求内积分，再求外积分，可得二重积分值。 四、算例、应用实例算例：计算二重积分。（1）若区域，试分别用复合辛普森公式（取n=4）及高斯求积公式（取n=4）求积分。（2）若区域用复合辛普森公式（取n=4）求积分。解：（1） = 对各个积分应用复合辛普森公式。 也可应用MATLAB中的函数进行计算，程序见附录2。 先将区域变换为区域，其中 ，等价于，有 。对于取时的高斯求积公式节点及系数，即， 用的

4、高斯积分公式计算积分I， （2） ，等分为4等份，对应值为的值，用节点应用辛普森公式对内积分求积，再用复合辛普森公式对外积分求积，也可用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)+.sub

5、s(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,(

6、B-b)/2);else %复合辛普森公式q=0;for i=0:n-1for j=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