【数学建模笔记】数学建模的微分方程建模

微分方程建模是数学建模的重要方法，大体可以按以下几步：

1. 根据实际要求确定要研究的量 (自变量、未知函数、必要参数)，确定坐标系；
2. 找出这些量所满足的基本规律；
3. 运用规律列出方程和定解条件。

微分方程的数值解

考虑一阶常微分方程 在区间 上的解。

所谓数值解法，就是求 在若干点 处的近似值 ，称 为步长，一般取等步长 。

为求数值解，首先要将微分方程离散化，方法有：

• 用差商近似导数：

若用向前差商 代替 ，代入微分方程，则 即 将 的近似值 代入，即 于是原方程组化为 通过 可逐次算出 。

• 数值积分方法：

右边的积分用矩形或梯形公式计算。

• 泰勒多项式近似：

将 在 处展开，即

代入近似值得 上式又称欧拉法，为了提高精度，可以使用梯形积分代替矩形积分，称改进欧拉法，即 然而由于每一步都会产生截断误差，欧拉法会随着递推误差越来越大。

龙格-库塔法

龙格库塔法的主要思想，就是在 点附近选取一些特定的点，将这些点的函数值进行线性组合，用组合值代替泰勒展开中的导数值。

根据微分中值定理有 又 ，于是有 记 为区间 的平均斜率。

欧拉法简单取 ​​，精度很低，而改进欧拉法则取

则提高了精度。

于是，在区间 内多取几个点进行加权平均作为 ，就有可能构造出精度更高的计算公式。

不妨在区间内仍取 2 个点，有 又 ，则可以化为 在 作泰勒展开，则有

于是有

可见为使误差 ，只须令 解不唯一，其中若 即为改进欧拉法。

若取 4 个点，就构成常用的四阶龙格-库塔法 RK4 其中

​ 分别为起点、中点、终点的斜率，通过一定的加权平均使得误差尽可能小。

微分方程建模实例

Malthus 模型

设 表示 时人口数，人口增长率 是常数，假设人口数量的变化封闭，即人口数量增加和减少只取决于人口个体的生育和死亡。

则 时刻到 时刻，人口增量为 于是得 解为 然而模型的预测结果远高于实际人口增长，误差的原因是对增长率 估计过高。

Logistic 模型

由于资源是有限的，随着人口数量增长，资源对人口增长的限制将越来越显著。人口较少时可以将增长率 看作常数，到达一定数量后 应随人口增加而减小。因此，将 表示为人口 的函数 ，且 为 的减函数。

设 为 的线性函数，。

设容纳最大人口为 ，当 时，。

于是有 得 解为 以 为变量，对某一数据集进行拟合，就可以得到结果。

Python 代码

微分方程求解

求微分方程 的解析解和数值解，代码如下：

#! /usr/bin/env python
# -*- coding: utf-8 -*-
# @ author: Koorye
# @ date: 2021-7-27
# @ function: 使用 scipy、sympy 求微分方程的数值解、解析解

# %%

import numpy as np
from scipy.integrate import odeint
from sympy import *

# %%

# 使用 scipy 求数值解
# 微分方程
dy = lambda y,x:-2*y + x**2 + 2*x

# 数值范围
x1 = np.linspace(1,10,20)

# 求数值解，y 的初始值为 2
y1 = odeint(dy, 2, x1)
y1

# %%

# 使用 sympy 求解析解

# 定义变量和函数
x = symbols('x', real=True)
y = Function('y')

# 定义方程和约束
eq = y(x).diff(x) + 2*y(x) - x**2 - 2*x
con = {
	y(1): 2,
}

# 求解
f = simplify(dsolve(eq, ics=con))
f

# %%

# 向解中代入不同的值
x2 = np.linspace(1,10,100)
y2 = []
for each in x2:
	y2.append(list(sorted(f.subs(x,each).evalf().atoms()))[1])

# %%

# 画图
import matplotlib.pyplot as plt

plt.scatter(x1,y1, label='x1', color='coral')
plt.plot(x2,y2, label='x2')
plt.legend()

输出如下：

得到解析解的公式，以及数值解 和解析解的曲线 ，发现数值解和解析解大致吻合。

Logistic 模型

对下列人口数据：

年份	人口
1790	3.9
1800	5.3
1810	7.2
1820	9.6
1830	12.9
1840	17.1
1850	23.2
1860	31.4
1870	38.6

使用 Logistic 模型拟合，求出 的值。

代码如下：

#! /usr/bin/env python
# -*- coding: utf-8 -*-
# @ author: Koorye
# @ date: 2021-7-27
# @ function: 使用 Logistic 模型拟合人口数据

# %%

import numpy as np
import pandas as pd
from scipy.optimize import curve_fit

# %%

# 源数据
df = pd.DataFrame({
    'year': [1790, 1800, 1810, 1820, 1830, 1840, 1850, 1860, 1870],
    'population': [3.9, 5.3, 7.2, 9.6, 12.9, 17.1, 23.2, 31.4, 38.6],
})

x0 = float(df['population'][0])
t0 = float(df['year'][0])

# %%

# Logistic 模型
def x(t, r, xm):
    return xm / (1 + (xm/x0-1)*np.exp(-r*(t-t0)))

# 拟合参数
popt, pcov = curve_fit(x,
                       df['year'].tolist(),
                       df['population'].tolist(),
                       bounds=((0, 1), (.1, np.inf)))
r, xm = popt[0], popt[1]
print('r =', r)
print('xm =', xm)

# 预测 1900 人口
print('population in 1900 =', x(1900, r, xm))

# %%

# 画出预测曲线
import matplotlib.pyplot as plt

year = np.linspace(1790,2000,21)
population = []
for each in year:
	population.append(x(each,r,xm))
plt.scatter(df['year'], df['population'], label='actual')
plt.plot(year, population, label='predict', color='coral')
plt.legend()

输出如下：

r = 0.031999710099329476
xm = 159.3028433863908
population in 1900 = 73.09203839061287

得到 ，并预测 1900 年人口为 ，得到预测曲线。