《数学实验 第4版》课件 第四章 数值分析

1第四章数值分析

实验4.1插值

实验4.2离散数据的曲线拟合数学实验实验4.3MATLAB数值积分与微分实验4.4常微分方程的数值解2实验4.1插值一、拉格朗日（Lagrange）插值法二、分段线性插值三、三次样条插值四、应用举例3实验4.1插值一、拉格朗日（Lagrange）插值法实际应用中，经常遇到下面的问题：通过实验、测量等方法得到一些离散数据，这些数据从数学的角度可以看作是由某个函数产生的，如何利用这些数据来寻找这个函数满足要求精度，且相对简单的近似表达式呢？我们介绍以下三种方法：多项式插值、分段线性插值、三次样条插值.4这时我们称为插值函数,为插值节点.作为的近似表达式,使满足的一组测量数据设函数,要寻求一个函数1.拉格朗日插值多项式在我们所学的函数类型中，多项式相对比较简单，用多项式作为插值函数是常用的方法，也称为多项式插值法.实验4.1插值5利用上述方法求n次插值多项式计算比较麻烦，由插值多项式的唯一性，无论用什么方法，由n+1个节点确定的n次插值多项式都是相同的.当n+1个节点互不相同时，由线性代数的知识，方程组(3)唯一确定

的一组系数

。由此可见,n+1个节点可以唯一确定一个n次插值多项式.实验4.1插值6比较方便的方法是构造一组基函数：（4）(5)令（6）就是满足条件(1)的n次插值多项式,称为n次Lagrange(拉格朗表示，即日)插值多项式，通常用（7）实验4.1插值72.多项式插值的误差估计关于可以证明有下列结论成立.(8)若

在

上连续，

在

内存在，若

是

上的

个互异节点，则插值多项式

，对任意点

，插值余项为其中

，且与

有关定理(9)实验4.1插值8实验4.1插值functiony=lagrange(x0,y0,x)n=length(x0);m=length(x);fori=1:mz=x(i);s=0.0;fork=1:np=1.0;forj=1:nifj~=kp=p*(z-x0(j))/(x0(k)-x0(j));endends=p*y0(k)+s;endy(i)=s;endy3.拉格朗日插值法的MATLAB实现把这个程序存盘起来，作拉格朗日插值计算时可以随时调用.编写拉格朗日插值法程序：注意:程序中n表示节点的个数，以数组x0,y0输入，m表示插值点的个数，以数组x输入.输出数组y为m个插值.9分别用n=8、10的等距分点进行多项式插值，例1绘制f(x)及插值多项式的图形.解x=[-1:0.25:1];y=1./(1+9*x.^2);x0=[-1:0.05:1];y0=lagrange(x,y,x0);plot(x0,y0,'r--')holdon,y1=1./(1+9*x0.^2);plot(x0,y1),holdoff↙在命令窗口输入：取区间[-1,1]的8等分点作为节点,实验4.1插值函数f(x)及插值多项式的图形如图4.1所示.10实验4.1插值图4.111取区间[-1,1]的10等分点作为节点,类似地可得f(x)及插值多项式,如图4.2所示.实验4.1插值图4.212一般总认为插值多项式的次数越高,逼近的精度越高.但事实上并非如此,由图4.1和图4.2可以看出，离零越远,近似效果越差,以致完全失真.这一现象被称为Runge现象.Runge现象表明,盲目采用提高插值多项式的次数的方法来减少误差是不可取的.因此,我们只有用缩短插值区间,分段进行插值来达到减少误差的目的.与f（x）的近似程度比较好，但当x实验4.1插值13实验4.1插值二、分段线性插值

分段线性插值,就是用连接彼此相邻两节点的直线段形成的折线作为插值函数.MATLAB对分段线性插值提供了插值函数,见下表.函数功能yi=interp1(x,y,xi)对节点(x,y)插值，求xi处的内插值.yi=interp1(x,y,xi,'method')采用指定的方法，对节点(x,y)插值，求xi处的内插值.14实验4.1插值用n=20的等距分点进行线插值，绘制f(x)及插值例2多项式的图形.在命令窗口输入：x=-1:0.1:1;y=1./(1+9*x.^2);xi=-1:0.1:1;yi=interp1(x,y,xi);plot(x,y,'r-',xi,yi,'*')↙解利用plot函数检验插值效果，由图4.3可见，在[-1,1]内线性插值函数收敛于函数

f(x).15

一般地，线性插值函数都有良好的收敛性.数学、物理中用的特殊函数表,计算机绘图都采用了分段线性插值的原理.实验4.1插值图4.316三、三次样条插值

分段线性插值虽然有良好的收敛性,但是由于在节点处不光滑,实用上受到了一定的限制.下面介绍一种在实际应用中,常用的提高光滑度的方法:三次样条插值.(1)在每一个小区间上是三次多项式且则称S(x)为三次样条插值函数.若函数S(x)满足下列条件：给定[a,b]上的n+1个节点满足实验4.1插值17实验4.1插值其中为待定常数.由条件（2）共4n-2个条件，还需要再给两个条件，才能唯一确定一个三次样条插值函数.最常用的是增设边界条件:也称为自然条件.18实验4.1插值在MATLAB中三次样条插值函数为：yi=interp1(x,y,xi,'spline'),

或yi=spline(x,y,xi)其中变量的含义与分段线性插值相同.19例3在一天24小时内,从零点开始每间隔2小时测得的环境温度为(摄氏度)12,9,9,10,18,24,28,27,25,20,18,15,13推测在每1s时的温度，并描绘温度曲线.在命令窗口输入：解t=0:2:24;T=[129910182428272520181513];plot(t,T,'*')ti=0:1/3600:24;T1i=interp1(t,T,ti);holdon,plot(ti,T1i)T2i=interp1(t,T,ti,'spline');holdon,plot(ti,T2i,'r-')↙实验4.1插值20输入不同的时刻ti便可得到相应的温度值Ti.利用线性插值方法描绘的温度曲线，如图4.4所示.通过图上两条曲线比较可以看出,三次样条插值函数在整个区间上有较好的收敛性;光滑性比分段线性插值有较大的提高,因此应用比较广泛.缺点是误差估计比较困难.这正像我们的人生在发扬优点修正错误中不断砥砺前行.图4.4实验4.1插值21四、应用举例数据加细问题在机械制造加工中,经常遇到数据加细的问题.例如，在现代机械工业中，用计算机程序控制加工机器零件时，根据设计可以给出零件外形曲线上某些点，加工时为控制每步刀的走向及步长就要算出零件外形曲线上其他点的函数值，加工出外表光滑的零件，这就涉及到数据加细的问题.例4

在飞机的机翼加工时，由于机翼尺寸很大，通常在图纸上只能标出部分关键点的数据.某型号飞机的机翼上缘轮廓线的部分数据如下：实验4.1插值x0.004.749.0519.0038.0057.0076.0095.00114.00133.00152.00171.00190.00y0.005.238.1011.9716.1517.1016.3414.6312.166.697.033.990.00对表中的数据进行细化，并画出机翼的上轮廓线.22解输入：x=[0.004.749.0519.0038.0057.0076.0095.00114.00133.00152.00171.00190.00];y=[0.005.238.1011.9716.1517.1016.3414.6312.169.697.033.990.00];xi=[0:0.001:190];yi=interp1(x,y,xi,'spline');plot(xi,yi)↙实验4.1插值机翼的上轮廓线，如图所示.图4.523实验4.1插值例5

天文学家在1914年8月份的7次观测中，测得地球与金星之间距离（单位：米），并取其常用对数值与日期的一组历史数据如下所示,试推断何时金星与地球的距离(米)的对数值为9.9352.日期:18202224262830距离对数:9.96189.95449.94689.93919.93129.92329.9150解由于对数值9.9352位于24和26两天所对应的对数值之间,所以对上述数据用三次样条插值加细为步长为1的数据,输入:x=[18:2:30];y=[9.96189.95449.94689.93919.93129.92329.9150];xi=[18:1:30];yi=interp1(x,y,xi,'spline')；A=[xi;yi]↙24经三次样条插值推断，25日时金星与地球的距离(米)的对数值为9.9352.A=18.000019.000020.000021.000022.000023.000024.00009.96189.95819.95449.95069.94689.94309.939125.000026.000027.000028.000029.000030.00009.93529.93129.92729.92329.91919.9150实验4.1插值25第四章数值分析

实验4.1插值

实验4.2离散数据的曲线拟合数学实验实验4.3MATLAB数值积分与微分实验4.4常微分方程的数值解26实验4.2离散数据的曲线拟合一、离散数据的多项式拟合二、曲线拟合的线性最小二乘法三、应用举例数学实验27实验4.2离散数据的曲线拟合一、离散数据的多项式拟合p=polyfit(x,y,n)用多项式拟合一组离散数据就是寻找一组多项式的系数使得多项式能够较好的拟合这组数据.它与实验4.1的插值法不同,数据不能保证都在拟合多项式曲线上，但能使整体拟合误差较小.在MATLAB中,多项式拟合可以通过polyfit函数来实现,该函数的调用格式为p是多项式系数按降幂排列得出的行向量.28利用poly2sym(p)可以得出相应多项式的表达式.如果要计算拟合多项式在x处的值y,输入y=polyval(p,x),就可输出y的值.求本章例4中，数据组的3次、6次和8次多项式并作图例6解输入：x=[0.004.749.0519.0038.0057.0076.0095.00114.00133.00152.00171.00190.00];y=[0.005.238.1011.9716.1517.1016.3414.6312.166.697.033.990.00];p3=polyfit(x,y,3);y1=polyval(p3,x);p6=polyfit(x,y,6);y2=polyval(p6,x);plot(x,y,'*',x,y1,'--',x,y2,'-.')↙实验4.2离散数据的曲线拟合29从图4.6可见，6次拟合多项式比3次多项式的拟合效果好.实验4.2离散数据的曲线拟合30在命令窗口输入：P6=polyfit(x,y,6);y2=polyval(p6,x);P8=polyfit(x,y,8);y3=polyval(p8,x);plot(x,y,'*'x,y2,'—.',x,y3)↙实验4.2离散数据的曲线拟合由图4.7所示，8次多项式比6次的拟合效果更好.可见，随着多项式次数的不断增加，拟合的效果也越来越好，当拟合多项式的次数就能得出较好的效果.那么利用多项式进行拟合时，是否多项式的次数越高拟合效果一定就越好呢？31设已知数据来自函数8次多项式拟合，并作图.例7解试用生成的数据进行3次、6次和在命令窗口输入：x0=-1:.01:1;y0=1./(1+9*x0.^2);p3=polyfit(x0,y0,3);y1=polyval(p3,x0);p6=polyfit(x0,y0,6);y2=polyval(p6,x0);p8=polyfit(x0,y0,8);y3=polyval(p8,x0);plot(x0,y0,'*',x0,y1,'--',x0,y2,'-.',x0,y3)↙32由图4.8可以看出,多项式拟合的效果并不一定总是很精确的.下面我们来介绍另一种方法—曲线拟合的线性最小二乘法.实验4.2离散数据的曲线拟合图4.833二、曲线拟合的线性最小二乘法最小,就是曲线拟合得最好.曲线拟合常用的方法是线性最小二乘法.

曲线拟合（curvefitting）是指选择适当的曲线类型来拟合观测数据，并用拟合的曲线方程分析两变量间的关系.实验4.2离散数据的曲线拟合已知某函数的一组测量数据根据这组数据寻求曲线逼近曲线因为测量时可能产生误差,所以我们不要求都经过这些点,只要与的距离最为接近,即34一般情况下,可以假设数据拟合曲线为将（2）代入(1),上述问题转化为根据多元函数极值的必要条件实验4.2离散数据的曲线拟合35可以证明方程组(4)的系数矩阵是可逆的,则方程组(4)有唯一解于是有这种求拟合曲线的方法称为曲线拟合的线性最小二乘法.实验4.2离散数据的曲线拟合36在MATLAB优化工具箱中，提供了函数lsqcurvefit,求解最小二乘曲线拟合问题.该函数的调用格式为a=lsqcurvefit(fun,a0,x,y)其中，fun是自定义函数的MATLAB表示，可以用inline('函数内容'，自变量列表)；a0是a的初始预测值，使用方法见例8.实验4.2离散数据的曲线拟合37三、应用举例用切削机床进行金属品加工时，为了适当的调整机床，需要测定刀具的磨损速度.在一定的时间测量刀具的厚度，得数据如下表所示：例8切削时间t/h012345678刀具厚度y/cm30.029.1

28.428.1

28.0

27.727.5

27.227.0切削时间t/h910111213141516刀具厚度y/cm26.8

26.5

26.326.125.725.324.824.0假设经验公式是试用最小二乘法确定实验4.2离散数据的曲线拟合38解定义函数f，在命令窗口输入:f=inline('a(1).*t.^3+a(2).*t.^2+a(3).*t+a(4)','a','t')确定经验公式的系数并作图，在命令窗口输入:t=[0:1:16];y=[30.029.128.428.128.027.727.527.227.026.826.526.326.125.725.324.824.0];a0=[0,0,0,0]a=lsqcurvefit(f,a0,t,y)y1=a(1).*t.^3+a(2).*t.^2+a(3).*t+a(4);plot(t,y,'*',t,y1)↙实验4.2离散数据的曲线拟合a=-0.00290.0678-0.713329.824939系数散点及拟合曲线图，实验4.2离散数据的曲线拟合如图4.9所示.40例9在加压实验中的应力-应变关系测试点的数据如下表所示：已知应力-应变关系可以用一条指数曲线来描述.即假设使用最小二乘法确定参数实验4.2离散数据的曲线拟合41解选取指数函数作拟合时,在拟合前需作变量代换,化为将(6)变形后取对数，得令式（7）化为在命令窗口输入：x=[500*1.0e-61000*1.0e-61500*1.0e-62000*1.0e-62375*1.0e-6];y=[3.100*1.0e+32.470*1.0e+31.953*1.0e+31.515*1.0e+31.217*1.0e+3];实验4.2离散数据的曲线拟合42w=[1.552.472.933.032.89];plot(x,w,'*')↙holdon,y1=exp(8.3020)*x.*exp(-495.4888*x);plot(x,y1,'r-'),holdoff↙实验4.2离散数据的曲线拟合z=log(y)↙z=8.03927.81207.57717.32327.1041a=polyfit(x,z,1)↙a=-495.48888.3020k1=exp(8.2959)↙k1=4.0319e+03所求的拟合曲线(见图4.10)所以43指数函数经常应用在预应力混凝土梁的分析中，作为应力－应变关系的数学模型.图4-10实验4.2离散数据的曲线拟合44在实际应用中常见的拟合曲线有直线多项式(一般n=2,3,不宜过高)双曲线（一支）指数曲线实验4.2离散数据的曲线拟合第四章数值分析

实验4.1插值

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

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

的基本思路，就是通过将积分区间[a,b]划分成若干小区间，在每个小区间上用简便易求的函数近似替代被积函数f(x)，并计算每个小区间上的近似函数所围成的面积之和来逼近定积分的值.实验4.3MATLAB数值积分与微分如自适应辛普森(Simpson)法、自适应洛巴托（Lobatto）法、高斯－勒让德(Gauss-Legendre)法、全局自适应求积法等都是经常采用的求数值积分的方法．50MATLAB提供了基于这些算法的相应函数：quad、quadl、quadgk和integral等函数integral和函数quad、quadl、quadgk功能基本相同，但前者更强大、更智能化，主要体现在：实验4.3MATLAB数值积分与微分（1）速度更快；（2）支持积分限为无穷大的积分计算以及含奇点的积分计算（quadgk函数也有此功能）；（3）如果是重积分，integral2和integral3还支持非矩形区域和非长方体区域上的积分.511．基于自适应求积法的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)或52例10

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

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

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

在MATLAB中，对于被积函数以数表的形式给出的定积分问题用trapz函数，调用格式为：其中向量X,Y为等长的两组向量，定义函数关系Y=f(X)．实验4.3MATLAB数值积分与微分I=trapz(X,Y)一般地，积分区间是55实验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.103356实验4.3MATLAB数值积分与微分注函数关系式已知的函数也可以用此命令求定积分的值，需要先生成X，Y的函数关系数据向量，这种函数求得的数值解比函数integral求得的数值精确度低.例13

用trapz函数计算定积分解>>X=1:0.01:2.5;Y=exp(-X.^2);%生成函数关系数据向量trapz(X,Y)↙ans=0.1390573.多重积分数值求解的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)，58实验4.3MATLAB数值积分与微分Zmin和zmax分别是z变量的积分下限和上限，这些可以是常数，也可以是函数句柄（即可以是x、y的函数）；ymin和ymax分别是y变量的积分下限和上限，这些可以是常数，也可以是函数句柄（即可以是x的函数）；fun是函数句柄；其中xmin和xmax分别是x变量的积分下限和上限，必须是有限或无限的实标量值；Name,Value的用法与函数integral完全相同.59integral2函数和integral3函数不仅可以求出矩形区域上二重积分和长方体区域上三重积分的数值解，也可以求出非矩形区域上二重积分和非长方体区域上三重积分的数值解，还支持含奇点的重积分，当奇异性位于积分边界上时，integral2和integral3的性能最佳.实验4.3MATLAB数值积分与微分60实验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]上的二重积分.61实验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).62例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数值积分与微分63实验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).64二、数值微分实际中常遇到仅给出了一系列离散点及相应函数值的列表型函数的求导问题，这就需要用这些离散点的函数值推算函数在某点的导数或高阶导数的近似值，这种方法称为数值微分.实验4.3MATLAB数值积分与微分对于难以求导的复杂函数，也可以用数值微分求导，不过需要先由函数表达式生成离散的数据列表.通常用以下三种思路建立数值微分公式：65实验4.3MATLAB数值积分与微分(1)差商近似数值微分：从导数定义出发，通过近似处理，得到数值微分；(2)插值型数值微分：利用本章实验4.1介绍的插值公式得到近似代替该函数的较简单函数，对其求导得到要求导数的近似值.(3)拟合型数值微分：利用本章实验4.2介绍的数据拟合的方法得到近似代替该函数的较简单函数，对其求导得到要求导数的近似值.66实验4.3MATLAB数值积分与微分1．差商近似数值微分（1）数值微分与微商导数定义为假设h>0，引进记号函数f

(x)在x点处以h

为步长的向前差分函数f

(x)在x点处以h

为步长的向前差商当步长h足够小时，有称为向前差商数值微分公式67实验4.3MATLAB数值积分与微分类似可得向后差商数值微分公式中心差商数值微分公式其中中心差商公式精度较高.68DX=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的调用格式为:69例18根据下表所示年份出生人口数，计算出生人口年增长率解差商近似求数值微分的程序如下：实验4.3MATLAB数值积分与微分年份193019351940194519501955196019651970人口/万650781914100514711861146824792801年份197519801985199019952000200520102015人口/万211418392043262116931379161715741655>>t=[1930:5:2015];p=[650781914100514711861146824792801211418392043262116931379161715741655];70dt=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数值积分与微分71实验4.3MATLAB数值积分与微分2．插值型数值微分插值型数值微分是差商近似法的推广.当函数可微性不太好时，利用样条插值进行数值微分要比多项式插值更适宜.仅就三次样条插值方法说明数值微分过程：离散数据三次样条插值函数pp的导数pppp在点xi的导数值fnder是对样条函数求导，fnval用来计算样条函数的函数值.72实验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];73实验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注插值型数值微分不但适用于求节点处的导数，还可以求非节点处的导数.74实验4.3MATLAB数值积分与微分3．拟合型数值微分如果离散点上的数据有不容忽视的随机误差，应该用曲线拟合代替函数插值，然后用拟合曲线的导数作为所求导数的近似值，这种做法可以起到减少随机误差的作用.仅就多项式拟合方法说明数值微分过程:离散数据多项式拟合函数导函数pppp在点xi的导数值polyder是对多项式求导.75实验4.3MATLAB数值积分与微分例20一底面面积为常数S的正圆柱体水塔的某一天

0~9

点的水位测量记录(这一时段没有给水塔充水)如下表所示,根据该表估计这一时段任何时刻从水塔流出的水流量.分析由于水塔截面积是常数，为简单起见，计算中将流量定义为单位时间流出水的高度，即水位对时间变化率的绝对值(水位是下降的).时间/h00.921.842.953.874.985.907.017.938.97水位/cm96894893191389888186985283982276实验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.973677实验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.290578实验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检验计算结果：79实验4.3MATLAB数值积分与微分在上述程序下继续运行>>y=0.1*trapz(x)%用数值积分计算给定时段的总用水量，

积分步长为0.1↙y=146.1815与用水量的绝对误差为146-146.1815=0.1815.80实验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;81实验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.973782实验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.450514.4966列64至7014.546914.601514.660214.723214.790214.861414.9366列71至7715.015815.098915.186015.277015.371815.470415.572783实验4.3MATLAB数值积分与微分列78至8415.678815.788615.902016.019016.139616.263716.3913列85至9116.522316.656716.794416.935517.079917.227517.3783我们仍然可以检验结果：w=0.1*trapz(s)%用数值积分计算给定时段的总用水量，

积分步长为0.1↙w=144.0565与用水量的绝对误差为146-144.0565=1.9435.84实验4.3MATLAB数值积分与微分注比较两种方法发现，拟合法误差更小，更精确.拟合法的不足之处在于：当给出等距节点时，不如差商近似法通用性强，拟合的多项式对这一组数据适用，对另一组数据可能就要用另外的多项式.因此，在实际问题中使用哪种方法要视具体问题而定，有时几种方法都要用到.实验4.3MATLAB数值积分与微分小结积分描述了一个函数的整体或者宏观性质，对函数的形状在小范围的改变不敏感，并且积分过程是对数据点进行求和，正的和负的随机误差倾向于相互抵消，因此，一般数值积分过程是稳定的，所得的解精确度也较高.微分则描述了一个函数在一点处的斜率，是函数的微观性质，它很敏感，一个函数小的变化，容易产生相邻点的斜率的大的改变.8586实验4.3MATLAB数值积分与微分差商近似就是用给定函数f(x)曲线上的割线斜率近似切线斜率，插值型数值微分和拟合型数值微分都是用近似多项式的曲线斜率近似给定函数f(x)的曲线斜率.无论是割线斜率，还是近似多项式的曲线斜率，都可能和给定函数f(x)的曲线斜率有很大不同，特别是当f(x)在给定区间内变化比较大时更是这样.同时，微分过程是对数据点进行相减，正的和负的随机误差倾向于相加，这都使得数值微分的解不稳定，并且精度也较差.第四章数值分析

实验4.1插值

实验4.2离散数据的曲线拟合数学实验实验4.3MATLAB数值积分与微分实验4.4常微分方程的数值解88实验4.4常微分方程的数值解一、几种求常微分方程数值解的方法二、应用举例实验4.4常微分方程的数值解89实验4.4常微分方程的数值解实验目的通过本实验了解常微分方程的数值解的概念，掌握利用MATLAB求常微分方程的数值解的方法.一、几种求常微分方程数值解的方法常微分方程是研究函数变化规律的有力工具，但是绝大多数变系数方程、非线性方程都是所谓“解不出来”的，于是常微分方程的数值解法就成为解常微分方程的主要手段．本实验只考虑初值问题.90实验4.4常微分方程的数值解常微分方程初值问题：保证方程的解存在且惟一在一系列离散点上，求的近似值，通常取等步长h，即因此数值解法得到的近似解是一个离散的函数表.91实验4.4常微分方程的数值解一、几种求常微分方程数值解的方法1.欧拉（Euler）方法欧拉方法是一种最古老、最简单而直观的解微分方程的数值方法，其基本想法是在小区间上用差商代替方程左端的导数,而方程右端已知函数中的x在小区间上的哪一点取值，则有以下不同的方法：92实验4.4常微分方程的数值解（1）向前欧拉公式中的x取小区间的左端点xn，以分别代替，就得到：向前欧拉公式（2）向后欧拉公式中的x取小区间的右端点xn+1，就得到：向后欧拉公式上式右端的yn+1未知，故称为隐式公式，无法用它直接计算yn+1．93实验4.4常微分方程的数值解（3）梯形公式将向前欧拉公式和向后欧拉公式加以平均，得到梯形公式94实验4.4常微分方程的数值解（4）改进的欧拉公式先由向前欧拉公式算出yn+1的预测值，再把它代入梯形公式右端，作为校正，即改进的欧拉公式95实验4.4常微分方程的数值解它还可写作上面4个公式中，我们常用的是便于计算的向前欧拉公式和改进的欧拉公式.96实验4.4常微分方程的数值解2.龙格—库塔（Runge-Kutta）公式在区间内多取几个点，就可以构造出精度更高的计算公式，这就是龙格—库塔公式，它是一类方法的总称.它是欧拉方法的一种推广，也是应用最广的求解常微分方程数值问题的方法.一般的龙格—库塔方法的形式为p阶龙格—库塔方法其中ai,bij,ci为待定参数97实验4.4常微分方程的数值解（1）二阶龙格—库塔公式其中为待定系数.满足：上式有4个未知数而只有3个方程，所以解不唯一，存在无穷多个解，可见二阶龙格—库塔方法是一族公式.

改进的欧拉公式98实验4.4常微分方程的数值解（2）四阶经典（标准）龙格—库塔公式实际应用中，用的最多的是四阶龙格—库塔公式四阶龙格—库塔公式也不止一个99实验4.4常微分方程的数值解3.龙格—库塔方法的MATLAB实现在MATLAB中，求微分方程数值解的函数调用格式为[t,x]=solver（’f’,ts,x0,options）自变量值函数值ode45、ode23、ode113、ode15s、ode23s中的一个由待解方程写成的m-文件名=[t0,tf],t0、tf为自变量的初值和终值函数的初值用于设定误差限options=odeset（'reltol',rt,'abstol',at）默认相对误差为10-3，绝对误差为10-6.相对误差绝对误差100实验4.4常微分方程的数值解在诸多的solver函数中，ode23、ode45都是基于龙格—库塔公式的函数:ode23是采用的2、3阶龙格-库塔组合算法，ode45是采用的4、5阶龙格-库塔组合算法.ode45比ode23精度更高，尤其成为大部分场合的首选算法.注(1)在解n个未知函数的方程组时，x0和x均为n维向量，m-文件中的待解方程组应以x的分量形式写成.(2)使用Matlab软件求数值解时，高阶微分方程必须等价地变换成一阶微分方程组.101例21

求微分方程的数值解，解先编写函数文件：functiondx=fun23(t,x)dx=t+x;再求方程的解实验4.4常微分方程的数值解已知精确解为>>ts=0:0.1:1;x0=1;[t,x]=ode45('fun23',ts,x0);y=2*exp(t)-t-1;[t,x,y]plot(t,x,'r-',t,y,'b-.'),gridon;title('SolutionofExample3');xlabel('timet');ylabel('solutionx');legend('x','y');↙ans=01.00001.0000102实验4.4常微分方程的数值解0.10001.11031.11030.20001.24281.24280.30001.39971.39970.40001.58361.58360.50001.79741.79740.60002.04422.04420.70002.32752.32750.80002.65112.65110.90003.01923.01921.00003.43663.4366

t的值x的值y的值例22

求微分方程组的数值解.实验4.4常微分方程的数值解103解先编写函数文件：functiondx=fu