目录
一、概念
1.1什么是微分方程建模
1.2使用场合
二、基于python的微分方程建模
2.1scipy.integrate.odeint() 函数
编辑2.2案例
编辑
三、基于MATLAB的微分方程建模
四、偏微分方程的求解
一、概念
1.1什么是微分方程建模
微分方程建模是数学模型的重要方法,因为许多实际问题的数学描述将导致求解微分方程的定解问题。所以微分方程求解就显得格外重要啦。
如何把实际问题转化为微分方程呢?
1.根据诗级要求确定要研究的晾(自变量、未知函数、必要的参数等)并且确定坐标系。
2.找出这些量所满足的基本规律(物理的、几何的、化学的或者生物学的)。
3.运用这些规律列出方程和定解条件。
列方程的方式
1.按规律直接列方程。
2.微元分析法与任意区域上取积分的方法。
3.模拟近似法。
1.2使用场合
1.2.1
例题
1.2.2例题2
1.2.3
1.2.4
二、基于python的微分方程建模
我们的选择是 Python 常用工具包三剑客:Scipy、Numpy 和 Matplotlib:
Scipy 是 Python 算法库和数学工具包,包括最优化、线性代数、积分、插值、特殊函数、傅里叶变换、信号和图像处理、常微分方程求解等模块。有人介绍 Scipy 就是 Python 语言的 Matlab,所以大部分数学建模问题都可以用它搞定。
Numpy 提供了高维数组的实现与计算的功能,如线性代数运算、傅里叶变换及随机数生成,另外还提供了与 C/C++ 等语言的集成工具。
Matplotlib 是可视化工具包,可以方便地绘制各种数据可视化图表,如折线图、散点图、直方图、条形图、箱形图、饼图、三维图,等等。
顺便说一句,还有一个 Python 符号运算工具包 SymPy,以解析方式求解积分、微分方程,也就是说给出的结果是微分方程的解析解表达式。很牛,但只能求解有解析解的微分方程,所以,你知道就可以了。
2.1scipy.integrate.odeint() 函数
SciPy 提供了两种方式求解常微分方程:基于 odeint
函数的 API 比较简单易学,基于 ode
类的面向对象的 API 更加灵活。
scipy.integrate.odeint(func, y0, t, args=(), Dfun=None, col_deriv=0, full_output=0, ml=None, mu=None, rtol=None, atol=None, tcrit=None, h0=0.0, hmax=0.0, hmin=0.0, ixpr=0, mxstep=0, mxhnil=0, mxordn=12, mxords=5, printmessg=0, tfirst=False)
odeint 的主要参数:
求解标准形式的微分方程(组)主要使用前三个参数:
2.2案例
# 1. 求解微分方程初值问题(scipy.integrate.odeint)
from scipy.integrate import odeint # 导入 scipy.integrate 模块
import numpy as np
import matplotlib.pyplot as plt
def dy_dt(y, t): # 定义函数 f(y,t)
return np.sin(t**2)
y0 = [1] # y0 = 1 也可以
t = np.arange(-10,10,0.01) # (start,stop,step)
y = odeint(dy_dt, y0, t) # 求解微分方程初值问题
# 绘图
plt.plot(t, y)
plt.title("scipy.integrate.odeint")
plt.show()
案例2:Scipy 求解一阶常微分方程组
# 2. 求解微分方程组初值问题(scipy.integrate.odeint)
from scipy.integrate import odeint # 导入 scipy.integrate 模块
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
# 导数函数, 求 W=[x,y,z] 点的导数 dW/dt
def lorenz(W,t,p,r,b): # by youcans
x, y, z = W # W=[x,y,z]
dx_dt = p*(y-x) # dx/dt = p*(y-x), p: sigma
dy_dt = x*(r-z) - y # dy/dt = x*(r-z)-y, r:rho
dz_dt = x*y - b*z # dz/dt = x*y - b*z, b;beta
return np.array([dx_dt,dy_dt,dz_dt])
t = np.arange(0, 30, 0.01) # 创建时间点 (start,stop,step)
paras = (10.0, 28.0, 3.0) # 设置 Lorenz 方程中的参数 (p,r,b)
# 调用ode对lorenz进行求解, 用两个不同的初始值 W1、W2 分别求解
W1 = (0.0, 1.00, 0.0) # 定义初值为 W1
track1 = odeint(lorenz, W1, t, args=(10.0, 28.0, 3.0)) # args 设置导数函数的参数
W2 = (0.0, 1.01, 0.0) # 定义初值为 W2
track2 = odeint(lorenz, W2, t, args=paras) # 通过 paras 传递导数函数的参数
# 绘图
fig = plt.figure()
ax = Axes3D(fig)
ax.plot(track1[:,0], track1[:,1], track1[:,2], color='magenta') # 绘制轨迹 1
ax.plot(track2[:,0], track2[:,1], track2[:,2], color='deepskyblue') # 绘制轨迹 2
ax.set_title("Lorenz attractor by scipy.integrate.odeint")
plt.show()
案例3:Scipy 求解高阶常微分方程
# 3. 求解二阶微分方程初值问题(scipy.integrate.odeint)
# Second ODE by scipy.integrate.odeint
from scipy.integrate import odeint # 导入 scipy.integrate 模块
import numpy as np
import matplotlib.pyplot as plt
# 导数函数,求 Y=[u,v] 点的导数 dY/dt
def deriv(Y, t, a, w):
u, v = Y # Y=[u,v]
dY_dt = [v, -2*a*v-w*w*u]
return dY_dt
t = np.arange(0, 20, 0.01) # 创建时间点 (start,stop,step)
# 设置导数函数中的参数 (a, w)
paras1 = (1, 0.6) # 过阻尼:a^2 - w^2 > 0
paras2 = (1, 1) # 临界阻尼:a^2 - w^2 = 0
paras3 = (0.3, 1) # 欠阻尼:a^2 - w^2 < 0
# 调用ode对进行求解, 用两个不同的初始值 W1、W2 分别求解
Y0 = (1.0, 0.0) # 定义初值为 Y0=[u0,v0]
Y1 = odeint(deriv, Y0, t, args=paras1) # args 设置导数函数的参数
Y2 = odeint(deriv, Y0, t, args=paras2) # args 设置导数函数的参数
Y3 = odeint(deriv, Y0, t, args=paras3) # args 设置导数函数的参数
# W2 = (0.0, 1.01, 0.0) # 定义初值为 W2
# track2 = odeint(lorenz, W2, t, args=paras) # 通过 paras 传递导数函数的参数
# 绘图
plt.plot(t, Y1[:, 0], 'r-', label='u1(t)')
plt.plot(t, Y2[:, 0], 'b-', label='u2(t)')
plt.plot(t, Y3[:, 0], 'g-', label='u3(t)')
plt.plot(t, Y1[:, 1], 'r:', label='v1(t)')
plt.plot(t, Y2[:, 1], 'b:', label='v2(t)')
plt.plot(t, Y3[:, 1], 'g:', label='v3(t)')
plt.axis([0, 20, -0.8, 1.2])
plt.legend(loc='best')
plt.title("Second ODE by scipy.integrate.odeint")
plt.show()
案例四求微分方程的通解
from sympy import *
x = symbols('x'); y = symbols('y',cls=Function)
eq1 = diff(y(x),x,2)-5*diff(y(x),x)+6*y(x)
eq2 = diff(y(x),x,2)-5*diff(y(x),x)+6*y(x)-x*exp(2*x)
print("齐次方程的解为:",dsolve(eq1,y(x)))
print("非齐次方程的解为:",dsolve(eq2,y(x)))
案例五
from sympy import *
x = symbols('x'); y = symbols('y',cls=Function)
eq1 = diff(y(x),x,2)-5*diff(y(x),x)+6*y(x)
eq2 = diff(y(x),x,2)-5*diff(y(x),x)+6*y(x)-x*exp(2*x)
print("初值问题的解为:{}".format(dsolve(eq1,y(x),ics={y(0):1,diff(y(x),x).subs(x,0):0})))
y2 = dsolve(eq2,y(x),ics={y(0):1,y(2):0})
print("边值问题的解为:{}".format(y2))
案例六
from sympy.abc import x
from sympy import Function, diff, dsolve, sin
y = Function('y')
eq = diff(y(x),x,2)+2*diff(y(x),x)+2*y(x)-sin(x) # 定义方程
con = {y(0): 0,diff(y(x), x).subs(x,0): 1} # 定义初值条件
y = dsolve(eq, ics=con)
print(y)
案例七求下述微分方程的解:
import sympy as sp
t = sp.symbols('t')
x1,x2,x3 = sp.symbols('x1,x2,x3', cls=sp.Function)
eq = [x1(t).diff(t)-2*x1(t)+3*x2(t)-3*x3(t),
x2(t).diff(t)-4*x1(t)+5*x2(t)-3*x3(t),
x3(t).diff(t)-4*x1(t)+4*x2(t)-2*x3(t),
]
con = {x1(0):1, x2(0):2, x3(0):3}
s = sp.dsolve(eq,ics=con)
print(s)
案例八
from matplotlib.font_manager import FontProperties
from scipy.integrate import odeint
from sympy.abc import t
import numpy as np
import matplotlib.pyplot as plt
my_font = FontProperties(fname=r"c:\windows\fonts\simsun.ttc", size=14)
def Pfun(y, x):
y1,y2 = y;
return np.array([y2, -2*y1-2*y2])
x = np.arange(0, 10, 1.0) # 创建时间点
sol1 = odeint(Pfun, [0.0, 1.0], x) # 求数值解
plt.plot(x, sol1[:, 0],'r*',label="数值解")
plt.plot(x, np.exp(-x)*np.sin(x),'g', label="符号解曲线")
plt.legend(prop=my_font) # 表示添加图例,并且用特有的 prop=my_font
plt.savefig("例八")
plt.show()
三、基于MATLAB的微分方程建模
人口模型
Matlab代码:
% 人口增长-指数增长模型
% 清空工作区和变量区
clear;clc;
% 定义迭代次数
n = 250;
% 定义并初始化参数
x = zeros(1,n);
x(1,1) = 5.42; % 初始化人口数(亿)
r = 0.018; % 人口增长率
% 开始迭代
for t = 2:n
x(1,t) = x(1,1)*exp(r*t);
end
% 绘图
plot(1:1:n,x);
legend('人口增长曲线');
xlabel('迭代次数');
ylabel('人口数(亿) ');
grid on;
% 人口增长-阻滞增长模型
% 清空工作区和变量区
clear;clc;
% 定义迭代次数
n = 500;
% 定义并初始化参数
x = zeros(1,n);
x(1,1) = 5.42; % 初始化人口数(亿)
r = 0.018; % 初始人口增长率
xm = 100; % 最大人口数
% 开始迭代
for t = 2:n
x(1,t) = xm/(1+(xm/x(1,1)-1)*exp(-r*t));
end
% 绘图
plot(1:1:n,x);
legend('人口增长曲线');
xlabel('迭代次数');
ylabel('人口数(亿) ');
grid on;
% 传染病SI模型
% 清空工作区和变量区
clear;clc;
% 定义迭代次数
n = 100;
% 定义并初始化状态E(s,i)
E = zeros(2,n);
E(1,1) = 0.99; % 初始健康者比例
E(2,1) = 0.01; % 初始感染者比例
% 初始化参数
L = 0.5; % 病人日接触率
% 开始迭代
for t = 1:n-1
E(2,t+1) = E(2,t) + L*E(2,t)*(1-E(2,t));
E(1,t+1) = 1 - E(2,t+1);
end
% 绘图
s = E(1,:); % 健康者比例数据
i = E(2,:); % 感染者比例数据
plot(s,'DisplayName','s');hold on;
plot(i,'DisplayName','i');
legend('健康者比例','感染者比例');
xlabel('迭代次数');
ylabel('比例');
grid on;
hold off;
SI模型的缺点在于,没有考虑病人可以被治愈的情况,导致最后所有健康者都会变为感染者,不太符合实际情况
% 传染病SIS模型
% 清空工作区和变量区
clear;clc;
% 定义迭代次数
n = 100;
% 定义并初始化状态E(s,i)
E = zeros(2,n);
E(1,1) = 0.99; % 初始健康者比例
E(2,1) = 0.01; % 初始感染者比例
% 初始化参数
L = 0.5; % 病人日接触率
M = 0.2; % 病人日治愈率
% 开始迭代
for t = 1:n-1
E(2,t+1) = E(2,t) + L*E(2,t)*(1-E(2,t))-M*E(2,t);
E(1,t+1) = 1 - E(2,t+1);
end
% 绘图
s = E(1,:); % 健康者比例数据
i = E(2,:); % 感染者比例数据
plot(s,'DisplayName','s');hold on;
plot(i,'DisplayName','i');
legend('健康者比例','感染者比例');
xlabel('迭代次数');
ylabel('比例');
grid on;
hold off;
SIS模型作为SI模型的升级版,考虑了病人被治愈的情况,但同时,SIS模型没有考虑病人被治愈后,可能会产生抗体,而不会被再次感染的情况,不太符合某些传染病的模拟
% 传染病SIR模型
% 清空工作区和变量区
clear;clc;
% 定义迭代次数
n = 100;
% 定义并初始化状态E(s,i,r)
E = zeros(3,n);
E(1,1) = 0.99; % 初始健康者比例
E(2,1) = 0.01; % 初始感染者比例
E(3,1) = 1 - E(1,1) - E(2,1); % 初始免疫者比例(假设一开始没有免疫者)
% 初始化参数
L = 0.5; % 病人日接触率
M = 0.1; % 日治愈率
% 开始迭代
for t = 1:n-1
E(1,t+1) = E(1,t) - L*E(1,t)*E(2,t);
E(2,t+1) = E(2,t) + L*E(1,t)*E(2,t) - M*E(2,t);
E(3,t+1) = 1 - E(1,t+1) - E(2,t+1);
end
% 绘图
s = E(1,:); % 健康者比例数据
i = E(2,:); % 感染者比例数据
r = E(3,:); % 免疫者比例数据
plot(s,'DisplayName','s');hold on;
plot(i,'DisplayName','i');
plot(r,'DisplayName','r');
legend('健康者比例','感染者比例','免疫者比例');
xlabel('迭代次数');
ylabel('比例');
grid on;
hold off;
可以看到,在SIR模型中,随着时间的推移,感染者人数先上升,达到一个峰值后,再下降,最后下降为0,所有人都成为免疫者
四、偏微分方程的求解
例题
例题
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
import numpy as np
plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']
def fun(t, w):
x = w[0]
y = w[1]
return [-x**3-y,-y**3+x]
# 初始条件
y0 = [1,0.5]
yy = solve_ivp(fun, (0,100), y0, method='RK45',t_eval = np.arange(0,100,1) )
t = yy.t
data = yy.y
plt.plot(t, data[0, :])
plt.plot(t, data[1, :])
plt.xlabel("时间s")
plt.show()