第二，启动MATLAB，新建脚本(Ctrl+N)，输入如下代码：
close all; clear all; clc
figure('Position',[50,50,600,500],'Name','double integration')
f=@(x,y) x*sin(y)-cos(x)+y-3;
Q1=dblquad(f,0,2*pi,-pi,pi,1.0e-3)
x=linspace(0,2*pi);y=linspace(-pi,pi);
ff=x.*sin(y)-cos(x)+y-3;
plot3(x,y,ff,'r','LineWidth',3)
grid on; view(35,35)
xlabel('x','FontSize',15); ylabel('y','FontSize',15)
其中Q1=dblquad(f,0,2*pi,-pi,pi,1.0e-3)采用默认方法quad计算二重积分，绝对计算精度设为1.0e-3。
plot3(x,y,ff,'r','LineWidth',3)是绘制被积分函数ff=x.*sin(y)-cos(x)+y-3的图像。

	

点击上方蓝字  关注我们
从本节开始，我们将进入matlab的数值微积分与方程求解模块，一起学习如何利用matlab去解决微积分问题。
对于本节内容，主要分为两个部分讲解，数值微分和数值积分，那么下面，就开始今天的学习吧！
一、数值微分
在正式开始之前，有几个新概念需要讲解一下(1) 数值差分与差商：微积分中 ，任意函数 f(x) 在x0 点的导数是通过极限定义的(如下图所示)：
如果去掉极限定义中h趋向于0的极限过程，得到函数在x0 点处以h( h>0 )为步长的向前差分、向后差分和中心差分公式(如下图)：
当步长h充分小时，得到函数在x0 点处以h( h>0 )为步长的向前差商、向后差商和中心差商公式(如下图)：
函数 f(x) 在 点x0的微分接近于函数在该点的差分，而f在点x的导数接近于函数在该点的差商。
(2)数值微分的实现：
MATLAB 提供了求向前差分的函数diff ，其调用格式有三种：
dx=diff(x) ：计算向量x的向前差分， dx(i)=x(i+1)- x(i) ， i=1，2， … ，n, n- 1 。
dx=diff( x,n) ) ：计算向量x的n阶向前差分。例如，diff(x,2)=diff(diff(x)) 。
dx=diff( A,n,dim) ) ：计算矩阵A的n阶差分， dim=1 时(默认状态)，按列计算差分；dim=2 ，按行计算差分。
注意 ：diff 函数计算的是向量元素间的差分，故差分向量元素的个数比原向量少了一个。同样，对于矩阵来说，差分后的矩阵比原矩阵少了一行或一列。
另外，计算差分之后，可以用 f(x) 在某点处的差商作为其导数的近似值。
例、设f(x)=ln(x)，采用上述方法绘制出[0，10]之间的其导函数的图像，求解出导函数在1处的值，对比理论值1之间的差值，看相差大不大。
>>  x=0:0.01:10;
y=log(x);
f1=diff(y)./diff(x);
plot(x(1:end-1),f1);
打开matlab的工作区查看f1在100处的值，即就是导函数在1处的值(如图)：
值为1.005，与理论值1相差0.005，可能是由于x=0:0.01:10;这一语句导致，发现当取点越多时，结果越接近1，再看图像，其符合1/x的函数图像。
二、 数值积分
1、数值积分基本原理在高等数学中，计算定积分依靠微积分基本定理，只要找到被积函数 f(x)的原函数F(x)，则可用 牛顿 — 莱布尼兹 ( Newton- - Leibniz ) 公式：
在有些情况下，应用牛顿—莱布尼兹公式有困难，例如，当被积函数的原函数无法用初等函数表示，或被积函数是用离散的表格形式给出的。这时就需要用数值解法来求定积分的近似值 。
求定积分的数值方法多种多样，如梯形法、辛普森( Simpson )法 、高斯求积公式等。它们的基本思想都是将积分区间 [a ， b] 分成n个子区间 [xi，xi+1] ，i=1 ，2 ，… ，n，其中x1 =a ，xn+1=b ，这样求定积分问题就分解为下面的求和问题。
在 每一个小的子区间上定积分的值可以近似求得, , 从而避免了 牛顿 — 莱布尼兹公式需要寻求原函数的困难 。
2、数值积分的实现(1)基于自适应辛普森方法：[ I,n ]=quad( filename,a,b,tol,trace)
(2)基于自适应 Gauss- - Lobatto方法：[ I,n ]= quadl ( filename,a,b,tol,trace)其中， filename 是被积函数名；a和b分别是定积分的下限和上限，积分限 [a ，b]必须是有限的，不能为无穷大( Inf )；tol 用来控制积分精度，默认时取tol =10^(- 6 )；
trace 控制是否展现积分过程，若取非0则展现积分过程，取0则不展现，默认时取 trace=0 ；返回参数I即定积分的值，n为被积函数的调用次数。
例、分别用 quad 函数和 quadl 函数求定积分的近似值，并在相同的积分精度下，比较被积函数的调用次数。求解在[1，10]上，sin(x)的定积分值。
>>  format long
f=@(x)cos(x);
[ I,n ]=quad(f,1,10,1e-8)
I =
  -1.385492095697597
n =
   229
>> [ I,n ]=quadl(f,1,10,1e-8)
I =
  -1.385492095687089
n =
   138
>> sin(10)-sin(1)
ans =
  -1.385492095697266
可以看出，两种方式求出的定积分值与理论值sin(10)-sin(1)几乎一致，但两种函数求解时的被积函数的调用次数却有差别。
(3)基于全局自适应积分方法：I=integral( filename,a,b )其中，I是计算得到的积分；filename 是被积函数；a和b分别是定积分的下限和上限，积分限可以为无穷大 。
(4)基于自适应高斯- - 克朗罗德方法[ I,err ]= quadgk( filename,a,b)其中 ，err返回近似误差范围，其他参数的含义和用法与quad函数相同。积分上下限可以是无穷大(−Inf 或 Inf)，也可以是复数。如果积分上下限是复数，则 quadgk函数在复平面上求积分。
上述四种积分的求解方法各有不同，同时也各有优劣，一、二两种方法要求积分上下限是有限的，不可取无穷值，三、四两种方法中积分上下限可以取无穷值，且第四种方法还可以在复平面上求积分。
下面介绍一种特殊的积分方法：梯形积分法I= trapz( x,y)其中，向量x、y 定义函数关系y=f(x)，值得注意的是，此处的x，y表示为向量形式。例、 5 设 x=1:10，y=[5,4,1,9,3,3,5,8,7,2], 用 trapz函数计算定积分并且画图。
>> x=1:10;y=[5,4,1,9,3,3,5,8,7,2];
>> plot(x,y);
>> I1=trapz(x,y)
I1 =
  43.500000000000000
三、 多重定积分的数值求解
现将多重积分的求解函数总结如下，它们的用法和上述的一重积分用法基本一致求二重积分的数值解：
I=integral2( filename,a,b,c,d)I=quad2d( filename,a,b,c,d)I= dblquad( filename,a,b,c,d,tol)求三重积分的数值解：
I=integral3( filename,a,b,c,d,e,f) )I= triplequad( filename,a,b,c,d,e,f,tol)
本节内容就讲述到这里，下节将讲述线性方程组求解，敬请期待！
关于MATLAB的学习：
大家可以关注我们的知乎专栏——数据可视化和数据分析中matlab的使用：
https://zhuanlan.zhihu.com/c_1131568134137692160
欢迎大家加入我们的MATLAB学习交流群：
953314432
▼往期精彩回顾▼Matlab数据统计与分析(四)——参数估计数据预处理小结MATLAB中的非线性空间滤波基础MATLAB中的图像变换之线性空间滤波数学建模与MATLAB--插值与拟合扫码关注我们更多精彩等待你发现出品：Asoul水云天课堂工作室好看你就点点我

二重积分的数值方法
《数值分析课程设计》
报  告
专业：
学号：
学生姓名：
指导教师：
一、题目
数值积分中二重积分探究。
二、理论
数值积分就是用数值方法近似计算定积分。其原理很简单，就是将积分核用插值多项式替代，用多项式的结果近似定积分的值。
一般常用的方法是，将积分区间等分为个子区间，即取步长
，
子区间端点为(k=0,1,…，n)，在每个子区间上套用插值积分公式，再将个区间的结果累加起来。
比较常用的有梯形公式，是在每个子区间上用1阶多项式(即直线段)近似并积分的结果：
另外一种在实际应用中很受欢迎的方法是，在每个子区间上用2阶多项式(即抛物线)近似并积分，得到著名的辛普森(Simpson)公式：
其中。
三、方法、算法与程序设计
Ⅰ.辛普森公式求二重积分
考虑二重积分，它是曲面与平面区域R围成的体积，对于矩形区域，可将它写成累次积分
。
若用复合辛普森公式，可分别将,分成N,M等份，步长
,，
先对积分，应用复合辛普森公式，令，，则
从而得
。对每个积分再分别用复合辛普森公式即可求得积分值。MATLAB程序见附录1，MATLAB中自带自适应辛普森公式dblquad()，对于变量区域同样适用。
对于变量区域，写成累次积分的形式：
进行数值计算的表达式为：
上面的表达式中、表示权重，取决于一维积分方法。我们常用复合辛普森公式，先对内积分进行计算，在计算外积分，与矩形区域情况基本一致。
Ⅱ高斯求积公式求二重积分
在高斯求积公式中，若取权函数，区间为，则得公式
。
勒让德多项式是区间上的正交多项式，因此，勒让德多多项式的零点就是求积公式的高斯点。
若取的零点做节点构造求积公式
；
若取的零点构造求积公式
；
当时，求积公式为
同样先用高斯求积公式求内积分，再求外积分，可得二重积分值。
四、算例、应用实例
算例：
计算二重积分。
(1)若区域，试分别用复合辛普森公式(取n=4)及高斯求积公式(取n=4)求积分。
(2)若区域用复合辛普森公式(取n=4)求积分。
解：
(1)
=
对各个积分应用复合辛普森公式。
也可应用MATLAB中的函数进行计算，程序见附录2。
先将区域变换为区域，其中
，等价于，有
。
对于取时的高斯求积公式节点及系数，即
，，，，
，，，
用的高斯积分公式计算积分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})+...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-1for j=0:m-1x=a+2*i*(A-a)/2/n;y=b+2*j*(B-b)/2/m;?????????

我有点困在一个函数上，我试图通过scipy，python进行数值积分。在
为了简单起见，我将函数定义为：integral f(x,y)= SUM[double integral(ax+by)dxdy]
a和b是常数，但对于每个积分的方程，它们是不同的。我已经将每个函数分别积分，然后将所有积分的结果求和，但是这需要花费大量的时间来计算，这对于我正在尝试实现的目标并不理想。在
是否有一种方法可以通过扩展和来同时集成整个函数：
^{pr2}$
然后将带有(a,b)元组列表的函数传递给scipy的dblquad函数？在
目前，我正在努力在文献中找到任何与此相关的东西。在
*编辑
我提供了一个示例代码，以更清楚地显示我想要实现的目标：import sys
import re
import math
from scipy.integrate import dblquad
def f((x,y),variables):
V=0
for v in variables:
a,b=v
V=V+ax+by
return (V)
def integral(x_max,y_max,variables):
return dblquad(f, 0, y_max, lambda x: 0, lambda x: x_max,args=variables)
def main():
variables=[(1,2),(3,4),(5,6)] #example variables. The length of this list can change with the code I am running.
x_max=y_max=1
integral(x_max,y_max,variables)
if __name__ == '__main__':
main()
因此，返回的错误是：Traceback (most recent call last):
File "integration_example.py", line 23, in 
main()
File "integration_example.py", line 19, in main
integral(x_max,y_max,variables)
File "integration_example.py", line 14, in integral
return dblquad(f, 0, y_max, lambda x: 0, lambda x: x_max,args=variables)
File "/usr/lib/python2.7/dist-packages/scipy/integrate/quadpack.py", line 435, in dblquad
return quad(_infunc,a,b,(func,gfun,hfun,args),epsabs=epsabs,epsrel=epsrel)
File "/usr/lib/python2.7/dist-packages/scipy/integrate/quadpack.py", line 254, in quad
retval = _quad(func,a,b,args,full_output,epsabs,epsrel,limit,points)
File "/usr/lib/python2.7/dist-packages/scipy/integrate/quadpack.py", line 319, in _quad
return _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit)
File "/usr/lib/python2.7/dist-packages/scipy/integrate/quadpack.py", line 382, in _infunc
myargs = (x,) + more_args
TypeError: can only concatenate tuple (not "list") to tuple
很明显，函数不喜欢我用我写这篇文章的方式传递一系列要放入积分中的值。有办法吗？
(抱歉，这可能是更好的措辞方式)。在