[个人总经验] 数模中平衡原理与机理模型

[复制链接]

字体大小: 正常放大

wenxinzi

6 主题	3 听众	51 积分

升级 48.42%

TA的每日心情

	开心 2016-11-7 00:15

签到天数: 7 天

[LV.3]偶尔看看II

电梯直达

1^#

发表于 2011-9-8 12:44 |只看该作者 |倒序浏览

|招呼Ta 关注Ta

§3.3  平衡原理与机理模型
一. 平衡原理
   自然界任何物质在其运动变化过程中一定受到某种平衡关系的支配。
二. 机理模型
   在一定的假设下，根据主要因素相互作用的机理，对它们之间的平衡关系的数学描述。
三. 连续模型
连续模型组建的微元法
在自变量的微小的区间内以简单的形式描述有关变量之间的平衡关系, 再利用微分学的思想进一步处理它, 得到以微分方程的形式描述的数学模型。

例1. 人口的自然增长.
建模描述一个地区内人口的自然增殖的过程。即考虑由于人口的生育和死亡所引起的人群数量变化的过程。
假设1. 人群个体同质。令N(t)表示t时刻的人口数。
假设2.  群体规模大。 N(t) 连续可微.
假设3. 群体封闭，只考虑生育和死亡对人口的影响。
平衡关系：人口数在区间[t，t+ t ]内的改变量等于这段时间内出生的个体数与死亡的个体数之差。令B(t, t, N), D(t, t , N) 分别表示在时间区间[t，t+ t ]内生育数和死亡数, 则有
  N(t+t)-N(t)=B(t,  t,N)-D(t,  t,N)
假设4. 从大群体的平均效应考虑生育和死亡对人口的影响。（生育率和死亡率）
生育率b(t, t, N) = B(t, t, N)/N, 死亡率d(t, t, N) = D(t, t, N)/N
记增长率为 R(t,  t,N)= b(t,  t,N)-d(t,  t,N) 则有 N(t+t)-N(t)=R(t,  t,N)N
将R(t, t,N)关于t展开. 由于R(t, h, N)|h=0=0，所以

N(t+ t ) - N(t)=r(t,N)N t +o(t).
两边除以t, 并令t →0,  得到 dN/dt=r(t, N)N
假设5. 群体增长恒定。（r与 t 无关） dN/dt=r(N) N
假设6. 个体增长独立。（r 与 N 无关）  dN/dt=r N
给定初值 N(0)=N0，可得人口增长的指数模型(Maithus  模型)  N(t)=N0ert
在离散时间点k=0, 1, 2, …, 上有 N(k+1) = er N(k )
Maithus:  “若我的两个假设是成立的,那么,我认为人口繁殖的能量是无限地大于自然界为人类提供资料的能量的。人口如果不受控制，它会以几何比率增长。而生活资料只能以算术比率增长。只要稍微看一下数字，就将明确第一种能量比之第二种能量是无比巨大的。” 《论人口原理》
总结对人口指数增长模型的假设,
1. 人群个体同质。
2. 群体规模大。
3.  群体封闭，只考虑生育和死亡对人口的影响。
4. 从大群体的平均效应考虑生育和死亡对人口的影响。（生育率和死亡率）
5. 群体增长恒定。
6. 个体增长独立。
由这些假设可分析这个模型的作用.

例2 池水含盐
池中有一定体积的盐水，从池的上部向池中注入一定浓度的盐水。混合后的盐水将从池的下部流出。建模描述池中盐水浓度的动态。
假设：盐水注入池中后迅速混合, 使得盐水浓度均匀。
变量、参量：
池中盐水体积 V(t), 池中盐水浓度 p(t);
池中原有盐水体积 V0, 原有盐水浓度p0;
流入盐水速度 rI(t), 流入盐水浓度 pI(t);
流出盐水速度rO(t), 流出盐水浓度 p(t).

平衡关系
在时间段[t+  t]内,
池中(纯)盐的改变量=这段时间内流入的(纯)盐的量与流出的(纯)盐的量之差。
池中盐的改变量：  p(t+t)V(t+  t)-p(t)V(t)
流入盐量：             流出盐量：

利用积分中值定理可得

在时间段[t+  t]内,
池中盐水体积的改变量=这段时间内流入盐水的体积与流出盐水体积之差；

令  t 0 得

模型

进一步问题：池中有水 2000 m3，含盐 2 kg，以 6m3 / 分的速率向池中注入浓度为 0.5 kg / m3 的盐水，又以 4 m3 / 分的速率从池中流出混合后的盐水。问欲使池中盐水浓度达到 0.2 kg / m3，需要多长时间？
此时 V(t)=2000+2*t.  dp/dt=3/V(t)-6*p(t)/V(t),  p(0)=0.001.
用MATLAB求 p(t)
求表达式（符号运算）
S=dsolve(‘Dx=(3-6*x)/(2000+2*t)’);
求数值解
建立M文件 fun . M, function  y=fun(t,x)
y=(3-6*x)/(2000+2*t);
t0=0; tf=200; x0=2; [t,x]=ode23(‘fun’,t0,tf,x0);
plot(t,x);

四. 离散模型
离散模型的组建
利用平衡原理，找出每一步对前一步或前几步的依赖关系，得到以差分方程的形式描述的数学模型。

例1. 买房贷款：银行可以向购房人提供个人住房贷款的业务。偿还贷款时要求借款人在借款期间内每月以相等的月均还款额偿还银行贷款本金和利息。试组建计算月均还款额的数学模型。
假设：1. 逐月偿还贷款；2. 每月还款金额相等；3. 按月计算利息；4. 每月月底还款。
参量、变量
贷款额：A(万元)，贷款期限：N年(n=12N月) , 月利率：r，月均还款额：x。
令Ck表示第k月月底还款后的欠款余额，记C0=A.
第n月的月底欠款应全部偿还完毕,  则有Cn=0
平衡关系:
本月月底还款后的欠款余额=上月欠款余额的本利和扣除本月还款后的金额。
模型：       Ck= (1+r)Ck-1-x
求解：递推可得 Ck = (1+r)kA-ki=0(1+r)ix = (1+r)kA-{[(1+r)k-1]/r}x
于是  0=(1+r)nA-x[(1+r)n-1]/r,  所以  x=A r (1+r)n/[(1+r)n-1]

例2. 兔子的繁殖I
   由一对兔子开始，一年可以繁殖成多少对兔子？假设兔子的生殖力是这样的：一对兔子每一个月可以生一对兔子，并且兔子在出生两个月以后就具有繁殖后代的能力。
n 1 2 3 4 5 6 7 8 9  10 11 12
a(n)  1 1 2 3 5 8  13  21  34  55 89  144
      a(n+1) = a(n) + a(n-1) 斐波那契数列(黄金数)
假设：
1. 每对兔子每一个月定生一对兔子。
2. 兔子出生两个月后都具有繁殖能力。
3. 兔子每经过一个月底就增加一个月令。
变量、参量:  月份：n，幼兔：a0(n)，成兔:  a1(n)
平衡关系
本月初(一月令)的幼兔是上月成兔繁殖的后代。
本月的成兔是上月的成兔和上个月(一月令)的幼兔发育结果的总和。
模型 I  a0(n) = a1(n-1) a1(n) = a0(n-1) + a1(n-1)
令 a(n) = (a0(n), a1(n))’,  则 a(n) = A a(n-1)

分析
1. 模拟. a0(1)=1,  a1(1)=0
      n 1 2 3 4 5 6 7 8 9  10  11  12
   a0(n)  1 0 1 1 2 3 5 8  13  21  34  55
   a1(n)  0 1 1 2 3 5 8  13  21  34  55  89
   a(n) 1 1 2 3 5 8  13  21  34  55  89  144
  2. 证明 a(n+1) = a(n) + a(n-1)
   因为 a0(n+1)=a1(n) a1(n+1)=a(n)=a0 (n)+a1 (n)
   所以  a(n+1)= a0 (n+1)+a1 (n+1) = a1(n)+a(n) = a(n-1)+a(n)
  3. 模型的作用机理：a11幼兔的繁殖能力， a12成兔的繁殖能力，a21幼兔的发育为成兔的比例，
a22成兔存活的比例。
4. 群体的渐近性质
      A 有主特征值 =1.618  相应的左特征向量 L=(0.382 0.618)’，
      于是，当n  时,       a(n)/a(n) L

例3. 录音机的运行
建模分析磁带录音机的运行规律(计数器的读数与运行时间的关系)。
数据：I. 读数与时间（秒）
t 1 2 3 4 5  10  15 20 25 30 31
n 9  18  28  37  47  97  151  211  280  362  382 385
   数据：II. 读数与转数
　k  2  4  10  14  18  22  26  31  35  41  60
　ｎ 1  2  5 7 9 11  13  15  17  20  29　　
时间与读数关系的散点图                            读数与转速的散点图

背景
1. 磁带盒内有二个磁带轮：送带轮和收带轮。放音时送带轮上的磁带减少，缠于收带轮上。
2. 计数器只记录某个磁带轮转动的转动情况。计数器的读数不刚好是磁带轮的转数。
3  磁带轮在放音时转动不是匀速的，送带轮加速，收带轮减速。
4  通过磁头时，磁带匀速运行。
假设
1. 计数器记录了送带轮的转数k。
2. 计数器的读数 n 与送带轮的转数 k成正比。
3. 磁带运行的线速度定常。
4. 磁带厚度均匀，缠绕松紧一致，无空隙。
5. 磁带缠绕一圈的周长等于缠绕的圆周长。
参量、变量：计数器读数：n，带轮转数：k (=从外向里第 k 圈),  运行时间：t(k)，磁带厚度：d，带芯轮半径：r，磁带速度：v ，从外向里第 k 圈磁带的半径：Rk，从外向里第k圈磁带长度：L(k) ，磁带最多圈数：N
平衡关系：运行k圈磁带的时间等于磁带的长度与运行速度之商。
分析：t=0时n=0，送带轮缠满磁带并开始转动。由假设5，Lk=2Rk 由假设4，Rk=r+(N-k+1)d，最外k圈磁带总长度
L(k)=I=1k 2 [r+(N-i+1)d]= 2[kr +kNd-(1+2+…+(k-1))d]= 2[kr +kNd-k(k-1)d/2]
=(2r +2Nd+d)k-k2d 
又由假设2,  k = c n．则有t(n)=c(2r+2Nd+d)n/v- c2n2d/v
模型：　t(n) = a n + b n2,  其中a= c(2r+2Nd+d)/v, b= - c2d/v
参数 a, b, c 的估计：
1. 用最小二乘法估计 a, b，得到  a = 0.11095，b = -7.7445•10-5
模型：
      t(n) = 0.11095 n – 7.4475• 10-5 n2
检验
n 9 18 28 37 47 97 151 211 280
t 0.99  1.97  3.05  4.00  5.04  10.03  14.99  19.96  24.99
t 1    2 3    4 5    10 15    20    25

2. 根据假设2：k = c n,  利用数据II可以给出参数 c 的最小二乘估计 c = 2.04。
又可测得 r = 1.1 cm  计算得N=385×2.04=785, 由a= c(2r+2Nd+d)/v, b= - c2d/v
可以求出d = 0.001628cm，v = 2.75m/min. 最后计算得到 L=2.75×31 =85.25m。

习题
P88 习题：  14， 15.
问题一. 根据录音机运行的数学模型及观测数据 I 给出模型参数最小二乘估计的正规方程组。问题二. 在数据 I、II 的基础上，使用MATLAB 给出录音机运行模型的数值分析。
问题三:  1997年11月8日电视正在播放十分壮观的长江三峡工程大江截流的实况。截流从8：55开始，当时龙口水面宽40m，水深60m。到11：50时，播音员报告宽为34。4m。到13：00时，播音员又报告水面宽31m。这时电视机旁的小明说，现在可以估算下午几点合龙。从8：55到11：50，进展的速度为每小时宽度减少1.9m。从11：50到13：00，每小时宽度减少2.9m。小明认为回填速度是越来越快的，近似地每小时速度加快1m。从下午1：00起，大约要5个多小时，即到下午6点多才能合龙。但到了下午3点28分，电视里传来了振奋人心的消息：大江截流成功！小明后来想明白了，他估算的方法不好。现在请你根据上面的数据设计一种合理的估算方法（建立一种合理的数学模型）进行估算，使你的计算结果更切合实际。

zan