微分方程，常微分方程，差分方程模型实例分析

一、微分方程建模

例：美国原子能委员会以往处理浓缩的放射性废料的方法， 一直是把它们装入密封的圆桶里，然后扔到水深为 90 多米的海底。生态学家和科学家们表示担心，怕圆桶下沉到海底时与海底碰撞而发生破裂，从而造成核污染。原子能委员会分辩说这是不可能的。为此工程师们进行了碰撞实验。发现当圆桶下沉速度超过 12.2 m/s 与海底相撞时，圆桶就可能发生碰裂。这样为避免圆桶碰裂，需要计算一下圆桶沉到海底时速度是多少?已知圆桶质量为 239.46 kg，体积为 0.2058m^3，海水密度为 1035.71kg/m3，如果圆桶速度小于 12.2 m/s 就说明这种方法是安全可靠的， 否则就要禁止使用这种方法来处理放射性废料。假设水的阻力与速度大小成正比例，其正比例常数 k = 0.6 。现要求建立合理的数学模型，解决如下实际问题：

（ 1）判断这种处理废料的方法是否合理?

（ 2）一般情况下v 大， k 也大；v 小， k 也小。当 v 很大时，常用 kv 来代替 k ，那么这时速度与时间关系如何? 并求出当速度不超过 12.2 m/s，圆桶的运动时间和位移应不超过多少? ( k 的值仍设为 0.6)

解：

1）

2）

二、常微分方程求解

2.1 解非刚性常微分方程

Matlab 的工具箱提供了几个解非刚性常微分方程的功能函数，如 ode45，ode23，ode113，称其为ode家族！！！其中 ode45 采用四五阶 RK 方法，是解非刚性常微分方程的首选方法， ode23采用二三阶 RK 方法， ode113 采用的是多步法，效率一般比 ode45 高。

Matlab 的工具箱中没有 Euler 方法的功能函数，但是可以自己编函数去解，这里不粘了，主要原因是我忘了欧拉方程是什么。

ode23，ode45，ode113的Matlab的函数形式是[t,y]=solver(‘F’,tspan，y0)。这里solver为ode45，ode23，ode113，输入参数 F 是用M文件定义的微分方程

y’= f (x, y) 右端的函数。tspan=[t0,tfinal]是求解区间，y0是初值。例：

y’= −2y + 2x^2 + 2x ， (0 ≤ x ≤ 0.5) ， y(0) = 1

解：编写函数文件 doty.m 如下：

在Matlab命令窗口输入：

即可求得数值解。
也可以利用Matlab的匿名函数，计算数值解，程序如下：

刚性方程是指一个微分方程，其数值分析的解只有在时间间隔很小时才会稳定，只要时间间隔略大，其解就会不稳定。目前很难去精确地去定义哪些微分方程是刚性方程，但是大体的想法是：这个方程的解包含有快速变化的部分。

Matlab的工具箱提供了几个解刚性常微分方程的功能函数，如ode15s， ode23s，

ode23t， ode23tb，这些函数的使用同上述非刚性微分方程的功能函数。

2.2 求解常微分方程的通解

在 Matlab 中，符号运算工具箱提供了功能强大的求解常微分方程的符号运算命令

dsolve。常微分方程在 Matlab 中按如下规定重新表达：

符号 D 表示对变量的求导。Dy 表示对变量 y 求一阶导数，当需要求变量的 n 阶导数时，用 Dn 表示， D4y 表示对变量 y 求 4 阶导数。由此，常微分方程 y’’+2y’= y 在 Matlab 中，将写成’D2y+2*Dy=y’。

无初边值条件的常微分方程的解就是该方程的通解。其使用格式为：

dsolve(‘diff_equation’)；dsolve(’ diff_equation’， ‘var’)；

式中 diff_equation 为待解的常微分方程， 第 1 种格式将以变量 t 为自变量进行求解，第 2 种格式则需定义自变量 var。

例：x^2 + y + (x − 2y) y’= 0

解:编写程序如下：

三、差分方程求解

规定t 只取非负整数。记 yt 为变量 y 在t 点的取值，则称 Δyt = y(t+1) − y(t) 为 y(t) 的一阶向前差分，简称差分，称 Δ^2 yt = Δ(Δyt ) = Δyt+1 − Δyt = yt+2 − 2yt+1 + yt 为 yt 的二阶差分。类似地，可以定义 yt 的 n 阶差分 Δ^n yt 。

由t、 yt 及 yt 的差分给出的方程称为 yt 的差分方程，其中含 yt 的最高阶差分的阶

数称为该差分方程的阶。差分方程也可以写成不显含差分的形式。例如，二阶差分方程Δ^2 y(t) + Δy(t) + y(t) = 0，也可改写成 y(t+2) − y(t+1) + y(t) = 0 。

满足一差分方程的序列 y(t) 称为差分方程的解。

结果：

我独昏昏：善忘是一种养生