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

ChoSeitaku

于 2024-08-22 13:18:35 发布

阅读量663

点赞数 12

分类专栏：数学建模文章标签：数学建模

本文链接：https://blog.csdn.net/CoderZzz6310/article/details/141425605

版权

数学建模专栏收录该内容

14 篇文章 0 订阅

订阅专栏

微分方程简介

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

微分方程的数学描述

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

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

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

微分方程的解

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

微分方程的初值问题

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

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

微分方程建模基本方法

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

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

微分方程建模基础

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

微分方程建模准则

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

微分方程建模过程

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

凶案时间推断

![[Pasted image 20240822080831.png]]

问题分析
该问题符合物理学的冷却现象，可以应用Newton冷却定律：“物体在介质中的冷却速度同该物体温度与介质温度之差成正比”来解决。
用速度来描述物体在某时刻内的变化率，涉及到导数概念，反映在数学模型上就可以运用微分方程来建模
模型建立
按一般情形考虑，记 $T_{t}$ 为时刻 $t$ 物体的温度， $T_{0}$ 为初始时刻 $t_{0}$ 物体的温度（本例中为受害人被害时的体温）， $T_{e}$ 为介质的温度（环境温度），则由Newton冷却定律可得一阶线性微分方程模型
$\left\{\begin{matrix} \frac{dT_{t}}{dt}=-\lambda(T_{t}-T_{e}) \\ T_{t_{0}}=T_{0} \end{matrix}\right.$
其中 $\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^{''} - 2 y^{'} - 3 y = 0$
$y (0) = 1, y^{'} (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 二阶非常系数微分方程
$\left\{\begin{matrix} x''(t)-(1-x^{2}(t))x'(t)+x(t)=0 \\ x(0)=3,\ x'(0)=0 \end{matrix}\right.$
输入

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
![[Pasted image 20240822100039.png]]

也可以写成：

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
求微分方程组的通解
$\left\{\begin{matrix} \frac{dx}{dt}=2x-3y+3z \\ \frac{dy}{dt}=4x-5y+3z \\ \frac{dz}{dt}=4x-4y+2z \end{matrix}\right.$
输入命令

[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)

微分方程数值解原理

常微分方程数值解的意义

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

建立数值解的一般方法

![[Pasted image 20240822104906.png]]

用差商代替导数
若步长h较小，则有如下变换
$y'(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}$
可得到公式
$\left\{\begin{matrix} y_{i+1}=y_{i}+hf(x_{i},y_{i}) \\ y_{0}=y(x_{0}) \end{matrix}\right. ,\ i=0,1,2,\dots,n-1$
欧拉公式
$\left\{\begin{matrix} \frac{dy}{dx}=y+ \frac{2x}{y^{2}} \\ y(0)=1 \end{matrix}\right.$
解析解
$y=\left( \frac{5}{3}e^{3x}-2x- \frac{2}{3} \right)^{1/3}$

使用数值积分方法

对方程 $y^{'} = 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}))]$
可以得到
$\left\{\begin{matrix} 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{matrix}\right.$
实际应用时与欧拉公式结合使用
$\left\{\begin{matrix} 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{matrix}\right.$
对给定的精度 $\varepsilon$ ，当满足
$y_{i+1}^{(k+1)}-y_{i+1}^{(k)}|<\varepsilon$
取
$y_{i+1}=y_{i+1}^{(k+1)}$
然后继续下一步 $y_{i+2}$ 的计算，这就是改进的欧拉方法

使用泰勒公式

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

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

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

![[Pasted image 20240822110919.png]]

微分方程数值解的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求解器简要说明

![[Pasted image 20240822112913.png]]

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

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

n阶微分方程
$y^{(n)}=f(t,y,\dot{y},\dots,y^{(n-1)})$
$y(0),\dot{y}(0),\dots,y^{(n-1)}(0)$
要降成一维
变量代换
$\begin{array}{} x_{1}=y \\ x_{2}=\dot{y} \\ \dots \\ x_{n}=y^{(n-1)} \end{array}$
n个变量，整理成变量x的一阶方程组
$\begin{array}{} \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^{'} (t) = f (t, y)$
变量属性必须一一对应

function dy = myfun(t, y)

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

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

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

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

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

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

微分方程数值解实例

求解Ver der Pol初值问题

![[Pasted image 20240822120559.png]]

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

先编写函数文件fun_test_1.m

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

再编写脚本文件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方程

![[Pasted image 20240822124354.png]]

令 $y_{1}=x(t),\ y_{2}=x'(t)$
$\left\{\begin{matrix} y_{1}'=y_{2} \\ y_{2}'=(1-y_{1}^{2})y_{2}-y_{1} \end{matrix}\right.$
$y_{1}(0)=3\qquad y_{2}(0)=0$

编写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)];

编写脚本文件如下：（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)\qquad蓝色曲线$
$y_{2}=x'(t)\qquad 红色曲线$

![[Pasted image 20240822125153.png]]
![[Pasted image 20240822125339.png]]

例3 Van der pol方程

![[Pasted image 20240822125415.png]]

令 $y_{1}=x,\ y_{2}=y_{1}'$
$\left\{\begin{matrix} y_{1}'=y_{2} \\ y_{2}'=1000(1-y_{1}^{2})y_{2}-y_{1} \\ y_{1}(0)=2,\ y_{2}(0)=0 \end{matrix}\right.$

建立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);

取 $t_{0}=0,t_{f}=3000$ ，输入命令（或编写脚本文件)

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

![[Pasted image 20240822130250.png]]

微分方程的初值敏感性

![[Pasted image 20240822130321.png]]

编写M函数文件(lorenz.m)；
编写脚本文件求解并画三维空间的相平面轨线；(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

![[Pasted image 20240822130648.png]]

初值1：
$x_{0} = [0, 0.1, 0.1]; \ [t_{0},t_{f}」= [0,30];$
解向量设为p
初值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})$

ChoSeitaku

关注

12
点赞
踩
5

收藏

觉得还不错? 一键收藏
0
评论
微分方程模型|凶案时间推断|欧拉法|Vanderpol方程

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

扫一扫