2019年数值积分的matlab实现

实验 10数值积分 实验目的 1．了解数值积分的基本原理； 2．熟练掌握数值积分的 MATLAB 实现； 3．会用数值积分方法解决一些实际问题。 实验内容 积分是数学中的一个基本概念， 在实际问题中也有很广泛的应用。 同微分一样， 在 微积分 中， 它也是通过极限定义的， 由于实际问题中遇到的函数一般都以列表形式给出，所以常常不能用来 直接进行积分。此外有些函数虽然有解析式，但其原函数不是初等函数，所sinx 1xd这时我们 一般考虑用数值方法计算其如不定积分。以仍然得不到积分的精确值， 近似值，称为数值积分。10.1 数值微分简介 *  x0 xfxy在设函数可导，则其导数为**xfxhf*limxf）（10.1 ，则其精确值无法求得，但可由h0hyfx以列表形式给出（见表 10-1） fxhfx*xf 10.2）（ 下式求得如果函数其近似值 ** h 表 10-1 x y n hyfx在（10.2）式右端项的分子称为函数一般的，步长越小， 所得结果越精确。 **xx 的差分， 所以右端项又称为差商。的差分，分母称为自变量在数值微分即用差商近似代替微商。常用的差 商公式为 fxhfxh00xf（10.3 ） 02h3y4yy210fx （10.4） y4y3ynn2n1xf（10.5） 三点公式。其误差均为 10.210.2数值微分的数值微分的 0h2． n2h2hO ，称为统称 MATLABMATLAB 实现实现 MATLAB 提供了一个指令求解一阶向前差分，其使用格式为 dxdiffx nx,xxx,x,x1n,这样基于两点的数值导 dx 维数组，为其中 x 是维数组 132n21数可通过指 令 diffx/h 实现。对于三点公式，读者可参考例 1 的 M 函数文件 diff3.m。 xfxxyf的值由下表给用三点公式计算处的导数值，在 1.0,1.2,1.41 1 例 出。 x xf 1.0 1.1 1.2 1.3 1.40.2500 0.2268 0.2066 0.1890 0.1736 解建立三点公式的 M 函数文件 diff3.m 如下 function fdiff3x,y nlengthx;hx2-x1; f1-3*y14*y2-y3/2*h; for j2n-1 fjyj1-yj-1/2*h; end fnyn-2-4*yn-13*yn/2*h; 在 MATLAB 指令窗中输入指令 x[1.0,1.1,1.2,1.3,1.4];y[0.2500,0.2268,0.2066,0.1890,0.1736];diff3x,y yfx所以，-0.0014。 ，运行得各点的导数值为-0.2470，-0.2170-0.1890， -0.1650 x1.0,1.2,1.4 处的导数值分别为-0.2470，-0.1890 和-0.0014 在。 对于高阶导数，MATLAB 提供了几个指令借助于样条函数进行求导，详细使用步骤如下 step1step1对给定数据点（x,y） ，利用指令 ppsplinex,y，获得三次样条函数数据 pp，供后面 ppval 等指令使用。 其中， pp 是一个分段多项式所对应的行向量， 它包含此多项式的阶数、 段数、 节点的横坐标值和各段多项式的系数。 step2step2对于上面所求的数据向量 pp，利用指令[breaks,coefs,m,n]unmkpppp进行处理，生 成几个有序的分段多项式 pp。 step3step3对各个分段多项式 pp 的系数，利用函数 ppval 生成其相应导数分段多项式的系数，再利 用指令 mkpp 生成相应的导数分段多项式 step4step4将待求点 xx 代入此导数多项式，即得样条导数值。 上述过程可建立 M 函数文件 ppd.m 实现如下 function dyppdpp [breaks,coefs,m]unmkpppp; for i1m coefsmi,polydercoefsi,; end dymkppbreaks,coefsm; 于是，如果已知节点处的值 x,y，可用下面指令计算 xx 处的导数 dyy ppsplinex,y,dyppdpp;dyyppvaldy,xx; ysinx的数据点，利用三点公式和三次样条插值分别求导，并与 基于正弦函数例例 2 2 解析所求得 的导数进行比较。 解编写 M 脚本文件 bijiao.m 如下 h0.1*pi;x0h2*pi;ysinx; dy1diff3x,y; ppsplinex,y;dyppdpp;dy2ppvaldy,x; zcosx; error1normdy1-z,error2normdy2-z br--,x,z,plotx,dy1,k,x,dy2,运行得结果为error1 0.0666， error2 0.0025，生成图形见图 10.1。 三点公式、三次样条插值与解析求导比较图图10.1 显然利用三次样条插值求导所得误差比三点公式求导小很多，同时由图 2.15 可知利用三次样条 插值求导所得曲线与解析求导曲线基本重合， 而三点公式在极值点附近和两个端点附近误差较大， 其它点吻合的较好。 10.3应用示例湖水温度变化问题应用示例湖水温度变化问题 问题问题湖水在夏天会出现分层现象，其特点是接近湖面的水的温度较高，越往下水的温度越低。 这种现象会影响水的对流和混合过程，使得下层水域缺氧，导致水生鱼类死亡。对某个湖的水温 进行观测得数据见表 10-2。 表 10-2 某湖的水温观测数据 深度（m） 温度（℃） 0 22.8 2.3 22.8 4.9 22.8 9.1 20.6 13.7 13.9 18.3 11.7 22.9 11.1 27.2 11.1 试找出湖水温度变化最大的深度。 1 1．问题的分析．问题的分析 湖水的温度可视为关于深度的函数，于是湖水温度的变化问题便转化为温度函数的导数问题，显 然导函数的最大绝对值所对应的深度即为温度变化最大的深度。对于给定的数据，但考虑到所给 从而得到温度变化最大的深度，可以利用数值微分计算各深度的温度变化值， 的数据较少，由此计算的深度不够精确，所以采用插值的方法计算加密深度数据的导数值，以得 到更准确的结果。 2．模型的建立及求解 ThT ThTh可（℃） ，相应的温度为，且有记湖水的深度为，并假定函数（m）导。 Th的插值导函数；然后将给定对给定的数据进行三次样条插值，并对其求导，得到的深度数据 加密，搜索加密数据的导数值的绝对值，找出其最大值及其相应的深度，相应的 MATLAB 指令如 下 h