计算方法实验指导书1.doc

第一章 绪论一、主要要求通过实验，认真理解和体会数值计算的稳定性、精确性与步长的关系。 二、主要结果回顾：1、算法：电子计算机实质上只会做加、减、乘、除等算术运算和一些逻辑运算，由这些基本运算及运算顺序规定构成的解题步骤，称为算法它可以用框图、算法语言、数学语言或自然语言来描述。用计算机算法语言描述的算法称为计算机程序。（如c语言程序，c+语言程序，Matlab语言程序等）。2、最有效的算法：应该运算量少，应用范围广，需用存储单元少，逻辑结构简单，便于编写计算机程序，而且计算结果可靠。3、算法的稳定性：一个算法如果输入数据有误差，而在计算过程中舍入误差不增长，则称此算法是数值稳定的，否则称此算法为不稳定的。换句话说：若误差传播是可控制的,则称此算法是数值稳定的,否则称此算法为不稳定的。4、控制误差传播的几个原则：1）防止相近的两数相减；2）防止大数吃小数；3）防止接近零的数做除数；4）要控制舍入误差的累积和传播；5）简化计算步骤,减小运算次数,避免误差积累。三、数值计算实验（以下实验都需利用Matlab软件来完成）实验1.1（体会数值计算精度与步长关系的实验）实验目的：数值计算中误差是不可避免的，要求通过本实验初步认识数值分析中两个重要概念：截断误差和舍入误差，并认真体会误差对计算结果的影响。问题提出：设一元函数f :RR，则f在x0的导数定义为：实验内容:根据不同的步长可设计两种算法,计算f在x0处的导数。计算一阶导数的算法有两种： （1） （2）请给出几个计算高阶导数的近似算法，并完成如下工作：1、对同样的h，比较(1)式和(2)式的计算结果；2、针对计算高阶导数的算法，比较h取不同值时(1)式和(2)式的计算结果。实验要求：选择有代表性的函数f(x)（最好多选择几个），利用Matlab提供的绘图工具画出该函数在某区间的导数曲线f(s)(x)，再将数值计算的结果用Matlab 画出来，认真思考实验的结果。实验分析：不论采用怎样的算法，计算结果通常都会有误差。比如算法(1)，由Taylor公式，知： 所以有 利用上式来计算f (x0)，其截断误差为： 所以误差是存在的，并且当步长h越来越小时，(1)的近似程度也越来越好。类似地可以分析(2)的截断误差为：上述截断误差的分析表明(2)是比(1)更好的算法，因为对步长h(0,n=1,2,当n=1时，得：当n2时，由分部积分可得： n=2,3,另外，还有：实验内容：由递推关系In=1-nIn-1，可得计算积分序列In的两种算法：算法一、直接使用递推公式得： In=1-nIn-1 n=2,3算法二、利用递推公式变形得：实验要求：用上述两种算法分别在计算中采用5位、6位和7位有效数字，请判断哪种算法给出的结果更精确。实验分析：两种算法的优劣可能与你的第一感觉完全不同。设算法一中I1的误差为e1，由I1递推计算In的误差为en算法二中IN的误差为N，由IN向前递推计算In (nN)的误差为n如果假定上述两种算法在后面的计算中都不再引入其他误差，则有：由此可见，算法一中的e1可能很小，但在计算中它的影响急剧扩散，结果真实的数据很快被淹没了；而算法二中的N可能相对比较大，但在计算中误差影响不扩散，某一步产生误差后，该误差对后面的影响不断衰减。注：误差扩散的算法是不稳定的，是我们所不期望的；误差衰减的算法是稳定的，是我们努力寻求的，也是贯穿本课程始终的目标。第二章 一元非线性方程的解法一、主要要求编写二分法、Newton迭代法和快速弦截法通用子程序。二、主要结果回顾1、二分法的基本思想：二分法就是将方程的有根区间对分,然后再选择比原区间缩小一半的有根区间,如此继续下去,直到得到满足精度要求的根为止的一种简单的区间方法。基本法原理： 给定方程f(x)=0，设f(x)在区间 a，b 连续,且f(a)f(b)0，则方程f(x)在(a，b)内至少有一根，为便于讨论,不妨设方程f(x)=0在(a，b)内只有一实根 x *采取使有根区间逐步缩小，从而得到满足精度要求的实根 x * 的近似值 x k 。 取a,b区间二等分的中点，若f(x0)=0,则x0是f(x)=0的实根； 若f(a)f(x0)0 成立,则 x * 必在区间(a, x0)内,取a1=a，b1= x0；否则 x *必在区间(x0,b)内，取a1= x0，b1=b, 这样,得到新区间(a1,b1),其长度为a,b的一半, 如此继续下去，进行k次等分后,得到一系列有根区间：，其中每一个区间长度都是前一个区间长度的一半，从而ak ,bk 的长度为如此继续下去,则有这些区间将收敛于一点，该点即为所求的根.当做到第k步时，有:（为给定的精度）此时即为所求方程的近似解。2、二分法算法：1）输入有根区间的端点a,b及预先给定的精度；2）计算(a+b)/2并将其赋值给变量c，记为(a+b)/2c；3）若f(a)f(c)0，则cb，转向4），否则ca，转向4）；4）若b-a0 输出失败信息 n=1y1*y0b=c y2=ya=c，y1=yb-aen=n+1输出xstopc=(a+b)/2，y=f(c)是是否否是否5、Newton迭代法的思想：给定方程f(x)=0，迭代公式为：，应用该公式来解方程的方法就称为牛顿迭代法。6、Newton迭代法的算法：1）给出初始近似根x0及精度；2）计算x1；3）若|x1-x0|或迭代次数达到预定指定的次数N ，则转向 4）；否则x1x0，转向2）；4）输出满足精度的根x1，结束。（程序框图略）7、收敛性判定定理：定理1：假定函数j (x) 满足下列条件：1、对任意xa,b 有 j(x) a,b 2、存在正数 L1，使对任意x1，x2a,b 有|j(x1) - j(x2)|L|x1-x2| (0L1)则迭代过程 xk+1= j (xk)对于任意初值x0a,b均收敛于方程x= j(x)的根x* ，且有如下的误差估计式：局部收敛性定义： 如果存在不动点x*的某个邻域U(x*,)，使得对于任意初值x0 U(x*,)，迭代公式xk+1= j(xk) （k=0,1,2.)产生的数列xk均收敛与x*，则称迭代公式xk+1= j(xk)是局部收敛的。定理2：设j(x)在x= j(x)的根x*邻近有连续的一阶导数， 且| j(x*)|1，则迭代公式xk+1=j(xk)具有局部收敛性。8、收敛速度定义定义 当k时，有 （C0,且为常数）则称迭代过程是p阶收敛的。特别地，当p=1，0C newton(x0) 即可求出 求方程fc1=0在x0附近的近似解 得出结果ans = 0.5018 3.1407 6.2832 9.4248 图3注：用Newton迭代法可求出多个根。具体做法：先用fplot命令作出函数的图形，再根据图形给出不同的初始值x0，最后通过调用Newton迭代法程序求出非线性方程的所有根。五、思考题思考题1、 在二分法中取区间中点，每次计算一次函数值。如果每次把区间三等份，计算两个函数值，是否可以找出方程的根？给出它的算法。和二分法比较它的效率如何？思考题2、目的：找出一维搜索的最佳途径。内容： 假设f(x)=0在a,b区间内只有一个根（可以是重根），求解该方程等价于在 a,b 区间找|f(x)|的极小值点。设计一种寻找极小值点的方法，使得计算f(x)的次数尽可能的少，并完成数值实验。你能从理论上证明你的搜索方法最好吗？思考题3、目的：了解牛顿收敛法的结构和“局部”收敛性。内容：牛顿法可以直接用来求解复数方程z3-1=0，在复平面上它的三个根是z1*=1 ，。选择中心位于坐标原点边长为2的正方形内的点为初始值，把收敛到三个不同根的初始值分别标上不同的颜色。只要计算足够多的点，你将得到关于牛顿收敛域的彩色图片。第三章 线性方程组的解法一、主要要求编写方程组求解的通用程序：Jacobi迭代、Seidel迭代以及SOR迭代程序 二、主要结果回顾1、迭代法： 它的基本思想是将线性方程组 AX=b 化为 ： X=BX+f 再由此构造向量序列X(k)： X(k+1)=BX (k)+f 若X(k)收敛至某个向量x *,则可得向量X *就是所求方程组 AX=b 的准确解. 线性方程组的迭代法主要有Jocobi迭代法、Seidel迭代法和超松弛(Sor)迭代法.2、收敛性判定定理：定理1：对任意初始向量X(0)及任意向量 f，由此产生的迭代向量序列X(k)收敛的充要条件是。定理2：若|B|1,则迭代公式X(k)=BX(k-1)+f 收敛.3、Jacobi迭代公式：其中：BJ = -D-1(U+L)，fJ=D-1b。4、Seidel迭代公式：其中：Bs= -(D+L)-1U，fs=(D+L)-1b。5、Jacobi迭代、Seidel迭代的算法：Jacobi迭代算法： 1)给出矩阵A，b，x02)计算B= D-1(-U-L), f=D-1b3) y=BX 0+f4)若 |y-x0|=1.0e-6，则转到5），否则，转到3）5）输出y和n。 Seidel迭代算法： 1)给出矩阵A，b，x02)计算B= (D-L)-1U，f=(D-L)-1b，则有：3) y=BX 0+f4)若 |y-x0|=1.0e-6x0=y y=B*x0+fn=n+1输出x结束 图4仿Jacobi迭代法算法和Seidel迭代算法可给出超松驰迭代的算法（略）三、数值计算实验（以下实验都需用Matlab软件来完成）实验3.1（求解线性方程组实验）实验目的：掌握各种迭代法，比较各种迭代法的渐进收敛速度.实验内容：令 1、计算A的条件数cond(A)；（可选用任何一种范数）2、上述方程组是否为病态方程组？若是，如何改变方程组的病态性？ 3、分别用Jacobi迭代、Seidel迭代与超松驰迭代求方程组AX=b的近似解；4、比较Jacobi迭代、Seidel迭代与超松驰迭代的渐进收敛速度。注：上述实验分两次完成。四、具体操作见下例（以下实验都需用Matlab软件来完成）例：用Jacobi方法求解下列方程组，设X(0)=0,精度为10-6解：1）先根据Jacobi算法编写M文件：Jacobi(a,b,X0)2）在Matlab命令窗口输入命令：a=10 -1 0;-1 10 -2;0 -1 10;b=9;7;6;Jacobi(a,b,0;0;0)3）输出结果：y = 0.9938 0.9381 0.6938 n = 9注： n=9为迭代次数。五、思考题思考题：目的：以Hilbert矩阵为例，研究处理病态问题可能遇到的困难。内容： 设Hilbert矩阵为 它是一个对称正定矩阵，而且cond(Hn)随着n的增加迅速增加。其逆矩阵Hn-1=(aij)，这里 1） 画出ln(cond(Hn) n之间的曲线（可以用任何一种范数）。你能猜出ln(cond(Hn) n之间有何种关系吗？提出你的猜测并想法验证。2） 设D是Hn的对角元素开方构成的矩阵。，不难看出它仍然是对称矩阵，而且对角线元素都是1。把Hn变成的技术称为预处理(preconditioning)。画出ln(cond()/ cond(Hn) n之间的曲线（可以用任何一种范数）。对于预处理你能得出什么印象？3） 对于4n12，给出不同的右端项b。分别用X1=Hn-1b，X2=D-1以及X3=Hnb(这是Matlab的一条命令)求解HnX=b，比较计算结果。4） 取不同的n并以Hn的第一列为右端向量b，同高斯塞德尔迭代求解HnX=b，观察其收敛性。最后你能对于有关Hilbert矩阵的计算得出哪些结论。第四章 插值与拟合 一、主要要求编写拉格朗日插值法、分段线性插值法、*三次样条插值通用程序(Matlab自带)。拉格朗日插值多项式：二、主要结果回顾插值法是函数逼近的重要方法之一，有着广泛的应用 。在生产和实验中，函数f(x)或者其表达式复杂不便于计算或者无表达式而只有函数在给定点的函数值(或其导数值) ，此时我们希望建立一个简单的而便于计算的函数j(x)，使其近似的代替f(x)，这就是所谓的插值法.有很多种插值法,其中以拉格朗日(Lagrange)插值和牛顿(Newton)插值为代表的多项式插值最有特点,常用的插值还有Hermit插值,分段插值和样条插值. 1、插值：求近似函数的方法:由实验或测量的方法得到所求函数 y=f(x) 在互异点x0 , x1, . , xn 处的值 y0 , y1 , , yn ,构造一个简单函数 j(x) 作为函数 y=f(x) 的近似表达式 ：y= f(x) j(x)使 j(x0)=y0 , j(x1)=y1 , , j(xn)=yn ,(*)这类问题称为插值问题。 f(x) 称为被插值函数，j(x) 称为插值函数， x0 , x1, . , xn 称为插值节点。(*)式称为插值条件。2、Lagrange插值公式其中x0 , x1,. , xn为插值节点，yj为插值节点xj处的函数值（j=1,2,n）。3、Lagrange插值的截断误差(插值余项)定理:设Ln(x)是过点x0 ，x1 ，x2 ，xn的 n 次插值多项式， f (n)(x)在a，b上连续，f (n+1)(x)在a，b上存在，其中a，b是包含点x0 ，x1 ，x2 ，,xn的任一区间，则对任意给定的xa，b，总存在一点x（a，b）（依赖于x）使其中：。4、拉格朗日插值法程序框图：（见图5）开始输入xi ,yii =0,1,2,.,nL=0,k=0T=1T=T(x-xi)/(xk-xi) i=0,1,2,.,k-1, k+1,.,nL=L+ykTk=nk=k+1输出L结束NY（图5）三、数值计算实验（以下实验都需利用Matlab软件来完成）实验4.1（观察龙格(Runge)现象实验）实验目的：观察拉格朗日插值的龙格(Runge)现象.实验内容：对于函数进行拉格朗日插值，取不同的节点数，在区间-5,5上取等距间隔的节点为插值点，把f(x)和插值多项式的曲线画在同一张图上进行比较。（a可以取任意值）具体步骤：1、 a=1时，1）n=4，作出f(x)和插值多项式的曲线图； 2）n=10，作出f(x)和插值多项式的曲线图；2、a=0.5时，1）n=4，作出f(x)和插值多项式的曲线图； 2）n=10，作出f(x)和插值多项式的曲线图；3、分析上述曲线图，你可以得出什么结论？四、具体操作例1、 给出f(x)=lnx的数值表，用Lagrange插值计算ln(0.54)的近似值。x0.40.50.60.70.8Lnx-0.916291-0.693147-0.510826-0.357765-0.223144解： 1）首先根据程序框图拉格朗日插值法的函数文件： lagrange(x0,y0,x)2）在Matlab命令窗口输入调用命令：x0=0.4: 0.1: 0.8;y0=-0.916291 -0.693147 -0.510826 -0.356675 -0.223144;lagrange(x0,y0,0.54) 3）输出结果为: -0.616143 (精确解: -0.616186)实验4.2 分段插值实验实验目的：分段插值计算，会使用一维插值函数 ,寻找最佳的插值方法。实验内容：在112的11小时内，每隔1小时测量一次温度，测得的温度依次为：5，8，9，15，25，29，31，30，22，25，27，24。1）试估计每隔1/2小时的温度值。2）分别用分段线性插值，立方插值，三次样条插值和最邻近插值估计其值。实验4.3 （插值与拟合实验）实验目的：求下列数据的多项式拟合函数x=0 ,1, 2 ,3, 4 ,5, 6 ,7 ,8, 9 ,10; y=1 ,2.3, 2.1 ,2, 4.6 ,4.7, 4.3, 8.1 ,9.2, 9.8, 10.3;实验内容：1、 做出原始数据的散点图；2、 做出原始数据的折线图（即分段线性插值函数图）；3、做出原始数据的分段二次拟合曲线图（见图6）；4、做出原始数据的分段三次拟合曲线图。（图6）五、思考题思考题1、目的：观察最小二乘多项式的数值不稳定性现象。内容：1、在-1,1区间上取n=20个等距节点，计算出以相应节点上ex的值做为数据样本，2、以1,x,x2,.,xt为基函数做l=3,5,7,9次的最小二乘多项式。3、画出ln(cond(A) n之间的曲线，其中A是确定最小二乘多项式系数的矩阵。4、计算出不同阶最小二乘多项式给出的最小偏差 思考题2、目的：观察对于非光滑函数进行多项式插值的可能性。内容：在0，1上取f(x)=|sinkx|，选择不同的k和n，用等距节点做拉格朗日多插值项式，观察误差大小和收敛情况。第五章 数值积分一、主要要求编写定步长复合梯形公式、定步长复合Simpson公式、变步长梯形法及龙贝格算法*通用子程序；二、主要结果回顾1、梯形公式：对于连续函数f (x)，有积分中值定理其中f ()是被积函数f (x)在积分区间上的平均值。因此，如果我们能给出求平均值f ()的一种近似方法，相应地就可以得到计算定积分的一种数值方法。如果取，即得下列梯形公式： 二、辛普生(Simpson)公式：先用某个简单函数近似逼近f (x)，然后用在a,b区间的积分值近似表示f (x)在a,b区间上的定积分，即取我们可以取为前面介绍的f (x)的插值多项式Ln(x)或拟合多项式P(x)进行近似计算。如取为插值多项式Ln(x)，则相应得到的数值积分公式称为插值型求积公式。如：若考虑过点A,B,C三点的抛物线段L代替曲线段y=f (x)（见图7）ABCaboxy（图7）L就是以a,b,c为节点构造二次插值多项式：取即得下列抛物线（辛普生）公式： 3、误差分析：梯形公式的截断误差(余项)为:辛普生(Simpson)公式的截断误差(余项)：注： 1)抛物线求积公式的计算精度高于梯形求积公式的计算精度。2)当积分区间较大时，利用上述公式计算积分产生的误差也较大。所以在建立积分公式时，应选择小区间和抛物线求积公式。 从余项的讨论看到，积分区间越小，也可使求积公式的截断误差变小。因此，我们经常把积分区间分成若干小区间，在每个小区间上采用次数不高的插值公式，如梯形公式或抛物线公式，构造出相应的求积公式，然后再把它们加起来得到整个区间上的求积公式，这就是复合求积公式的基本思想。 常用的复合求积公式是复合梯形公式和复合辛普生(Simpson)公式4、复合梯形公式：把区间a,b n等分，取节点xi=a+ih，i=0,1,.n, h=(b-a)/n，对每个小区间xi,xi+1用梯形求积公式，再累加起来得： 5、复合辛普生(Simpson)公式：把区间a,b n等分，取节点xi=a+ih，i=0,1,.n, h=(b-a)/n，对每个小区间xi,xi+1用抛物线求积公式，再累加起来得：（其中：为节点xi与xi+1的中点，i=0,1,.n）6、变步长复合积分法基本思想：是在求积过程中，通过对计算结果精度的不断估计，逐步改变步长（逐次分半），直到满足精度要求为止。即按照给定的精度实现步长的自动选取。下面以梯形公式为例来说明这种方法。利用梯形公式可得积分近似值：记为T0，即：此时步长h0=b-a，再将a,b二等分，即取步长使用n=2时的梯形公式可得积分近似值：再将a,b四等分，即取步长使用n=4时的梯形公式可得积分近似值：再将a,b八等分，即取步长使用n=8时的梯形公式可得积分近似值：依次类推可得变步长梯形公式：即： 若预定精度为，以公式(*)计算积分的近似值，直到|Tk+1-Tk|时终止计算,并以当前值Tk+1为所求近似值.对Simpson公式也可类似构造出相应的变步长积分公式。（略）7、算法：a) 复合梯形公式： b) 复合辛普生(Simpson)公式：算法：1）令2）对k=1,2,.,n-1，计算3）。c）变步长梯形公式公式：，这里。算法：1）输入精度，n=1，；2）s=0；3）对k=0,1,2,.,n-1计算：；4），n=2n；5）如果|T2-T1|，则结束；否则执行6）；6），T1=T2，转2）。三、数值计算实验（以下实验都需利用Matlab软件来完成）实验5.1 数值积分实验实验内容：用数值积分法求解下面的积分方程： 1、 用梯形积分公式求解上述积分方程；2、 用抛物线积分公式求解上述积分方程；3、比较不同数值积分公式对方程解的影响。 实验5.2 数值微分实验 实验目的：熟悉MATLAB中实现两点公式的函数diff，称为微分或差分函数。其主要调用格式有如下三种：i）diff（S） 对由findsym返回的自变量，求符号表达式S的微分。ii）diff（S,v）或ddif（S,sym（v） 对自变量v，求表达式S的微分。iii）diff（S,n） 对整数n，对表达式S微分n次。diff（S,v,n）或diff（S,v,n） 这两种格式都可以被识别。如：求sinx2的一阶导数和t6的6阶导数。输入syms x,diff(sin(x2)syms t,diff(t6,6)ans = 2*cos(x2)*x ans = 720x1.21.31.41.51.61.7y=f (x)1.2441.4061.6021.8372.1212.465实验内容：已知函数表试求f (x)在各节点的一、二阶导数的近似值。四、思考题思考题、目的：比较不同数值积分公式对方程解的影响。内容：用n=4，利用复合梯形公式和复合Simpson公式，重复解实验中的积分方程，比较其结果。第六章 微分方程数值解一、主要要求编写Euler法、改进的Euler法、四阶Runge-Kutta法(Matlab自带)通用子程序；二、主要结果回顾解一阶常微分方程初值问题 将区间a,b作n等份，取步长显式Euler公式：隐式Euler公式：改进的Euler公式：或表示为： 三、数值计算实验（以下实验都需利用Matlab软件来完成）实验目的：比较几种数值解法的精度实验内容：用Euler法和Runge-Kutta法求解初值问题： 四、具体操作例：设有微分方程： 1） 求出上述微分方程的解析表达式（精确解）；2） 用Euler法解微分方程（步长取h）。3） 将两方法的结果显示在一张图上，比较其精确度。给出数据：x0 , T, y0 , h输出(显示器)：数值解。如 给出：0 0.4 1 0.1输出 x y0.1 1.100000.2 1.220000.3 1.362000.4 1.52820解：1）由微分方程可得精确解为：2）先编写M文件euler(a,b,x0,y0,h,n) 再根据题目要求调用函数Euler(0.1,0.4,0,1,0.1,4)结果：ans =1.1000 1.2100 1.3310 1.46413）再Matlab命令窗口输入 0.1 0.2 0.3 0.4； 得 1.1103 1.2428 1.3997 1.5836而用Euler法计算得 =1.1050 1.2210 1.3492 1.4909作图：在Matlab命令窗口输入：plot(x,y)hold on plot(x,y1,-r) 从图中可看出，用Euler法计算出的结果与精确值有较大误差，还需改进。五、思考题思考题1、目的：用数值微分的数值解解决实际问题：内容：1. 小型火箭初始质量为900千克，其中包括600千克燃料。火箭竖直向上发射时燃料以15千克/秒的速率燃烧掉，由此产生30000牛顿的恒定推力。当燃料用尽时引擎关闭。设火箭上升的整个过程中，空气阻力与速度平方成正比，比例系数为0.4（千克/米）。重力加速度取9.8米/秒2.A. 建立火箭升空过程的数学模型（微分方程）；B. 求引擎关闭瞬间火箭的高度、速度、加速度，及火箭到达最高点的时间和高度。思考题2、目的：用数值微分的数值解解决实际问题：内容： 种群的数量（为方便起见以下指雌性）因繁殖而增加，因自然死亡和人工捕获而减少。记xk(t)为第t年初k岁（指满k-1岁，未满k岁，下同）的种群数量，bk为k岁种群的繁殖率（1年内每个个体繁殖的数量），dk为k岁种群的死亡率（1年内死亡数量占总量的比例），hk为k岁种群的捕获量（1年内的捕获量）。今设某种群最高年龄为5岁（不妨认为在年初将5岁个体全部捕获），b1=b2=b5=0，b3=2，b4=4，d1=d2=0.3，d3=d4=0.2，h1=400，h2=200，h3=150，h4=100。A. 建立xk(t+1)与xk(t)的关系（k=1,2,5, t=0,1,），如。为简单起见，繁殖量都按年初的种群数量xk(t)计算，不考虑死亡率。B. 用向量表示t年初的种群数量，用bk和dk定义适当的矩阵L，用hk定义适当的向量h，将上述关系表成的形式。C. 设t=0种群各年龄的数量均为1000，求t=1种群各年龄的数量。又问设定的捕获量能持续几年。D. 种群各年龄的数量等于多少，种群数量x(t)才能不随时间t改变。E. 记D的结果为向量x*, 给x* 以小的扰动作为x(0)，观察随着t的增加x(t)是否趋于x*, 分析这个现象的原因。案例分析 估计水塔的水流量该题是1991年美国大学生数学建模竞赛A题（数据略有改动）1、问题的提出某地区用水管理机构需要对居民的用水速度（单位时间的用水量）和日总用水量进行估计。现有一居民区，某自来水是由一个圆柱行水塔提供，水塔高12.2m，塔的 直径为17.4m。水塔是由水泵根据水塔中的水位自动加水，一般水泵每天工作两次。按照设计，当水塔中的水位降至最低水位，约8.2m，水泵自动启动加水；当水位升高到最高水位，约10.8m，水泵停止工作。表3.4给出的是某一天的测量记录数据，测量了28个时刻的数据，当由于水泵正向水塔供水，有3个时刻无法侧到水位（表中为）。表3.4 水塔中水位原始数据时刻(t)/h00.9211.8432.9493.8714.9785.900水位/m9.6779.4799.3089.1258.9828.8148.686时刻(t)/h7.0067.9288.9679.98110.92510.95412.032水位/m8.5258.3888.22010.82010.500时刻(t)/h12.95413.87514.98215.90316.92617.93119.037水位/m10.2109.9369.6539.4099.1808.9218.662时刻(t)/h19.95920.83922.01522.95823.88024.98625.908水位/m8.4338.22010.82010.59110.35410.180试建立数学模型，来估计居民的用水速度和日总用水量。2、 问题的求解1)水塔中水的体积的计算由问题的要求，解决问题的关键是求出居民的用水速度，即单位时间内所用水的体积。假设水塔是标准圆柱体，因此 （*）式中，D为水塔的直径（17.4m）。h为水塔的水位高度。 计算水的流量，首先需要计算出水塔中水的体积，然后再由式（*）可计算出不同时刻水塔中水的体积V。输入t= 0 0.921 1.843 2.949 3.871 4.978 5.900. 7.0067.928 8.967 9.981 10.925 10.954 12.032 12.954 13.87514.98215.90316.82617.93119.037.19.95920.83922.01522.95823.88024.98625.908;h=9.6779.4799.3089.1258.9828.8148.6868.5258.3888.220 10.82010.50010.2109.9369.6539.4099.1808.9218.6628.4338.22010.82010.59110.35410.180;D=17.4;V=pi/4*D2*h;无法得到的水位高度处用代替。得到水塔中水的体积，其中结果列在表3.5中（负值舍去）表3.5 水塔中水的体积时刻(t)/h00.9211.8432.9493.8714.9785.900体积(V)/m32031.12254.02213.32169.82135.82095.92065.4时刻(t)/h7.0067.9288.9679.98110.92510.95412.032体积(V)/m32027.11994.61954.62572.92496.8时刻(t)/h12.95413.87514.98215.90316.82617.93119.037体积(V)/m32427.82362.72295.42237.32182.92121.32059.7时刻(t)/h19.95920.83922.01522.95823.88024.98625.908体积(V)/m32005.31954.62572.92518.42462.02420.72)水塔中水流速度的估计水流速度应该是水塔中水的体积对时间导数，但由于没有水体积的具体数学表达式，只能用均差（差商）近似导数。可以选择的是Matlab提供的函数gradient()。由于在两个阶段水泵向水塔供水的水位无法确定，因此在计算水塔中水流速度时要分三段计算。第一段，从08.967时刻，第二段，10.95420.839时刻，第三段，22.95825.908时刻。输入t1=t(1:10);t2=t(13:23);t3=t(25:28);V1=V(1:10);V2=V(13:23);V3=V(25:28);dV=-gradient(V1,t1) gradient(V2,t2) gradient(V3,t3);得到导数的近似值，并将它们列在表中。见表3.6。表3.6 水塔中水的体积时刻(t)/h00.9211.8432.9493.8714.9785.900流速/（m3h-1)51.120447.609041.507232.224236.447434.689533.8858时刻(t)/h7.0067.9288.9679.98110.92510.95412.032流速/（m3h-1)34.941136.983738.448770.586272.5251时刻(t)/h12.95413.87514.98215.90316.82617.93119.037流速/（m3h-1)72.768365.309461.791860.994257.219055.709557.2190时刻(t)/