数学建模（五） 微分方程，常微分方程，差分方程模型

这部分其实很多教程中没有拎出来讲，但是写论文常用，为了回忆起微分方程是啥，我试图从实验室找出一本高数书，并没有找到，鉴于买一本高数书比较贵，我就从网上下载了部分19年的张宇考研视频，看得特别开心，几乎到了边看边笑的程度。然我觉得并没讲到数学建模要用的重点，直到我看到一本数学建模书里这三个章节的内容！瞬间有一种踏破铁鞋无觅处，得来全不费工夫的感觉！！！

1.微分方程建模

前两个例子是多级火箭和人口模型，倒是看过，举另一例：
美国原子能委员会以往处理浓缩的放射性废料的方法， 一直是把它们装入密封的圆桶里，然后扔到水深为 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）

clc,clear
syms m V rho g k
s=dsolve('m*D2s-m*g+rho*g*V+k*Ds','s(0)=0,Ds(0)=0');
s=subs(s,{m,V,rho,g,k},{239.46,0.2058,1035.71,9.8,0.6});
s=vpa(s,10) %求位移函数
v=dsolve('m*Dv-m*g+rho*g*V+k*v','v(0)=0');
v=subs(v,{m,V,rho,g,k},{239.46,0.2058,1035.71,9.8,0.6});
v=vpa(v,7) %求速度函数
y=s-90;
tt=solve(y) %求到达海底 90 米处的时间
vv=subs(v,tt) %求到底海底 90 米处的速度

2）

clc,clear
syms m V rho g k t real
v=dsolve('m*Dv-m*g+rho*g*V+k*v^2','v(0)=0');
v=subs(v,{m,V,rho,g,k},{239.46,0.2058,1035.71,9.8,0.6});
v=vpa(v,7) %求速度函数
v1=real(v),v2=imag(v) %求速度函数的实部和虚部
tt=solve(v1-12.2) %求时间的临界值 T
s=int(v1,0,tt) %求位移的临界值

数学部分的过程不粘了

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 如下：

function f=doty(x,y);
f=-2*y+2*x^2+2*x;

在Matlab命令窗口输入：

[x,y]=ode45('doty',[0,0.5],1)

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

f=@(x,y)-2*y+2*x^2+2*x;
[x,y]=ode45(f,[0,0.5],1)

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

syms x y
diff_equ='x^2+y+(x-2*y)*Dy=0';
dsolve(diff_equ,'x')
 
ans =
 
 x/2 + ((4*x^3)/3 + x^2 + C1)^(1/2)/2
 x/2 - ((4*x^3)/3 + x^2 + C1)^(1/2)/2

3.差分方程求解

规定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) 称为差分方程的解。


结果：

yk =
 
 (-1)^k - (3*0^k)/2 - (-2)^k
     0^k - (-1)^k + 2*(-2)^k