22 一些数值方法的建模问题.docx

1、22 一些数值方法的建模问题2.2 一些数值方法的建模问题2.2.1. 排水渠道的设计问题1) 问题的提出 土木工程师和环境工程师在设计一条排水渠道时必须考虑渠道的宽度、深度、内壁光滑度等参数及水流的速度、流量、水深等物理量之间的关系。现在要设计一条横断面为矩形的水渠，其宽度为20米,水流量为5立方米/秒，水渠的斜度系数为0.0002,Manning粗糙系数为0.03，试确定渠道中水的深度。2） 符号约定H：渠道中水的深度，单位为米B：渠道宽度为，单位为米，本题中B=20米Q：水流量，单位为立方米/秒，本题中Q=5立方米/秒S: 水渠的斜度系数,无量刚量, 本题中S=0.0002N：Manni

2、ng粗糙系数,是一个和水渠内壁材料的光滑度有关的无量刚量,本题中N =0.03U：水的流速，单位为米/秒3） 问题分析和建立模型 查阅水工学使用的有关资料，获得一个有关流速的关系式： （1）而质量守衡定律告诉我们Q=UBH (2)将式（1）代入式（2），同时将S,N,B,Q用给定的数值替代，有 （3）令，于是问题归结为求方程的正根问题。4） 模型求解由于方程没有求根公式，故采用数值方法求根，这里采用Newton迭代法来求根。为找到迭代初始值x0，先画出的图形观察：In1:= fh_=Sqrt0.0002/0.03*(20*h)(5/3)*(20+2h)(-2/3)-5; In2:= Plotf

3、h,h,0,5输出图形为Out2= -Graphics-从图中可以看到函数在1附近有根，选取迭代初始值x0=1，根的精度为10-6做Newton迭代，对应的程序为Clearx,f,g;fx_=Sqrt0.0002/0.03*(20*x)(5/3)*(20+2x)(-2/3)-5;gx_=Dfx,x;x0=1；eps1=0.1;nmax=500;eps=0.000001;Dou1=gx0; IfAbsu1/Neps1,Print迭代法失效;Break; x=Nx0-fx0/u1,10; u1=Absx-x0/N; PrintH=,x, n=,n, eps=,u1; Ifu1eps,Print迭代

4、失败 执行程序后得输出结果为H=0.7292262368 n=1 eps=0.270774H=0.7025912396 n=2 eps=0.026635H=0.702293294 n=3 eps=0.000297946H=0.7022932563 n=4 eps=3.777210 -8从输出结果说明所求的水渠深度约为0.7022932563，其误差为3.777210 -82.2.2男大学生的身高问题1） 问题的提出 有关统计资料表明,我国大学生男性群体的平均身高约为170cm，且该群体中约有99.7的人身高在150cm至190cm之间。如果将150cm，190cm等分成20个高度区间，试问该群

5、体身高在每一高度区间的分布情况怎样？特别地，身高中等（165cm至175cm之间）的人占该群体的百分比会超过60吗？解：2）问题分析和建立模型 因为一个人的身高涉及很多因素，通常它是一个服从正态分布N(,)的随机变量。正态分布的概率密度函数(x)为：于是密度函数(x)在区间a,b上的定积分的值正好代表大学生身高在区间acm,bcm的分布。在密度函数(x)中的两个参数、分别为正态分布的均值和标准差。根据题目中我国大学生男性群体的平均身高约为170cm，可选取正态分布的均值=170cm，而由“该群体中约有99.7的人身高在150cm至190cm之间”和正态分布N(,)的“3规则”，有-3=150c

6、m+3=190cm于是可以得到=,故其密度函数(x)为将150cm，190cm等分成20个高度区间后，得到高度区间为150，152，152，154，188，190对应的分布为 （1）身高在165cm至175cm之间的人占该群体的百分比为 （2）如上式（1）和（2）的定积分是不能用定积分基本公式方法求出的，但用计算方法中的数值积分可以算出。3） 模型求解选用数值积分中的复化梯形公式求积方法编程可以计算出求误差小于10-4的式（1）的定积分值从而得出相应分布依次为150cm，152cm的分布为0.0021152cm，154cm的分布为0.0047154cm，156cm的分布为0.0097156cm

7、，158cm的分布为0.0181158cm，160cm的分布为0.0309160cm，162cm的分布为0.0483162cm，164cm的分布为0.0690164cm，166cm的分布为0.0902166cm，168cm的分布为0.1078168cm，170cm的分布为0.1179170cm，172cm的分布为0.1179172cm，174cm的分布为0.1078174cm，176cm的分布为0.0902176cm，178cm的分布为0.0690178cm，180cm的分布为0.0483180cm，182cm的分布为0.0309182cm，184cm的分布为0.0181184cm，186cm

8、的分布为0.0097186cm，188cm的分布为0.0047188cm，190cm的分布为0.0021对（2）式有0.5467246173说明身高中等（165cm至175cm）的大学生约为54.67，不足60。4）Mathematica编写的通用程序为Cleara,b,x,n,f;a=Inputa=b=Inputb=fx_=3/20*Exp-9(x-170)2/800/Sqrt2Pieps=0.0000001;n=1;h=b-a;t1=(fa+fb)*h/2;h=h/2;t2=t1/2+h*fa+h;er=t2-t1/N;WhileAbsereps,Printn=,2n,定积分值为,Nt2,

9、10;Print误差=,er;h=h/2;t1=t2;n=n+1;t2=t1/2+h*Sumfa+k*h,k,1,2n,2;er=t2-t1/N;Printn=,2n,定积分值为,Nt2,10;Print误差=,er说明 本程序从梯形公式T1开始，用复合梯形求积公式求a,b上满足精度小于eps要求的定积分近似值。程序执行后，按要求通过键盘输入积分下限a、积分上限b、被积函数f(x)和精度要求eps 后，计算机则给出满足精度要求的定积分近似值及中间计算值和误差。程序中变量说明a：存放积分下限b: 存放积分上限fx: 存放被积函数f(x)eps: 存放求积精度h: 存放节点步长x：为函数fx:提供

10、变量t1: 存放复合梯形值Tnt2: 存放复合梯形值T2n和定积分近似值n: 存放复合梯形公式的节点加密次数er:存放误差注：为提高计算效率，计算T2n时使用了公式2.2.3湖水温度变化问题1) 问题的提出 湖水在夏天会出现分层现象，其特点为接近湖面的水温度较高，越往下温度变低。这种上热下冷的现象影响了水的对流和混合过程，使得下层水域缺氧，导致水生鱼类的死亡。下面是某个湖的观测数据深度（m）02.34.99.113.718.322.927.2温度(C)22.822.822.820.613.911.711.111.1请问1 湖水在10m处的温度是多少？2 湖水在什么深度温度变化最大？2）问题的分

11、析和假设本问题只给出了有限的实验数据点，可以想到用插值和拟合的方法来解决题目的要求。假设湖水深度是温度的连续函数，引入符号如下：h: 湖水深度,单位为mT: 湖水温度,单位为C,它是湖水深度的函数：T=T(h)这里多项式拟合的方法来求出湖水温度函数T(h)。然后利用求出的拟合函数就可以解决本问题了。3）模型的建立将所给数据作图，横轴代表湖水深度，纵轴代表湖水温度，用Mathematica 数学软件画出散点图，操作的命令为：In1:= d=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.1In2:=

12、q=ListPlotd,PlotStyle-PointSize0.04得出如下散点图Out2=-Graphics-观察散点图的特点，并通过实验选取不同的基函数类进行实验发现用4次多项式拟合比较好，相应的命令为：In3:= q1=Fitd,1,h,h2,h3,h4, h5,hOut3=22.711+0.0280756 h+0.0865552 h 2 -0.0235702 h 3 +0.00132018 h 3 -0.0000218133 h 5In4:= q2=Plotq1,h,0,27.2Out4= -Graphics-In5:= Showq,q2Out5= -Graphics-因此有T(h)

13、 q1=22.711+0.0280756 h+0.0865552 h 2 -0.0235702 h 3 +0.00132018 h 3 -0.0000218133 h 5为求出湖水在10m处的温度，键入命令：In6:= q1/.h-10Out6= 19.0975于是得出湖水在10m处的温度约为19.0975 C。为求湖水在什么深度温度变化最大，要求出函数T(h)的导函数T(h)的绝对值最大值点。为此对所求拟合函数q1关于h的导数并找出最大值点，键入如下命令：In7:= q2=Dq1,h Out7= 0.0280756+0.17311h-0.0707107h 2 + 0.00528073 h 3

14、 - 0.000109067 h 4In8:= Plotq2,x,0,27.2 输出图形为 Out8= -Graphics-从导函数图形上观察到其在10附近可以取得绝对值最大值，键入命令：In9:= FindMinimumq2,x,10 Out9= -1.21357, h - 11.9312得出 h=11.9312 是导函数的绝对值最大值点，于是可以知道湖水在深度为11.9312m 时温度变化最大。2.2.4储量计算问题1) 问题的提出 某地区为估计某矿物的储量，在该地区内进行勘探，得到如下数据（见表一）。请估计出此地区内(1X4，1Y5)该矿物的储量编号12345678910x坐标(KM)1

15、111122222y坐标(KM)1234512345矿物体厚度H(M)13.7225.808.4725.2722.3215.4721.3314.4924.8326.19编号11121314151617181920x坐标(KM)3333344444y坐标(KM)1234512345矿物体厚度H(M)23.2826.4829.1412.0414.5819.9523.7315.3518.0116.29 表一 地区勘探数据表2）问题分析和建立模型 由表一可以知矿物体的厚度H是坐标x和坐标y的二元函数H=H(x,y)，根据二重积分的知识有所求该矿物的储量就是二重积分的值。由于函数H(x,y)没有给出具体

16、表达式，而只给出一些离散点上的函数值，下面采用数值积分的方法来近似计算该矿物的储量。 在计算方法中有很多计算定积分的公式，根据本题特点，这里采用计算定积分的复化梯形公式来计算二重积分。为得到计算二重积分的复化梯形公式，注意到计算定积分的复化梯形公式为这里节点是等距节点xk=a+khx,k=0,1,.,n, hx是相邻节点间的距离。对二重积分令，则有，于是二重积分计算可以由g(x)的复化梯形公式计算，而其中的也是用复化梯形公式计算的，于是可以利用逐次复化梯形求积公式近似计算出来。通过具体推导，可以得到计算二重积分的复化梯形公式： （1）其中是下面矩阵A的元素：3）模型求解 表一的数据是网格化的（

17、数据点位置见图1）图1数据点位置示意图由上图可得：h x= h y =1,n=3，m=4，相应的矩阵利用公式(1) 编程计算有：In1:= da=13.72/4,15.47/2,23.28/2,19.95/4, 25.8/2, 21.33, 26.48,23.73/2, 8.47/2, 14.49, 29.14,15.35/2, 25.27/2, 24.83, 12.04,18.01/2, 22.32/4,26.19/2,14.58/2,16.29/4; daa=4da; Sumdaai,j,i,1,5,j,1,4/4Out1= 244.455考虑到单位问题，得出该矿物的储量约为0.24445

18、510 9(立方米)。本题还可以直接用Mathematica中二元函数分段插值命令求得，具体处理命令为In1:= data = 1, 1, 13.72, 1, 2, 25.80, 1, 3, 8.47, 1, 4, 25.27, 1, 5,22.32, 2, 1, 15.47, 2, 2, 21.33, 2, 3, 14.49, 2, 4,24.83, 2, 5, 26.19, 3, 1, 23.28, 3, 2, 26.48, 3, 3, 29.14, 3, 4, 12.04, 3, 5, 14.58, 4, 1, 19.95,4, 2,23.73, 4, 3, 15.35, 4, 4, 18.01, 4, 5, 16.29; ws=InterpolationdataOut1= InterpolatingFunction1, 4, 1, 5, In2:= Plot3Dwsx,y,x,1,4,y,1,5输出图形为In3:= NIntegratewsx,y,x,1,4,y,1,5Out3= 250.958考虑到单位的关系得出该矿物的储量约为0.25095810 9(立方米)。