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

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

《22 一些数值方法的建模问题.docx》由会员分享，可在线阅读，更多相关《22 一些数值方法的建模问题.docx（12页珍藏版）》请在冰豆网上搜索。

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.0002

N：

Manning粗糙系数,是一个和水渠内壁材料的光滑度有关的无量刚量,本题中N=0.03

U：

水的流速，单位为米/秒

3） 问题分析和建立模型

查阅水工学使用的有关资料，获得一个有关流速的关系式：

（1）

而质量守衡定律告诉我们

Q=UBH

（2）

将式

（1）代入式

（2），同时将S,N,B,Q用给定的数值替代，有

（3）

令

，于是问题归结为求方程

的正根问题。

4） 模型求解

由于方程

没有求根公式，故采用数值方法求根，这里采用Newton迭代法来求根。

为找到迭代初始值x0，先画出的图形观察：

In[1]:

=f[h_]=Sqrt[0.0002]/0.03*（20*h）^（5/3）*（20+2h）^（-2/3）-5;

In[2]:

=Plot[f[h],{h,0,5}]

输出图形为

Out[2]=-Graphics-

从图中可以看到函数在1附近有根，选取迭代初始值x0=1，根的精度为10-6做Newton迭代，对应的程序为

Clear[x,f,g];

f[x_]=Sqrt[0.0002]/0.03*（20*x）^（5/3）*（20+2x）^（-2/3）-5;

g[x_]=D[f[x],x];

x0=1；

eps1=0.1;

nmax=500;

eps=0.000001;

Do[u1=g[x0];

If[Abs[u1//N]

x=N[x0-f[x0]/u1,10];

u1=Abs[x-x0]//N;

Print["H=",x,"n=",n,"eps=",u1];

If[u1

{n,1,nmax}];

If[u1>eps,Print["迭代失败"]]

执行程序后得输出结果为

H=0.7292262368n=1eps=0.270774

H=0.7025912396n=2eps=0.026635

H=0.702293294n=3eps=0.000297946

H=0.7022932563n=4eps=3.7772⨯10-8

从输出结果说明所求的水渠深度约为0.7022932563，其误差为3.7772⨯10-8．

2.2.2．男大学生的身高问题

1）问题的提出

有关统计资料表明,我国大学生男性群体的平均身高约为170cm，且该群体中约有99.7﹪的人身高在150cm至190cm之间。

如果将[150cm，190cm]等分成20个高度区间，试问该群体身高在每一高度区间的分布情况怎样？

特别地，身高中等（165cm至175cm之间）的人占该群体的百分比会超过60﹪吗？

解：

2）问题分析和建立模型

因为一个人的身高涉及很多因素，通常它是一个服从正态分布N（μ,σ）的随机变量。

正态分布的概率密度函数ϕ（x）为：

于是密度函数ϕ（x）在区间[a,b]上的定积分

的值正好代表大学生身高在区间[acm,bcm]的分布。

在密度函数ϕ（x）中的两个参数μ、σ分别为正态分布的均值和标准差。

根据题目中我国大学生男性群体的平均身高约为170cm，可选取正态分布的均值μ=170cm，而由“该群体中约有99.7﹪的人身高在150cm至190cm之间”和正态分布N（μ,σ）的“3σ规则”，有

μ-3σ=150cm

μ+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.0021

[152cm，154cm]的分布为0.0047

[154cm，156cm]的分布为0.0097

[156cm，158cm]的分布为0.0181

[158cm，160cm]的分布为0.0309

[160cm，162cm]的分布为0.0483

[162cm，164cm]的分布为0.0690

[164cm，166cm]的分布为0.0902

[166cm，168cm]的分布为0.1078

[168cm，170cm]的分布为0.1179

[170cm，172cm]的分布为0.1179

[172cm，174cm]的分布为0.1078

[174cm，176cm]的分布为0.0902

[176cm，178cm]的分布为0.0690

[178cm，180cm]的分布为0.0483

[180cm，182cm]的分布为0.0309

[182cm，184cm]的分布为0.0181

[184cm，186cm]的分布为0.0097

[186cm，188cm]的分布为0.0047

[188cm，190cm]的分布为0.0021

对

（2）式有

≈0.5467246173

说明身高中等（165cm至175cm）的大学生约为54.67﹪，不足60﹪。

4）Mathematica编写的通用程序为

Clear[a,b,x,n,f];

a=Input["a="]

b=Input["b="]

f[x_]=3/20*Exp[-9（x-170）^2/800]/Sqrt[2Pi]

eps=0.0000001;

n=1;

h=b-a;

t1=（f[a]+f[b]）*h/2;

h=h/2;

t2=t1/2+h*f[a+h];

er=t2-t1//N;

While[Abs[er]>eps,

Print["n=",2^n,"定积分值为",N[t2,10]];

Print["误差=",er];

h=h/2;

t1=t2;

n=n+1;

t2=t1/2+h*Sum[f[a+k*h],{k,1,2^n,2}];

er=t2-t1//N];

Print["n=",2^n,"定积分值为",N[t2,10]];

Print["误差=",er]

说明本程序从梯形公式T1开始，用复合梯形求积公式求[a,b]上满足精度小于eps要求的定积分近似值。

程序执行后，按要求通过键盘输入积分下限a、积分上限b、被积函数f（x）和精度要求eps后，计算机则给出满足精度要求的定积分近似值及中间计算值和误差。

程序中变量说明

a：

存放积分下限

b:

存放积分上限

f[x]:

存放被积函数f（x）

eps:

存放求积精度ε

h:

存放节点步长

x：

为函数f[x]:

提供变量

t1:

存放复合梯形值Tn

t2:

存放复合梯形值T2n和定积分近似值

n:

存放复合梯形公式的节点加密次数

er:

存放误差

注：

为提高计算效率，计算T2n时使用了公式

2.2.3湖水温度变化问题

1）问题的提出

湖水在夏天会出现分层现象，其特点为接近湖面的水温度较高，越往下温度变低。

这种上热下冷的现象影响了水的对流和混合过程，使得下层水域缺氧，导致水生鱼类的死亡。

下面是某个湖的观测数据

深度（m）

0

2.3

4.9

9.1

13.7

18.3

22.9

27.2

温度（︒C）

22.8

22.8

22.8

20.6

13.9

11.7

11.1

11.1

请问

1． 湖水在10m处的温度是多少？

2． 湖水在什么深度温度变化最大？

2）问题的分析和假设

本问题只给出了有限的实验数据点，可以想到用插值和拟合的方法来解决题目的要求。

假设湖水深度是温度的连续函数，引入符号如下：

h:

湖水深度,单位为m

T:

湖水温度,单位为︒C,它是湖水深度的函数：

T=T（h）

这里多项式拟合的方法来求出湖水温度函数T（h）。

然后利用求出的拟合函数就可以解决本问题了。

3）模型的建立

将所给数据作图，横轴代表湖水深度，纵轴代表湖水温度，用Mathematica数学软件画出散点图，操作的命令为：

In[1]:

=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.1}}

In[2]:

=q=ListPlot[d,PlotStyle->PointSize[0.04]]

得出如下散点图

Out[2]=-Graphics-

观察散点图的特点，并通过实验选取不同的基函数类进行实验发现用4次多项式拟合比较好，相应的命令为：

In[3]:

=q1=Fit[d,{1,h,h^2,h^3,h^4,h^5},h]

Out[3]=22.711+0.0280756h+0.0865552h2-0.0235702h3+0.00132018h3-0.0000218133h5

In[4]:

=q2=Plot[q1,{h,0,27.2}]

Out[4]=-Graphics-

In[5]:

=Show[q,q2]

Out[5]=-Graphics-

因此有

T（h）≈q1=22.711+0.0280756h+0.0865552h2-0.0235702h3+0.00132018h3-0.0000218133h5

为求出湖水在10m处的温度，键入命令：

In[6]:

=q1/.h->10

Out[6]=19.0975

于是得出湖水在10m处的温度约为19.0975︒C。

为求湖水在什么深度温度变化最大，要求出函数T（h）的导函数T'（h）的绝对值最大值点。

为此对所求拟合函数q1关于h的导数并找出最大值点，键入如下命令：

In[7]:

=q2=D[q1,h]

Out[7]=0.0280756+0.17311h-0.0707107h2+0.00528073h3-0.000109067h4

In[8]:

=Plot[q2,{x,0,27.2}]

输出图形为

Out[8]=

-Graphics-

从导函数图形上观察到其在10附近可以取得绝对值最大值，键入命令：

In[9]:

=FindMinimum[q2,{x,10}]

Out[9]={-1.21357,{h->11.9312}}

得出h=11.9312是导函数的绝对值最大值点，于是可以知道湖水在深度为11.9312m时温度变化最大。

2.2.4储量计算问题

1）问题的提出

某地区为估计某矿物的储量，在该地区内进行勘探，得到如下数据（见表一）。

请估计出此地区内（1≤X≤4，1≤Y≤5）该矿物的储量．

编号

1

2

3

4

5

6

7

8

9

10

x坐标（KM）

1

1

1

1

1

2

2

2

2

2

y坐标（KM）

1

2

3

4

5

1

2

3

4

5

矿物体厚度H（M）

13.72

25.80

8.47

25.27

22.32

15.47

21.33

14.49

24.83

26.19

编号

11

12

13

14

15

16

17

18

19

20

x坐标（KM）

3

3

3

3

3

4

4

4

4

4

y坐标（KM）

1

2

3

4

5

1

2

3

4

5

矿物体厚度H（M）

23.28

26.48

29.14

12.04

14.58

19.95

23.73

15.35

18.01

16.29

表一地区勘探数据表

2）问题分析和建立模型

由表一可以知矿物体的厚度H是坐标x和坐标y的二元函数H=H（x,y），根据二重积分的知识有所求该矿物的储量就是二重积分

的值。

由于函数H（x,y）没有给出具体表达式，而只给出一些离散点上的函数值，下面采用数值积分的方法来近似计算该矿物的储量。

在计算方法中有很多计算定积分的公式，根据本题特点，这里采用计算定积分的复化梯形公式来计算二重积分。

为得到计算二重积分的复化梯形公式，注意到计算定积分的复化梯形公式为

这里节点是等距节点xk=a+khx,k=0,1,...,n,hx是相邻节点间的距离。

对二重积分

令

，则有

，于是二重积分计算可以由g（x）的复化梯形公式

计算，而其中的

也是用复化梯形公式计算的，于是

可以利用逐次复化梯形求积公式近似计算出来。

通过具体推导，可以得到计算二重积分的复化梯形公式：

（1）

其中

是下面矩阵A的元素：

3）模型求解

表一的数据是网格化的（数据点位置见图1）

图1数据点位置示意图

由上图可得：

hx=hy=1,n=3，m=4，相应的矩阵

利用公式

（1）编程计算有：

In[1]:

=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;

Sum[daa[[i,j]],{i,1,5},{j,1,4}]/4

Out[1]=244.455

考虑到单位问题，得出该矿物的储量约为0.244455⨯109（立方米）。

本题还可以直接用Mathematica中二元函数分段插值命令求得，具体处理命令为

In[1]:

=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=Interpolation[data]

Out[1]=InterpolatingFunction[{{1,4},{1,5}},<>]

In[2]:

=Plot3D[ws[x,y],{x,1,4},{y,1,5}]

输出图形为

In[3]:

=NIntegrate[ws[x,y],{x,1,4},{y,1,5}]

Out[3]=250.958

考虑到单位的关系得出该矿物的储量约为0.250958⨯109（立方米）。