机理分析建模讲解

当一个数值公式的截断误差可表示为O（hk+1） （其中k为正整数，h为步长）时，称它是一个k阶公 式． k越大，则数值公式的精度越高．
•欧拉法是一阶公式，改进的欧拉法是二阶公式． •龙格-库塔法有二阶公式和四阶公式． •线性多步法有四阶亚当斯外插公式和内插公式．
1.3 用MATLAB软件求常微分方程的数值解
在实际问题中， “改变”、“变化”、“增加”、 “减少”等关键词提示我们注意什么量在变化；关键词 “速率”、“增长” “衰变” ，“边际的” ，常涉及 到导数。这些都是建立微分方程模型的关键。
(一) 微分方程的建立
建立常微分方程模型的常用方法：
运用已知物理定律 利用平衡与增长式 运用微元法 运用分析法
dt
M
假设1
市场余额
假设2
销售速度因广告作用增大, 同时 又受市场余额的限制。
(二) 微分方程的求解
求解常微分方程模型的常用方法： 微分方程的数值解 微分方程的定性分析
1、常微分方程的数值解
1.1 常微分方程数值解的定义：
在生产和科研中所处理的微分方程往往很复杂,且 大多得不出一般解。而实际问题中对初值问题的求 解，一般是要求得到在若干个点上满足规定精确度 的近似值，或者得到一个满足精确度要求的便于计 算的表达式。
即：对常微分方程

y' y ( x0
f )
(x,yy0 )，其数值解是指由初始点x0
开始
的若干离散的x处的值，即对x0 x1 x2 xn，求出准确值y(x1 ),
y(x2 ), , y(xn ) 的相应近似值y1, y2 , , yn .
1.2 建立数值解法的一些途径： 用差商代替导数 使用数值积分 使用泰勒公式 数值公式的精度
例1.3 一个高为2米的球体容器里盛了一半的水，水 从它的底部小孔流出，小孔的横截面积为1平方厘米。 试求放空容器所需要的时间。
对孔口的流速做两条假设 ：
1．t 时刻的流速v 依赖于此
2米
刻容器内水的高度h(t)。
2 ．整个放水过程无能量损失。
分析：放空容器
容器内水的体积为零 容器内水的高度为零
模型建立：由水力学知：水从孔口流出的流量Q为 通过“孔口横截面的水的体积V对时间t 的变化率”，即
此类建模方法的关键是分析并正确描述基本模型的 右端，使平衡式成立。
例1.2(战斗模型) 两方军队交战，希望为这场战斗 建立一个数学模型，应用这个模型达到如下目的：
1. 预测哪一方将获胜？
2. 估计获胜的一方最后剩下多少士兵？ 3. 计算失败的一方开始时必须投入多少士兵才能 赢得这场战斗？ 模型建立
机理分析建模常用方法： 常微分方程 偏微分方程 逻辑方法 比例方法 代数方法
常微分方程建模 微分方程的建立 微分方程的求解 逻辑方法建模
目录
一 微分方程建模
当实际问题需寻求某个变量y 随另一变量 t 的变化 规律 y=y(t)，且直接求很困难时，可以建立关于未知变 量、未知变量的导数以及自变量的方程(即变量满足的 微分方程)。
1、运用已知物理定律
建立微分方程模型时应用已知物理定律，可事半功 倍。
例1.1 一个较热的物体置于室温为180C的房间内， 该物体最初的温度是600C，3分钟以后降到500C 。想知 道它的温度降到300C 需要多少时间？10分钟以后它的 温度是多少？
牛顿冷却(加热)定律：将温度为T的物体放入处于 常温 m 的介质中时，T的变化速率正比于T与周围介质 的温度差。
dT
k(T
m)
dt
T (0) 60
其中参数k >0，m=18，求得一般解为
ln(T－m)=－k t+c 或 T m cekt (t 0)
代入条件，求得c=42 ，
k


1 3
ln
16 21
,
最后得
1 ln 16 t
T (t ) 18 42e 3 21 (t 0)
发射导弹，导弹头始终对准乙舰．如果乙舰以最大的速度
v0(常数)沿平行于y轴的直线行驶，导弹的速度是5v0，求导
弹运行的曲线方程．乙舰行驶多远时，导弹将它击中？
解法一（解析法）
y12 ) y2

y1

y1
(0)

2,
y2
(0)

0
1．建立M文件vdp1000．m如下： function dy=vdp1000(t,y) dy=zeros(2,1); dy(1)=y(2);
dy(2)=1000*(1-y(1)^2)*y(2)-y(1); 2．取t0=0，tf=3000，输入命令：
设: x(t) — t 时刻X方存活的士兵数； y(t) — t 时刻Y方存活的士兵数；
假设： 1) 双方所有士兵不是战死就是活着参加战斗， x(t)与 y(t)都是连续变量。
2) Y方军队的一个士兵在单位时间内杀死X 方军队 a 名士兵；
3) X 方军队的一个士兵在单位时间内杀死Y方军队 b 名士兵；
Q dV 0.62S 2gh dt
S—孔口横截面积(单位：平方厘米) h(t) —水面高度(单位：厘米) t—时间(单位：秒) 当S=1平方厘米，有
dV 0.62 2ghdt (1)
r1
h(t) r2
h+Δh
在[t，t+Δt ]内，水面高度 h(t) 降至h+Δh(Δh<0), 容 器中水的体积的改变量为
[T,Y]=ode15s('vdp1000',[0 3000],[2 0]); plot(T,Y(:,1),'-')
3．结果如图：
2 1.5
1 0.5
0 -0.5
-1 -1.5
-2 -2.5
0
500
ห้องสมุดไป่ตู้1000
1500
2000
2500
3000
例1.5 导弹追踪问题
设位于坐标原点的甲舰向位于x轴上点A(1, 0)处的乙舰
h t0 100
积分后整理得
35
t (700000 1000h2 3h2 )
4.65 2g
(0≤h≤100)
令 h=0，求得完全排空需要约2小时58分。
4、分析法
基本思想：根据对现实对象特性的认识，分析其因 果关系, 找出反映内部机理的规律。
例1.4(独家广告模型) 广告是调整商品销售的强有 力的手段, 广告与销售量之间有什么内在联系？如何评 价不同时期的广告效果？
2．使用MATLAB软件求数值解时，高阶微分方程必须 等价地变换成一阶微分方程组．
d 2x 例： dt 2
1000(1
x2
)
dx dt

x

0
x(0) 2, x' (0) 0
解：令y1 x, y2 y1'
则微分方程变成一阶微分方程组：


y1' y 2'
y2 1000(1
[t，x]=solver(’f’,ts,x0,options)
自变 函数 量值 值
ode45 ode23 ode11 3ode1 5sode 23s
由待解
方程写 成的M 文件名
ts=[t0,t f]，t0、
tf为自变
量的初值 和终值
函数 的初 值
ode23：组合的2/3阶龙格–库塔–费尔贝格算法
ode45：运用组合的4/5阶龙格–库塔–费尔贝格算法
ΔV=V(h)－V(h+Δh)=－πΔh[3(r12+r22)＋o(Δh)] ≈－πr2Δh＋o(Δh)
记 r 1002 (100 h)2 200h h2
令Δt 0, 得 dV=－πr2 dh， (2)
比较(1)、(2)两式得微分方程如下：
0.62 2ghdt (200h h2 )dh
在很短的时间段Δt 内，关于P(t)变化的一个最简单 的模型是：
{Δt时间内的人口增长量} ={Δt内出生人口数}－{Δt内死亡人口数}
+ {Δt内迁入人口数}－{Δt内迁出人口数} 更一般地
{Δt时间内的净改变量} ={Δt时间内输入量}－{Δt时间内输出量} 不同的输入、输出情况对应不同的差分或微分方程。 输入量：含系统外部输入及系统内部产生的量； 输出量：含流出系统及在系统内部消亡的量。
用差商代替导数（欧拉法）
设 xi1 xi h, i 0,1, 2,
微分方程
y ' f (x, y)

y(
x0
)

y0
, n 1, 则可用以下离散化方法求解
若步长h较小，则有
y'(x) y(x h) y(x) h
故有公式：

yi1 yi hf y0 y(x0 )
结果：
T (10)

18
1 ln 1610
42e 3 21

39.3(0 C )
该物体温度降至300C 需要8.17分钟。
2、利用平衡与增长式
许多研究对象在数量上常常表现出某种不变的特性， 如封闭区域内的能量、货币量等。
利用变量间的平衡与增长特性，可分析和建立有关 变量间的相互关系.
续 人口增长模型 对某地区时刻t的人口总数P(t)，除考虑个体的出生、 死亡，再进一步考虑迁入与迁出的影响。
平衡式：
{Δt 时间内X军队减少的士兵数 } = {Δt 时间内Y军队消灭对方的士兵数}
即有：Δx =－ayΔt ，同理：Δy =－bxΔt
令Δt 0，得到微分方程组：
dx