微分方程模型

代码部分可见Sympy

初值问题

这类问题一般是知道一个函数的导数解，并该函数在某一点的求解，且无法直接解出该函数式，只能通过数值方法来代替。

[注]: y是一个单独的函数，f(x,y)为x和y(x)的函数

欧拉差商法

1. 向前差商

2. 向后差商

3. 梯形公式

4. 改进的欧拉公式

局部截断误差估计

• 局部截断误差
• p阶精度

其中，

• 向前欧拉公式的精度为1阶
• 向后欧拉公式的精度为1阶
• 梯形公式的局部截断误差的精度为2阶
• 改进欧拉公式的局部截断误差的精度为2阶

龙格—库塔方法

改进的欧拉公式可以写成

• 二阶龙格－库塔公式

• 三阶龙格－库塔公式

• 四阶龙格－库塔公式

Scipy求解

• 通用函数scipy.integrate.solve_ivp，这里有一个翻译版本

  scipy.integrate.solve_ivp（fun，t_span，y0，method = 'RK45'，t_eval = None，density_output = False，events = None，vectorized = False）

  对于该函数，公式如下
  dy / dt = f(t, y)
  y(t0) = y0

  参数：

  • fun: 要求解的函数
  • t_span: 积分间隔(t0，tf),求解器从t = t0开始并进行积分，直到达到t = tf
  • y0: y的初始值
  • method:
    对于非刚性问题，应使用显式Runge-Kutta方法(‘RK23’，‘RK45’，‘DOP853’)，对于刚性问题应使用隐式方法(‘Radau’，‘BDF’)。在Runge-Kutta方法中，建议使用‘DOP853’进行高精度求解(rtol和atol的值较低)。如果不确定，请首先尝试运行‘RK45’。如果异常地进行了许多次迭代，发散或失败，则您的问题可能很棘手，应使用‘Radau’或‘BDF’。 ‘LSODA’也可以是一个很好的通用选择，但由于它包装了旧的Fortran代码，因此使用起来可能不太方便。
  • t_eval: 指定求解的点，若无则用求解器分的点 返回值：
  • t： 取的点
  • y： 求解值
  • sol: 求解器，可以使用return.sol(t)来获得t点处的解值

例题

1. 简单的微分初值问题

from scipy.integrate import solve_ivp
import numpy as np
import matplotlib.pyplot as plt
%matplotlib widget

def fun(x,y):
    return -2*x*y

t_ev1 = np.linspace(0, 1.8, 100)

fig = plt.figure()
ax1 = fig.add_subplot(1,1,1)
plt.xlim([0,1.8])

sol2 = solve_ivp(fun, t_span=[0,1.8], y0=[1], method='RK23', t_eval=t_ev1)
plt.plot(sol2.t, sol2.y[0],'r-', label='RK23')
plt.legend(loc = 1)

sol1 = solve_ivp(fun, t_span=[0,1.8], y0=[1], method='RK45', t_eval=t_ev1)
plt.plot(sol1.t, sol1.y[0], 'c.',label='RK45', markersize=3)
plt.legend(loc = 1)

绘图如下

2. Van der Pol
   Van der pol又称范德波尔振荡器，来源于物理学中的一个公式，这里给出简单的公式并进行求解 可将以上公式转化成以下两个公式:

#y1为y的1次导数，返回[y2, y1]为求解微分值的列表
def van_del_pol(y, z):
    y, y1 = z
    return [y1, (1-y**2)*y1 - y]

# 为了便于观察只取100个点，太密了不好观察
t_ev1 = np.linspace(0,20,100)

sol2 = solve_ivp(van_del_pol, t_span=[0,20], y0=[3,0], method='RK23', t_eval=t_ev1)
plt.plot(sol1.t, sol1.y[0], 'c-',label='RK23', markersize=2)
plt.legend(loc = 1)

sol1 = solve_ivp(van_del_pol, t_span=[0,20], y0=[3,0], method='RK45', t_eval=t_ev1)
plt.plot(sol1.t, sol1.y[0], 'ro',label='RK45', markersize=2)
plt.legend(loc = 1)

3. 也可以使用odenit求解，这个其实比较老，官方建议用上面那个函数

from scipy.integrate import odeint
def Pfun(y, x):
    y1, y2 = y
    return [y2, -2*y1-2*y2]
x = np.arange(1,10,0.1)
sol = odeint(Pfun, [0,1], x)

Lorenz模型的混沌效应

Lorenz简介，该模型已经成混沌领域的经典模型，方程为

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# 图片大小
WIDTH, HEIGHT, DPI = 1000, 750, 100

# Lorenz参数 和 初值
sigma, beta, rho = 10, 2.667, 28
u0, v0, w0 = 0, 1, 1.05

# 区间0~tmax 和 点的数量n
tmax, n = 100, 10000

def lorenz(t, X, sigma, beta, rho):
    """The Lorenz equations."""
    u, v, w = X
    up = -sigma*(u - v)
    vp = rho*u - v - u*w
    wp = -beta*w + u*v
    return up, vp, wp

# Integrate the Lorenz equations.
soln = solve_ivp(lorenz, (0, tmax), (u0, v0, w0), args=(sigma, beta, rho),
                 dense_output=True)

# t为取值点
t = np.linspace(0, tmax, n)

# sol(t)为函数在t点的解值
x, y, z = soln.sol(t)

# 绘制三维 Lorenz吸引子(attractor)
fig = plt.figure(facecolor='w', figsize=(WIDTH/DPI, HEIGHT/DPI))
ax = fig.gca(projection='3d')
ax.set_facecolor('w')
fig.subplots_adjust(left=0, right=1, bottom=0, top=1)

# 使用不同颜色绘制
s = 10
cmap = plt.cm.winter
for i in range(0,n,s):
    ax.plot(x[i:i+s+1], y[i:i+s+1], z[i:i+s+1], color=cmap(i/n), alpha=0.4)

# 移除所有的坐标轴，只留下曲线
ax.set_axis_off()
plt.savefig('lorenz.png', dpi=DPI)
plt.show()

混沌系统对于不同的初值变化是否敏感，下面为一个初值仅差0.0001的两个解的偏差曲线(右)

def lorenz(w, t):
    sigma, rho, beta = 10, 28, 8/3
    x, y, z = w
    return np.array([sigma*(y-x), rho*x-y-x*z, x*y-beta*z])
    t = np.arange(0, 50, 0.01)
    
sol1 = odeint(lorenz, [0, 1, 0], t)
sol2 = odeint(lorenz, [0, 1.0001, 0], t)
ax1 = plt.subplot(121, projection='3d')
ax1.plot(sol1[:, 0], sol1[:, 1], sol1[:, 2], 'r')
ax2 = plt.subplot(122, projection='3d')
ax2.plot(sol2[:, 0] - sol1[:, 0], sol2[:, 1] - sol1[:, 1], sol2[:, 2] - sol1[:, 2], 'g')
plt.show()

掠食者—猎物方程

Lotka-Volterra equation是一个二元一阶非线性微分方程组。经常用来描述生物系统中，掠食者与猎物进行互动时的动态模型，也就是两者族群规模的消长

• t 是时间
• y 是掠食者的数量(狼)
• x 是猎物的数量(兔子)
• dy/dt 与 dx/dt 表示上述两族群相互对抗的时间之变化
• α, β, γ 与 δ 表示与两物种互动有关的系数，皆为正实数
• 第一式所表达的是猎物族群的增值速度 此模型假设猎物所接受的食物供给已经达到最极限，且除非遭遇掠食者的捕食，否则繁殖数量的增加以指数方式成长，其指数成长的情形，则以上述方程中的 αx 表现。此外并假设猎物遭遇捕食的比例，和猎物遭遇掠食者的机会成常数比，以上述方程中的 βxy 表现。如果 x 或 y 其中一个为零，则皆有可能是没有捕食行为出现。
  由上述的方程可知：猎物族群规模的改变，源于本身受到捕食而产生的成长衰减。
• 第二式所表达的是掠食者族群的增值速度 此方程中的 δxy 表示掠食者族群的成长（可能会与掠食者与猎物的数量比例相似，但是掠食者与猎物的数量比例是以不同的常数表示，且不一定与族群的成长相等。） γy 表示掠食者的自然死亡，为指数衰减。 由上述的方程可知：掠食者族群规模的改变，是猎食者族群的成长，减去其自然死亡的部分。

下面是官方给的一个例子

def lotkavolterra(t, z, a, b, c, d):
    x, y = z
    return [a*x - b*x*y, -c*y + d*x*y]

sol = solve_ivp(lotkavolterra, [0, 15], [10, 5], args=(1.5, 1, 3, 1),
                dense_output=True)
t = np.linspace(0, 15, 300)
z = sol.sol(t)
plt.plot(t, z.T)
plt.xlabel('t')
plt.legend(['x', 'y'], shadow=True)
plt.title('Lotka-Volterra System')
plt.show()

指数增长模型——Malthus 模型

1. 模型假设
   • 设x(t)表示t时刻的人口数，且x(t)连续可微
   • 人口增长率r为常数(r = 出生率 - 死亡率)
   • 人口数量增减取决于个体的生育和死亡，且每个个体生育能力相同
2. 建模与求解
   解为$ x(t)=x_0e^{rt} $
3. 模型评价 该模型对人口数量没有限制，长时间后人口数量将过高。

阻滞增长模型——Logistic 模型

地球所能容纳的人口数量是有限的，当人口达到一定数量后，增长率r应该为0，因此r为一个关于x的减函数，该模型使用了Logistic函数

1. 模型假设
   • 设r(x)为x的线性函数，$r(x)=r-sx$(首先使用线性)
   • 设最大人口数为$x_m$，当$x=x_m$时，增长率$r(x_m)=0$
2. 建模与求解 根据上述假设，得$r(x)=r(1-\frac{x}{x_m})$
   上式为一个可分离变量的方程，解得
   
   然后可以使用该函数对已知数据进行拟合来预测未来人口数量

传染病模型

1. 模型假设

   • 符号约定

     符号	意义
     n	总人数
     s(t)	t时刻易感染者数量
     i(t)	t时刻感染者数量
     r(t)	t时刻免疫者的数量

   • 易感人数的变化率 和 感染人数成正比，比例系数为$\lambda$
   • 免疫着人数变化率 和 感染人数成正比，比例系数为$\mu$
   • 每个人都为这三类人数中的一种

2. 根据以上假设建立如下模型
   该模型又称Kermack-Mckendrick方程