《科学计算引论》实验课：第三章 数值微积分

1、科学计算引论科学计算引论实验实验第三章第三章 数值微积分数值微积分本章内容简介 3.1 3.1 数值微分数值微分 3.2 3.2 机械求积公式机械求积公式 3.3 3.3 Newton-Cotes 公式及其复合求积法公式及其复合求积法 3.4 3.4 变步长求积法变步长求积法 3.5 3.5 Gauss 求积公式求积公式数值微分qMatlab 中差分的命令：中差分的命令：diff调用格式为：调用格式为：DX=diff(X) %计算向量计算向量X的向前差分，的向前差分，DX(i)=X(i+1)-X(i)，i=1,2,n-1。DX=diff(X,n) %计算计算X的的n阶向前差分。例如，阶向前差分

2、。例如，diff(X,2)=diff(diff(X)。DX=diff(A,n,dim) %计算矩阵计算矩阵A的的n阶差分，阶差分，dim=1时时(缺省状态缺省状态)，按列计算，按列计算差分；差分；dim=2，按行计算差分。，按行计算差分。数值微分例例：生成以向量：生成以向量V=1,2,3,4,5,6为基础的范得蒙矩阵，按为基础的范得蒙矩阵，按列进行差分运算。列进行差分运算。解解：命令如下：命令如下： V=vander(1:6) DV=diff(V) %计算计算V的一阶差分的一阶差分 DV = 31 15 7 3 1 0 211 65 19 5 1 0 781 175 37 7 1 0 2101

3、 369 61 9 1 0 4651 671 91 11 1 0数值微分例例：设：设 用不同的方法求函数用不同的方法求函数f(x)的数值导数，并在同一个坐标系中做的数值导数，并在同一个坐标系中做出出f (x)的图像。的图像。解解：用三种方法：用三种方法：1、用一个、用一个5次多项式次多项式p(x)拟合函数拟合函数f(x)，并对，并对p(x)求一般意义下求一般意义下的导数的导数dp(x)，求出，求出dp(x)在假设点的值；在假设点的值；2、直接求、直接求f(x)在假设点的数值导数；在假设点的数值导数；3、求出、求出g(x)=f (x)的表达式的表达式 然后求然后求f (x)在假设点的导数值。在假

4、设点的导数值。326( )212552f xxxxxx232563411( )522126 (5)xxfxxxxx 数值微分程序如下：程序如下：f=inline(sqrt(x.3+2*x.2-x+12)+(x+5).(1/6)+5*x+2);g=inline(3*x.2+4*x-1)./sqrt(x.3+2*x.2-x+12)/2+1/6./(x+5).(5/6)+5);x=-3:0.01:3;p=polyfit(x,f(x),5); %1用用5次多项式次多项式p拟合拟合f(x)dp=polyder(p); %对拟合多项式对拟合多项式p求导数求导数dpdpx=polyval(dp,x); %求

5、求dp在假设点的函数值在假设点的函数值dx=diff(f(x,3.01)/0.01; %2直接对直接对f(x)求数值导数求数值导数gx=g(x); %3求函数求函数f的导函数的导函数g在假设点的导数在假设点的导数plot(x,dpx,g.,x,dx,b.,x,gx,r-); %作图作图数值微分图如下：图如下：数值积分：复合梯形法r Matlab 梯形积分函数：梯形积分函数：lZ=trapz(x,y) %计算计算y 对对x的梯形积分，其中的梯形积分，其中x、y定义函数关系定义函数关系y=f(x)。例例: 用梯形法求下面定积分。用梯形法求下面定积分。20sinVx dx解解：梯形法求定积分：梯形法

6、求定积分：x=0:pi/1000:pi;%细分积分区间细分积分区间y=sin(x.2); %求相应函数值求相应函数值v=trapz(x,y) %用用Matlab专门的梯形积分函数专门的梯形积分函数% 可以看出矩形法和梯形法都是利用可以看出矩形法和梯形法都是利用离散点列离散点列来求相应的定来求相应的定积分的，它们更适合实际问题中的积分的，它们更适合实际问题中的离散的表格函数离散的表格函数。数值积分：自适应Simpson法r 自适应自适应Simpson法法： l I,n=quad(fname,a,b,tol )ufname是被积函数名是被积函数名, a和和b分别是定积分的下限和上限。分别是定积分的

7、下限和上限。u tol用来控制积分精度，缺省时取用来控制积分精度，缺省时取tol=0.001。u I为定积分值为定积分值, n为被积函数的调用次数。为被积函数的调用次数。l 注意注意：调用调用quad函数时，先要建立一个描述被积函数函数时，先要建立一个描述被积函数f(x)的函数文件、的函数文件、内联函数或匿名函数。内联函数或匿名函数。r 例例 求求sin(x2)在在0,pi上的积分。上的积分。 解解： fun=inline(sin(x.2),x); % 内联函数内联函数 quad(fun,0,pi) %默认误差精度默认误差精度0.001, 也可自由设定也可自由设定: Q,n=quad(fun,

8、0,pi,1e-6)r 方法方法2 2：用匿名函数：用匿名函数 fun=(x)sin(x.2) fun=(x)sin(x.2) Q,nQ,n=quad(fun,0,pi)=quad(fun,0,pi)数值积分：自适应Lobatto法r 自适应自适应Lobatto法法： l I,n=quadl(fname,a,b,tol )ufname是被积函数名是被积函数名, a和和b分别是定积分的下限和上限。分别是定积分的下限和上限。u tol用来控制积分精度，缺省时取用来控制积分精度，缺省时取tol=10(-6)。u I为定积分值为定积分值, n为被积函数的调用次数。为被积函数的调用次数。l 同样：调用同

9、样：调用quadl函数时，先要建立一个描述被积函数函数时，先要建立一个描述被积函数f(x)的函数文件、内的函数文件、内联函数或匿名函数。自适应联函数或匿名函数。自适应Lobatto法比自适应法比自适应Simpson法要更高效、精确。法要更高效、精确。r例例：求：求sin(x2)在在0,pi上的积分。上的积分。 fun=inline(sin(x.2),x); % 内联函数内联函数 Q,n=quadl(fun,0,pi) %方法方法2：用匿名函数：用匿名函数 fun=(x)sin(x.2) % 匿名函数匿名函数 quadl(fun,0,pi)数值积分：二重积分r 求二重积分：求二重积分： l I=

10、dblquad (fname,a,b,c,d,tol)u fname是被积函数名是被积函数名u a、b是内积分的下上限是内积分的下上限, c、d是外积分的下上限。是外积分的下上限。u tol用来控制积分精度，缺省时取用来控制积分精度，缺省时取tol=10(-6)。u I为定积分值。为定积分值。数值积分：二重积分r 例例：求二重积分：求二重积分：解解：r %方法方法1：用内联函数：用内联函数r fun=inline(x+y2,x,y); r q=dblquad(fun,-1, 1,0,2)r %方法方法2: 用匿名函数用匿名函数r fun=(x,y)x+y2; r q=dblquad(fun,-

11、1, 1,0,2)12120()dxxydyMatlab积分函数计算符号积分计算符号积分 更精确更精确int(int(f,v,a,bf,v,a,b) ) int(int(f,a,bf,a,b) ) int(int(f,vf,v) )， int(int(f f) )计算计算f f关于关于变量变量v v的在的在 a,ba,b 上定积分上定积分计算计算默认变量默认变量x x的定积分的定积分计算计算f f关于关于变量变量v v或默认变量或默认变量 的不定积分的不定积分计算数值积分计算数值积分z=z=trapz(x,ytrapz(x,y) )梯形积分法：其中：梯形积分法：其中：x x表示积分区间的离散化

12、变量表示积分区间的离散化变量y y为相应于为相应于x x的被积函数的值的被积函数的值此方法，精度低，适用于数值函数和光滑性不好的函此方法，精度低，适用于数值函数和光滑性不好的函数数z=z=quad(fun,a,b,tolquad(fun,a,b,tol) )quad()quad()使用的是自适应步长使用的是自适应步长SimpsonSimpson积分法积分法, ,它比它比quadlquadl()()的精度要低，不能用于广义积分，对一些假奇的精度要低，不能用于广义积分，对一些假奇异积分也不能直接求解，需要适当整理变换。异积分也不能直接求解，需要适当整理变换。funfun表示被积函数名，表示被积函数名，a,ba,b是积分上下限，是积分上下限，toltol表示精度表示精度z=z=quadl(fun,a,b,tolquadl(fun,a,b,tol) )quadl()quadl()使用的是自适应步长使用的是自适应步长LobattoLobatto积分法积分法z=z=