数学建模-模型详解（2）

微分模型

当谈到微分模型时，通常指的是使用微分方程来描述某个系统的动态行为。微分方程是描述变量之间变化率的数学方程。微分模型可以用于解决各种实际问题，例如物理学、工程学、生物学等领域。

微分模型可以分为两类：常微分方程和偏微分方程。常微分方程描述的是只有一个自变量的函数的导数，而偏微分方程描述的是多个自变量的函数的导数。

常微分方程的一些常见类型包括：一阶线性常微分方程、一阶非线性常微分方程、高阶线性常微分方程等。偏微分方程的一些常见类型包括：热传导方程、波动方程、扩散方程等。

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.032 , x m = 159.3 r=0.032,x_m=159.3r=0.032,x =159.3，并预测 1900 年人口为 73.09 73.0973.09，得到预测曲线。

Python 代码
微分方程求解
求微分方程:


# %%

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