微分方程模型|凶案时间推断|欧拉法|Vanderpol方程

微分方程简介

在研究某些实际问题时，希望能得到问题内部特征在数量方面的函数关系，利用此关系对实际问题的规律进行研究。但是在实际中往往很难或者无法直接得到各内部特征之间的联系，问题的特性反而表现为相关变量的变化率之间的关系。利用这些关系，我们可以建立相应变量之间的微分方程模型
在自然界以及工程技术领域中，微分方程模型是以客观规律为基础，已经渗透到经济学、人口问题以及其他社会科学领域，影响非常广泛。

微分方程的数学描述

微分方程：联系自变量、未知函数及未知函数导数的方程
$$\frac{dx}{dt}=f(x,t)$$

• t，自变量
• x，未知函数，因变量
• f(x,t)，函数关系

1. 若自变量个数为一个，称为常微分方程(ODE）
2. 若自变量个数为两个或两个以上，称为偏微分方程（PDE）
3. 未知函数导数实际出现的最高阶数n称为微分方程的阶
   $$\frac{d^{2}x}{d{t}^2}=a二阶微分方程$$
   $$\frac{\partial u}{\partial t}=a^2\frac{{\partial }^{2}u}{\partial x^2},其中u(x,t)（热传导方程）$$

微分方程的解

若微分方程中含有任意常数，且独立的任意常数个数与方程的阶数相等，则称这个解为微分方程的通解
$$x=\frac{1}{2}a{t}^2+C_{1}t+C_2是微分方程\frac{d^{2}x}{d{t}^2}=a的通解$$
$$\{\begin{array}{c}\frac{dx}{dt}=f(x,t)\\ x(t_0)=x_0\end{array},称为常微分方程的初值问题$$
$$x(t_0)=x_0此条件称为初值条件$$

微分方程的初值问题

初值问题主要讨论的问题有：
初值问题是否存在解；解的存在域有多大；解是否唯一；当初值变化时，解如何变化；等等。

• 定量
  解有初等函数表达式（解析解）
  数值解：计算机求解
• 定性
  稳定性

微分方程建模基本方法

微分方程模型是一类应用十分广泛而且最常见的数学模型，其建模方法在数学建模应用中占有重要的地位

1. 微分方程模型是机理分析建模方法的最佳体现；
2. 微分方程模型是物理、生物进化等定律最精确的定量描述

微分方程建模基础

1. 根据对象的特征和研究目的，对问题进行必要的简化和假设，将假设用精确的语言给出适合的数学语言表达。然后根据假设分析对象的因果关系，利用对象的内在规律建立各个量之间的等式关系或其他数学结构
2. 基于各种已有的确定性的、一般性的机理和理论，结合实际情况建立微分方程模型

微分方程建模准则

1. 确定影响变化的主要因素：确定各因素之间的因果关系，也就是确定一种函数关系，根据内在的机理分析确定各种量之间的平衡关系
2. 改变量=输入量－输出量：只要明确什么在变，如何变，理清楚输入量输出量有哪些，如果输入量和输出量在数值上不相等，微分模型就派上用场

微分方程建模过程

1. 转化翻译：
   问题的表面分析。有许多表示导数的常用词，如：速率，增长、衰变、边际、弹性等。改变、变化、增加、减少这些词可能是一种暗示信号，只需清楚什么在变，随什么而变，这时可考虑用微分方程建模
2. 机理分析：
   基于基本的原理或物理定律，作合理的推导与类比，如：能量不会消失，只是互相转化.要遵循物质守恒，能量守恒等已知公认的定律
   建立平衡式：
   变化率 = 输入率 - 输出率
3. 微分方程模型：
   微分方程是在任何时刻必须正确的瞬时表达式。如看到了表示导数的关键词，就要寻找$y^{\prime }(t)$与$y(t)$，$t$的关系。首先将注意力集中在文字形式的关系式：变化率 = 输入率 - 输出率，然后准确填好式中的所有项.
4. 单位：
   采用统一的物理单位确保等式平衡；采用物理学中约定熟成的数学符号
5. 定解条件：
   系统在某一特定时刻的信息，独立于微分方程而成立。利用它们来确定有关的常数，包括比例系数、原微分方程的其它参数以及解中的系数.

凶案时间推断

1. 问题分析
   该问题符合物理学的冷却现象，可以应用Newton冷却定律：“物体在介质中的冷却速度同该物体温度与介质温度之差成正比”来解决。
   用速度来描述物体在某时刻内的变化率，涉及到导数概念，反映在数学模型上就可以运用微分方程来建模
2. 模型建立
   按一般情形考虑，记$T_t$为时刻$t$物体的温度，$T_0$为初始时刻$t_0$物体的温度（本例中为受害人被害时的体温），$T_e$为介质的温度（环境温度），则由Newton冷却定律可得一阶线性微分方程模型
   $$\{\begin{array}{c}\frac{d{T}_t}{dt}=-\lambda (T_t-T_e)\\ T_{t_0}=T_0\end{array}$$
   其中$\lambda >0$为比例系数，由物体和介质的性质决定，负号则表示温度是下降的

• 尸体的温度随时间的变化率与温差成正比
• 同时被害人刚刚被害的时候，体温是37度

模型求解

对上述模型采用分离变量法求解，可得
$$T_t=(T_0-T_e)e^{-\lambda (t-t_0)}+T_e$$

• $T_t$，任意时刻尸体的温度
  这就是物体冷却过程中，物体温度随时间变化的函数关系，再根据物体与介质的性质决定入值之后，利用$T_0,T_e,T_t$等已知条件，可以得到便于应用的形式
  $$t-t_0=\frac{1}{\lambda }\ln \frac{T_0-T_e}{T_t-T_e}$$
• $t-t_0$，时差

微分方程解析解及MATLAB实现

微分方程的解析解

Matlab软件求解微分方程（组）的解析解命令格式：
自变量

dsolve('方程1','方程2',...,'方程n','初始条件','自变量')

$$\frac{dx}{dt}=x^3+1$$

Dx = x^3 + 1

在Matlab中，字母D表示求1阶微分，D2、D3分别表示求二、三阶微分。
D后面跟随的字母为因变量。自变量可以指定或者缺省.缺省情况下默认为t

例1 一阶常系数微分方程
$$\frac{dy}{dx}=1+y^2,y(0)=1$$
输入：

y1 = dsolve('Dy = 1 + y^2', 'x')
y2 = dsolve('Dy = 1 + y^2', 'y(0) = 1', 'x')

输出

y1 = tan(x - C1)    %通解
y2 = tan(x - PI/4)  %特解

例2 二阶常系数微分方程
$$y^{\prime \prime }-2y^{\prime }-3y=0$$
$$y(0)=1,y^{\prime }(0)=0$$
输入

y1 = dsolve('D2y - 2*Dy - 3*y = 0', 'x')
y2 = dsolve('D2y - 2*Dy - 3*y = 0', 'y(0) = 1, Dy(0)=0', 'x')

输出

y1 = C1 * exp(-x) + C2 * exp(3*x)  %通解
y2 = 3/4 * exp(-x) + 1/4 * exp(3*x) %特解

例3 二阶非常系数微分方程
$$\{\begin{array}{c}{x}^{\prime \prime }(t)-(1-x^2(t))x^{\prime }(t)+x(t)=0\\ x(0)=3,x^{\prime }(0)=0\end{array}$$
输入

x = dsolve('D2x - (1 - x^2) * Dx + x = 0', 'x(0)=3,Dx(0)=0', 't')

输出
无解析解表达式

微分方程组的解析解

dsolve的输出结果个数只能是一个或者与方程组个数相等
![[Pasted image 20240822095644.png]]
输入Matlab命令：

S = dsolve('Df = f + g','Dg = g - f','Df(0) = 1,Dg(0) = 1')

输出

s = g:[1 x 1 sym]
	f:[1 x 1 sym]

输入：

simplify(S.f)

输出

ans = exp(t) * sin(t)

输入

simplify(S.g)

输出

ans = exp(t) * cos(t)

例5

也可以写成：

r = dsolve('Dx + 5 * x = 0', 'Dy - 3*y = 0', 'x(0) = 1', 'y(0) = 1′,'t')

这里r是一个结构类型的数据，可以用以下方式查看其结构内部数据：

r.x    %查看解函数x(t)
r.y    %查看解函数y(t)

只有很少一部分微分方程（组）能求出解析解
大部分微分方程（组）只能采用数值方法求得数值解

例6
求微分方程组的通解
$$\{\begin{array}{c}\frac{dx}{dt}=2x-3y+3z\\ \frac{dy}{dt}=4x-5y+3z\\ \frac{dz}{dt}=4x-4y+2z\end{array}$$
输入命令

[x,y,z] = dsolve('Dx = 2*x - 3*y + 3*z', 
				'Dy = 4*x - 5*y + 3*z', 
				'Dz = 4*x - 4*y + 2*z)

结果

x = C3 * exp(2 * t) + exp(-t) * C1
y = C2 * exp(-2 * t) + C3 * exp(2 * t) + exp(-t) * C1 
Z = C2 * exp(-2 * t) + C3 * exp(2 * t)

微分方程数值解原理

常微分方程数值解的意义

在生产和研究中所处理的微分方程往往很复杂，且很多情况下得不到解析解。而实际生活中的初值问题，往往要求得到解在若干个点上满足一定精确度的近似值，或者得到一个满足精确度要求且便于计算的表达式。因此，研究常微分方程的数值解法是十分必要的
对常微分方程：
$$\{\begin{array}{c}{y}^{\prime }=f(x,y)\\ y(x_0)=y_0\end{array}$$
而言，其数值解是指：
在初始点$x_0$开始的若干个离散的x值处，即对一系列x值：$x_0<x_1<\cdots <x_n$，求解准确的值：$y(x_1),y(x_2),\dots ,y(x_n)$的相应近似值，$y_1,y_2,\dots ,y_n$

建立数值解的一般方法

1. 用差商代替导数
   若步长h较小，则有如下变换
   $$y^{\prime }(x_i)\approx \frac{y(x_{i+1}-y(x_i))}{x_{i+1}-x_i}=\frac{y(x_{i+1}-y(x_i))}{h}$$
   可得到公式
   $$\{\begin{array}{c}{y}_{i+1}=y_i+hf(x_i,y_i)\\ y_0=y(x_0)\end{array},i=0,1,2,\dots ,n-1$$
   欧拉公式
   $$\{\begin{array}{c}\frac{dy}{dx}=y+\frac{2x}{y^2}\\ y(0)=1\end{array}$$
   解析解
   $$y={\left(\frac{5}{3}{e}^{3x}-2x-\frac{2}{3}\right)}^{1/3}$$
   

使用数值积分方法

对方程$y^{\prime }=f(x,y)$两边由$x_i$到$x_{i+1}$积分，并利用梯形公式
$$y(x_{i+1})-y(x_i)={\int }_{x_i}^{x_{i+1}}f(t,y(t))dt\approx \frac{x_{i+1}-x_i}{2}[f(x_i,y(x_i))+f(x_{i+1},y(x_{i+1}))]$$
可以得到
$$\{\begin{array}{c}{y}_{i+1}=y_i+\frac{h}{2}[f(x_i,y_i)+f(x_{i+1},y_{i+1})]\\ y_0=y(x_0)\end{array}$$
实际应用时与欧拉公式结合使用
$$\{\begin{array}{c}{y}_{i+1}^{(0)}=y_i+hf(x_i,y_i)\\ y_{i+1}^{(k+1)}=y_i+\frac{h}{2}[f(x_i,y_i)+f(x_{i+1},y_{i+1}^{(k)})],k=0,1,2\dots \end{array}$$
对给定的精度$\epsilon$，当满足
$$|y_{i+1}^{(k+1)}-y_{i+1}^{(k)}|<\epsilon$$
取
$$y_{i+1}=y_{i+1}^{(k+1)}$$
然后继续下一步$y_{i+2}$的计算，这就是改进的欧拉方法

使用泰勒公式

以此方法为基础，有龙格-库塔法、线性多步法等方法。
数值公式的精度：当一个数值公式的截断误差可表示为$O(h^{k+1})$时（k为正整数，h为步长），称它是一个k阶公式。
k越大，则数值公式的精度越高

• 欧拉法是一阶公式，改进的欧拉法是二阶公式。
• 龙格-库塔法有二阶公式和四阶公式
• 线性多步法有四阶阿达姆斯外插公式和内插公式

欧拉法与龙格库塔(R-K)法误差比较

微分方程数值解的Matlab实现

常微分方程组数值解介绍

[T, Y] = solver('f', ts, Y0, options)

• T，自变量值
• Y，函数值
• solver，ode45 ode23 ode113 ode15s ode23s
  ode23：组合的$\frac{2}{3}$阶龙格-库塔-芬尔格算法
  ode45：运用组合的$\frac{4}{5}$阶龙格-库塔-芬尔格算法
  ode113：多步法
  ode15s：多步法；Gear’s反向数值微分，精度中等
  ode23s:一步法，2阶Rosebrok算法，低精度
• f，由待解方程写成的m文件名
• ts，$ts=[t_0,t_f],t_0,t_f$为自变量区间的初值和终值
• Y0，函数的初值
• options
  用于设定误差限（缺省时设定相对误差$10^{-3}$，绝对误差$10^{-6}$)
  命令格式为：

options = odeset('reltol', 'rt', 'abstol', at)

其中：rt，at：分别为相对误差、绝对误差

Matlab数值解求解方法详解

求解命令格式一般写成：

[T, Y] = solver('f', ts, Y0)

输入部分：
f，微分公式；ts，求解区间，Y0，初值条件；
solver：即ODE求解器（可用的求解器如：ode45、ode23、ode113 ode15s、ode23s、0de23t、0de23tb)
输出部分：
Matlab在求解时，自动将求解区间分割成离散点，T(向量）返回分割点的值（自变量），Y（向量）返回解函数在对应分割点上的求得的数值.

一般没有一种算法可以有效地解决所有的ODE问题，因此MATLAB提供了多种ODE求解器，对于不同的微分方程，可以调用不同的求解器

ODE求解器简要说明

刚性问题与积分步长的关系

有一类常微分方程（组），在求数值解时遇到相当大的困难，这类问题
的解的分量有的变化很快，有的变化很慢，也就是说：
变化快的分量很快趋于它的稳定值，而变化慢的分量缓慢地趋于它的稳定值
以数值解的角度看，当解变化快时，积分要小步长；当解的变化缓慢时分量趋于稳定，步长应该放大。
但是理论和实践表明，无法同时满足分量的这种需求，出现数值不稳定现象，误差急剧增加，无法继续求解，常微分方程的这个性质叫刚性。
刚性问题的解决关键在于步长的选择。

n阶微分方程
$$y^{(n)}=f(t,y,\dot{y},\dots ,y^{(n-1)})$$
$$y(0),\dot{y}(0),\dots ,y^{(n-1)}(0)$$
要降成一维
变量代换
$$\begin{array}{c}{x}_1=y\\ x_2=\dot{y}\\ \dots \\ x_n=y^{(n-1)}\end{array}$$
n个变量，整理成变量x的一阶方程组
$$\begin{array}{c}{\dot{x}}_1=x_2\\ {\dot{x}}_2=x_3\\ \dots \\ {\dot{x}}_n=f(t,x_1,x_2,\dots ,x_n)\end{array}$$
n个一阶方程

数值求解过程

方程转换形式
$$y^{\prime }(t)=f(t,y)$$
变量属性必须一一对应

function dy = myfun(t, y)

• dy，必须是列向量，长度为方程组的个数，通常与y的长度相同
• 如果是常微分方程组，y就是列向量
  函数中的输入参数和输出参数是形参，名字可以任意取，但必须满足上述条件。
  即输入参数有两个，第一个表示自变量，第二个是由因变量组成的列向量，输出参数必须是列向量。

Matlab软件求解时，需要建立两个文件：

1. 存放函数的.m文件
2. 计算的脚本文件
   一般的求解过程如下：
3. 建立m文件函数，文件名fun.m

function f = fun(t, x)
f = [x2(t); x3(t);... f(t, x1(t), x2(t), ..., xn(t)];

2. 数值计算（执行以下命令）写脚本文件

[t, x1, x2, ..., xn] = ode45('fun', [t0, tf], [x1(0), x2(0), ..., xn(0)])

微分方程数值解实例

求解Ver der Pol初值问题

注意：如果需求解的问题是高阶常微分方程，则需将其化为一阶常微分方程组，此时需用函数文件来定义该常微分方程组。
令
$$x_1=y,x_2=\frac{dy}{dt}$$
则原二阶微分方程可以转化成2个一阶微分方程组：
$$\{\begin{array}{c}\frac{d{x}_1}{dt}=x_2\\ \frac{d{x}_2}{dt}=7(1-x_1^2)x_2-x_1\\ x_1(0)=1,x_2(0)=0\end{array}$$

1. 先编写函数文件fun_test_1.m

function xprime = fun_test_1 (t, x)
xprime = [x(2); 7*(1 - x(1)^2) * x(2) - x(1)];

2. 再编写脚本文件test1.m，在命令窗口直接运行该文件。

clear;
y0 = [1; 0];
[t,x] = ode45('fun_test_1', [0, 40], y0); 
plot(t, x(:,1), 'r-');

• t，存放自变量的离散值
• x，存放变量代换之后的x1和x2
• ode45，45阶
• funtest1，存放微分方程组文件的文件名
• 0，40，自变量t的取值范围，ode45自动会在这个范围之内把t离散化
• y0，初值
  

Van der pol方程

令$y_1=x(t),y_2=x^{\prime }(t)$
$$\{\begin{array}{c}{y}_1^{\prime }=y_2\\ y_2^{\prime }=(1-y_1^2)y_2-y_1\end{array}$$
$$y_1(0)=3y_2(0)=0$$

1. 编写M文件（文件名为fun_test_2.m):

function dy = fun_test_2 (t, y);
	dy = zeros(2, 1); 
	dy(1) = y(2);
	dy(2) = (1 - y(1)^2) * y(2) - y(1);
%或用该表达式：
	dy = [y(2); (1 - y(1)^2) * y(2) - y(1)];

2. 编写脚本文件如下：（test2.m）

[t, y] = ode23('fun_test_2', [0, 20], [3, 0]);
	y1 = y(:, 1);     %原方程的解
	y2 = y(:, 2);
	plot(t, y1, t, y2, '--')    %y1(t),y2(t)曲线图
	pause,
	plot(y1, y2), grid    %相轨迹图，即y2(y1)曲线

计算结果
$$y_1=x(t)蓝色曲线$$
$$y_2=x^{\prime }(t)红色曲线$$

例3 Van der pol方程

令$y_1=x,y_2=y_1^{\prime }$
$$\{\begin{array}{c}{y}_1^{\prime }=y_2\\ y_2^{\prime }=1000(1-y_1^2)y_2-y_1\\ y_1(0)=2,y_2(0)=0\end{array}$$

1. 建立m文件fun_test_3.m如下：

function dy = fun_test_3(t, y)
	dy = zeros(2, 1); 
	dy(1) = y(2);
	dy(2) = 1000 * (1 - y(1)^2) * y(2) - y(1);

2. 取$t_0=0,t_f=3000$，输入命令（或编写脚本文件)

[T, Y] = ode15s('fun_test_3', [0 3000], [2 0]); 
plot(T, Y(:, 1), '-')

微分方程的初值敏感性

1. 编写M函数文件(lorenz.m)；
2. 编写脚本文件求解并画三维空间的相平面轨线；(Ltest.m)
   (1）编写M函数文件：

function xdot = lorenz(t, x)
xdot = [-8/3, 0, x(2); 0, -10, 10; -x(2), 28, -1] * x;

(2）编写脚本文件：Ltest.m

x0 = [0 0 0.1];
[t, x] = ode45('lorenz', [0,10], x0);
plot(t, x(:,1), '-', t, x(:,2), '*', t, x(:,3), '+’)
pause
plot3(x(:,1), x(:,2), x(:,3)), grid on

1. 初值1：
   $$x_0=[0,0.1,0.1];[t_0,t_f」=[0,30];$$
   解向量设为p
2. 初值2：
   $$x_0=[0.01,0.11,0.11];[t_0,t_f]=[0,30];$$
   解向量设为q
   2个初值同时演化，离散之间对应的值相差：
   $$p-q=(p_1-q_1,p_2-q_2,p_3-q_3)$$