《数学实验 第4版》课件 4.3 MATLAB数值积分与微分

第四章数值分析

实验4.1插值

实验4.2离散数据的曲线拟合数学实验实验4.3MATLAB数值积分与微分实验4.4常微分方程的数值解2实验4.3MATLAB数值积分与微分一、数值积分二、数值微分实验4.3MATLAB数值积分与微分3一、数值积分实验4.3MATLAB数值积分与微分实验目的通过本实验了解数值积分和数值微分的方法，会用MATLAB进行数值积分和数值微分.数值积分也称为数值求积，是求定积分的近似值的数值方法.数值积分算法的发展与完善可以追溯到17世纪，从最早的牛顿-柯特斯公式到如今的自适应积分法和高维积分方法，经历了多个阶段不断的改进.4实验4.3MATLAB数值积分与微分

定积分是微积分中的基本计算方法，但在很多实际问题中,经常会遇到被积函数的原函数不能用初等函数表示；或虽然能找到原函数但因其很复杂而难以给出最后的积分结果；或被积函数以数表的形式给出，因此求定积分的数值解在实际中应用显得特别重要．数学家们通过不断的努力和创新，经过了几个世纪的发展与完善，提高了数值积分算法的精度和效率，为科学计算和工程应用提供了重要的支持.5用数值方法近似求定积分

的基本思路，就是通过将积分区间[a,b]划分成若干小区间，在每个小区间上用简便易求的函数近似替代被积函数f(x)，并计算每个小区间上的近似函数所围成的面积之和来逼近定积分的值.实验4.3MATLAB数值积分与微分如自适应辛普森(Simpson)法、自适应洛巴托（Lobatto）法、高斯－勒让德(Gauss-Legendre)法、全局自适应求积法等都是经常采用的求数值积分的方法．6MATLAB提供了基于这些算法的相应函数：quad、quadl、quadgk和integral等函数integral和函数quad、quadl、quadgk功能基本相同，但前者更强大、更智能化，主要体现在：实验4.3MATLAB数值积分与微分（1）速度更快；（2）支持积分限为无穷大的积分计算以及含奇点的积分计算（quadgk函数也有此功能）；（3）如果是重积分，integral2和integral3还支持非矩形区域和非长方体区域上的积分.71．基于自适应求积法的MATLAB实现基于自适应求积法，MATLAB给出了integral函数来求定积分.该函数的调用格式为：I=integral(fun,a,b,Name,Value)fun是函数句柄．其中a和b分别是定积分的下限和上限．tol用来控制积分精度，默认值为tol=0.001．Name,Value是用于指定积分选项的名称-值对参数，例如‘AbsTol’和‘RelTol’用于控制绝对和相对误差容限，默认值分别为和.实验4.3MATLAB数值积分与微分I=integral(fun,a,b)或8例10

求定积分解>>I=integral(f10,0,3*pi)↙I=0.9008实验4.3MATLAB数值积分与微分注

integral函数也支持无穷区间，并且能够处理端点包含奇点的情况.9实验4.3MATLAB数值积分与微分例11

求定积分解>>f11=@(x)exp(-x.^2).*(log(x).^2);I=integral(f11,0,inf)↙I=1.9475注本题积分区间端点0为奇点.102．梯形积分法的MATLAB实现

在MATLAB中，对于被积函数以数表的形式给出的定积分问题用trapz函数，调用格式为：其中向量X,Y为等长的两组向量，定义函数关系Y=f(X)．实验4.3MATLAB数值积分与微分I=trapz(X,Y)一般地，积分区间是11实验4.3MATLAB数值积分与微分例12已知某次物理实验测得如下表所示的两组样本点：x1.381.562.213.975.517.799.1911.1213.39y3.353.965.128.9811.4617.6324.4129.8332.21现已知变量x和变量y满足一定的函数关系，但此关系未知，设y=f(x)，求积分的数值.解>>X=[1.38,1.56,2.21,3.97,5.51,7.79,9.19,11.12,13.39];Y=[3.35,3.96,5.12,8.98,11.46,17.63,24.41,29.83,32.21];I=trapz(X,Y)↙I=217.103312实验4.3MATLAB数值积分与微分注函数关系式已知的函数也可以用此命令求定积分的值，需要先生成X，Y的函数关系数据向量，这种函数求得的数值解比函数integral求得的数值精确度低.例13

用trapz函数计算定积分解>>X=1:0.01:2.5;Y=exp(-X.^2);%生成函数关系数据向量trapz(X,Y)↙ans=0.1390133.多重积分数值求解的MATLAB实现使用MATLAB提供的integral2函数和integral3函数可以求出矩形区域上二重积分和长方体区域上三重积分的的数值解．这两个函数的调用格式分别为：I=integral2(fun,xmin,xmax,ymin,ymax)，实验4.3MATLAB数值积分与微分q=integral3(fun,xmin,xmax,ymin,ymax,zmin,zmax)，或I=integral2(fun,xmin,xmax,ymin,ymax,Name,Value)，q=integral3(fun,xmin,xmax,ymin,ymax,zmin,zmax,Name,Value)，14实验4.3MATLAB数值积分与微分Zmin和zmax分别是z变量的积分下限和上限，这些可以是常数，也可以是函数句柄（即可以是x、y的函数）；ymin和ymax分别是y变量的积分下限和上限，这些可以是常数，也可以是函数句柄（即可以是x的函数）；fun是函数句柄；其中xmin和xmax分别是x变量的积分下限和上限，必须是有限或无限的实标量值；Name,Value的用法与函数integral完全相同.15integral2函数和integral3函数不仅可以求出矩形区域上二重积分和长方体区域上三重积分的数值解，也可以求出非矩形区域上二重积分和非长方体区域上三重积分的数值解，还支持含奇点的重积分，当奇异性位于积分边界上时，integral2和integral3的性能最佳.实验4.3MATLAB数值积分与微分16实验4.3MATLAB数值积分与微分例14计算二次积分解>>f14=@(x,y)exp(-x.^2/2).*sin(x.^2+y);I=integral2(f14,-2,2,-1,1)↙I=1.5745注这是函数

在矩形区域[-2,2]×[-1,1]上的二重积分.17实验4.3MATLAB数值积分与微分例15计算二次积分解>>f15=@(x,y)1./(sqrt(x+y).*(1+x+y).^2);ymax=@(x)1-x;%定义y的上限为1-xI=integral2(f15,0,1,0,ymax)↙I=0.2854注这是函数

在三角形区域

上的二重积分.积分边界上含奇点(0,0).18例16计算三次积分解>>f16=@(x,y,z)y.*sin(x)+z.*cos(x);Q=integral3(f16,0,pi,0,1,-1,1)↙Q=

2.0000注这是函数在长方体区域[0,π]×[0,1]×[-1,1]上的三重积分.实验4.3MATLAB数值积分与微分19实验4.3MATLAB数值积分与微分例17计算三次积分解>>f17=@(x,y,z)1./(1+x+y+z);ymax=@(x)(1-x);%定义y的上限为1-xzmax=@(x,y)(1-x-y);%定义z的上限为1-x-yq=integral3(f17,0,1,0,ymax,0,zmax)↙q=0.0966注这是函数在四面体区域

上的三重积分.积分边界上含奇点(0,0,0).20二、数值微分实际中常遇到仅给出了一系列离散点及相应函数值的列表型函数的求导问题，这就需要用这些离散点的函数值推算函数在某点的导数或高阶导数的近似值，这种方法称为数值微分.实验4.3MATLAB数值积分与微分对于难以求导的复杂函数，也可以用数值微分求导，不过需要先由函数表达式生成离散的数据列表.通常用以下三种思路建立数值微分公式：21实验4.3MATLAB数值积分与微分(1)差商近似数值微分：从导数定义出发，通过近似处理，得到数值微分；(2)插值型数值微分：利用本章实验4.1介绍的插值公式得到近似代替该函数的较简单函数，对其求导得到要求导数的近似值.(3)拟合型数值微分：利用本章实验4.2介绍的数据拟合的方法得到近似代替该函数的较简单函数，对其求导得到要求导数的近似值.22实验4.3MATLAB数值积分与微分1．差商近似数值微分（1）数值微分与微商导数定义为假设h>0，引进记号函数f

(x)在x点处以h

为步长的向前差分函数f

(x)在x点处以h

为步长的向前差商当步长h足够小时，有称为向前差商数值微分公式23实验4.3MATLAB数值积分与微分类似可得向后差商数值微分公式中心差商数值微分公式其中中心差商公式精度较高.24DX=diff(X)：计算向量X的向前差分，DX(i)=X(i+1)-X(i)，i=1,2,…,n-1DX=diff(X,n)：计算X的n阶向前差分．例如，diff(X,2)=diff(diff(X))DX=diff(A,n,dim)：计算矩阵A的n阶差分，dim=1时(默认状态)，按列

计算差分；dim=2，按行计算差分．实验4.3MATLAB数值积分与微分（2）差分的MATLAB实现在MATLAB中，没有直接提供求数值导数的函数，只有计算向前差分的函数diff，利用它可以求出微商，从而得到要求的近似导数.diff的调用格式为:25例18根据下表所示年份出生人口数，计算出生人口年增长率解差商近似求数值微分的程序如下：实验4.3MATLAB数值积分与微分年份193019351940194519501955196019651970人口/万650781914100514711861146824792801年份197519801985199019952000200520102015人口/万211418392043262116931379161715741655>>t=[1930:5:2015];p=[650781914100514711861146824792801211418392043262116931379161715741655];26dt=diff(t);%求时间t的差分dp=diff(p);%求人口p的差分q=dp./dt↙%利用差商求数值导数,即出生人口增长率列1至626.200026.600018.200093.200078.0000-78.6000列7至12202.200064.4000-137.4000-55.000040.8000115.6000列13至17-185.6000-62.800047.6000-8.600016.2000差商近似法是最简单的数值微分方法，在实际中十分常用，但其精确度不高,误差较大.实验4.3MATLAB数值积分与微分27实验4.3MATLAB数值积分与微分2．插值型数值微分插值型数值微分是差商近似法的推广.当函数可微性不太好时，利用样条插值进行数值微分要比多项式插值更适宜.仅就三次样条插值方法说明数值微分过程：离散数据三次样条插值函数pp的导数pppp在点xi的导数值fnder是对样条函数求导，fnval用来计算样条函数的函数值.28实验4.3MATLAB数值积分与微分例19某液体冷却时，温度随时间的变化数据如下表所示.试分别计算t=2,3,4min及t=1.5,2.5,4.5min时的降温速率.t/min012345T/oC92.085.379.574.570.267.0分析前者是计算节点处的一阶导数，后者是计算非节点处

的一阶导数.解三次样条插值函数求数值微分的程序如下：>>t=[0:5];T=[92.0,85.3,79.5,74.5,70.2,67.0];29实验4.3MATLAB数值积分与微分p=spline(t,T);%生成三次样条插值函数pp=fnder(p);%生成三次样条插值函数的导函数t1=[2,3,4,1.5,2.5,4.5];dT=fnval(pp,t1);%计算导函数在t1处的导数值disp('相应时间时的降温速率：')disp([t1;dT])↙相应时间时的降温速率：2.00003.00004.00001.50002.50004.5000-5.3722-4.6722-3.8389-5.7972-4.9889-3.2222注插值型数值微分不但适用于求节点处的导数，还可以求非节点处的导数.30实验4.3MATLAB数值积分与微分3．拟合型数值微分如果离散点上的数据有不容忽视的随机误差，应该用曲线拟合代替函数插值，然后用拟合曲线的导数作为所求导数的近似值，这种做法可以起到减少随机误差的作用.仅就多项式拟合方法说明数值微分过程:离散数据多项式拟合函数导函数pppp在点xi的导数值polyder是对多项式求导.31实验4.3MATLAB数值积分与微分例20一底面面积为常数

S的正圆柱体水塔的某一天

0~9

点的水位测量记录(这一时段没有给水塔充水)如下表所示,根据该表估计这一时段任何时刻从水塔流出的水流量.分析由于水塔截面积是常数，为简单起见，计算中将流量定义为单位时间流出水的高度，即水位对时间变化率的绝对值(水位是下降的).时间/h00.921.842.953.874.985.907.017.938.97水位/cm96894893191389888186985283982232实验4.3MATLAB数值积分与微分解方法一多项式拟合求数值微分的程序如下：>>t=[0,0.92,1.84,2.95,3.87,4.98,5.90,7.01,7.93,8.97];h=[968,948,931,913,898,881,869,852,839,822];A=polyfit(t,h,3);%3次多项式拟合B=polyder(A);%对拟合多项式求导tp=0:0.1:9;x=-polyval(B,tp)↙%求tp时刻的水流量x=列1至722.107921.838521.573921.313921.058720.808220.5623列8至1420.321220.084919.853219.626219.404019.186418.973633实验4.3MATLAB数值积分与微分列

15至2118.765518.562118.363418.169517.980217.795717.6158列22至2817.440717.270317.104616.943616.787316.635816.4889列29至3516.346816.209416.076715.948715.825415.706815.5929列36至4215.483815.379415.279615.184615.094315.008814.9279列43至4914.851714.780314.713514.651514.594214.541614.4937列50至5614.450614.412114.378414.349314.325014.305414.290534实验4.3MATLAB数值积分与微分列

57至6314.280314.274814.274114.278014.286714.300114.3182列64至7014.341014.368514.400714.437714.479314.525714.5768列71至7714.632614.693114.758314.828214.902914.982215.0663列78至8415.155115.248515.346715.449715.557315.669615.7867列85至9115.908516.034916.166116.302016.442716.588016.7380我们可以用给定时段的用水量

968-822=146检验计算结果：35实验4.3MATLAB数值积分与微分在上述程序下继续运行>>y=0.1*trapz(x)%用数值积分计算给定时段的总用水量，

积分步长为0.1↙y=146.1815与用水量的绝对误差为146-146.1815=0.1815.36实验4.3MATLAB数值积分与微分方法二差商近似求数值微分的程序如下：>>t=[0,0.92,1.84,2.95,3.87,4.98,5.90,7.01,7.93,8.97];h=[968,948,931,913,898,881,869,852,839,822];dt=diff(t);%求时间t的差分dh=diff(h);%求水位h的差分q=dh./dt↙%利用差分求数值导数，即得已知时刻水流量q=-21.7391-18.4783-16.2162-16.3043-15.3153-13.0435-15.3153-14.1304-16.3462>>u=[0,0.92,1.84,2.95,3.87,4.98,5.90,7.01,7.93];A=polyfit(u,q,3);%用3次多项式拟合流量函数的系数t=0:0.1:9;37实验4.3MATLAB数值积分与微分s=-polyval(A,t)↙%输出t时刻的水流量值s=列1至721.339821.048520.763920.486020.214819.950219.6921列8至1419.440619.195618.957118.724918.499218.279818.0666列15至2117.859717.659117.464617.276217.094016.917816.7476列22至2816.583316.425116.272716.126115.985415.850415.7212列29至3515.597615.479715.367515.260815.159615.063914.973738实验4.3MATLAB数值积分与微分列36至4214.888814.809414.735314.666414.602814.544414.4912列43至4914.443114.400214.362214.329314.301314.278314.2601列50至5614.246814.238314.234614.235514.241214.251514.2665列57至6314.286014.310014.338514.371414.408814.4