Python的Scipy库解微分方程

微分方程：

初始值：

问题：

求解其他三个参数：

代码实现：

import numpy as np
from numpy import zeros, linspace, arange
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3Ddef read_data(path):with open(path, 'r') as f:data = f.read()time, num = [], []data = data.strip().split('\n')for items in data:t, n = items.split()time.append(float(t))num.append(float(n))return time, num# dS/dt
def funcS(S, R, A, r, Nmax, ds, beta, Emax, MICs, H):Es = 1 - Emax * np.power(A, H) / (np.power(MICs, H) + np.power(A, H))return r * (1 - (R + S) / Nmax) * Es * S - ds * S - beta * S * R / (S + R)# dR/dt
def funcR(S, R, A, r, Nmax, dr, beta, Emax, MICr, H, gamma):Er = 1 - Emax * np.power(A, H) / (np.power(MICr, H) + np.power(A, H))return r * (1 - gamma) * (1 - (R + S) / Nmax) * Er * R - dr * R - beta * S * R / (S + R)# dA/dt
def funcA(A, da):return -da * A# function that need some parameters
def func_all(y, t, parameters):S, R, A = yr, Nmax, ds, dr, Emax, H, MICs, MICr, beta, gamma, da = parametersfunc1 = funcS(S, R, A, r, Nmax, ds, beta, Emax, MICs, H)func2 = funcR(S, R, A, r, Nmax, dr, beta, Emax, MICr, H, gamma)func3 = funcA(A, da)return np.array([func1, func2, func3])path = './amrdata.txt'
time, num = read_data(path)
print(time)
print(num)# set up a parameters
r = 0.5
Nmax = 1.0e7
ds = 0.025
dr = 0.025
Emax = 2
H = 2
MICs = 8
MICr = 2000# guess
beta = 1.0e-11
gamma = 0.05
da = 0.001# set up initial conditions
S0 = 9 * 1.0e6
R0 = 1.0e5
A0 = 5.6
y0 = [S0, R0, A0]# beta_size = 10
# beta_range = []
# for i in range(beta_size):
#     beta_range.append(beta)
#     beta *= 10
# print(beta_range)
bmin, bmax = 1.0e-6, 1.0e-4
gmin, gmax = 0.01, 0.2
amin, amax = 0.001, 0.03# grid
beta_range = arange(bmin, bmax, (bmax - bmin) / 30)
print(beta_range)
gamma_range = arange(gmin, gmax, (gmax - gmin) / 30)
print(gamma_range)
da_range = arange(amin, amax, (amax - amin) / 30)
print(da_range)MIN_RMSE = 1000# find parameters
best_beta = beta
best_gamma = gamma
best_da = da# grid search
for b in beta_range:for g in gamma_range:for a in da_range:beta = bgamma = gda = aparameters = [r, Nmax, ds, dr, Emax, H, MICs, MICr, beta, gamma, da]# Solve the equations using built-in SciPy ODE solver odeintsol = odeint(func_all, y0, time, (parameters,))S, R, A = sol[:, 0], sol[:, 1], sol[:, 2]  # sensitive bacteria, resistant bacteria and antibioticN = 100 * R / (S + R)  # equal to 100*R/(S + R)res = 0for i in range(len(num)):res += (num[i] - N[i]) ** 2RMSE = np.sqrt(res) / len(num)# print(RMSE)# print("beta= %f, gamma= %f, da= %f" % (beta, gamma, da))if RMSE < MIN_RMSE:MIN_RMSE = RMSEbest_beta = betabest_gamma = gammabest_da = daprint("---------MIN_RMSE-------\n RMSE= %f" % MIN_RMSE)print("beta= %f, gamma= %f, da= %f" % (beta, gamma, da))parameters = [r, Nmax, ds, dr, Emax, H, MICs, MICr, best_beta, best_gamma, best_da]
sol = odeint(func_all, y0, time, (parameters,))
S, R, A = sol[:, 0], sol[:, 1], sol[:, 2]  # sensitive bacteria, resistant bacteria and antibiotic
N = 100 * R / (S + R)  # equal to 100*R/(S + R)# fig = plt.figure()
# ax = Axes3D(fig)
# print(sol)
# ax.plot(sol[:, 0], sol[:, 1], sol[:, 2])
# ax.set_xlabel("S concentration ")
# ax.set_ylabel("R concentration ")
# ax.set_zlabel("A concentration ")plt.plot(time, num)
plt.plot(time, N)
plt.xlabel("time/H")
plt.ylabel("the number of resistant strains")
plt.show()

拟合情况：

可以看到两条曲线很接近，效果很好。

Python的Scipy库解微分方程相关推荐

Python学习-Scipy库稀疏矩阵的建立（面向列的稀疏矩阵、基于坐标格式的稀疏矩阵）
Python学习-Scipy库稀疏矩阵的建立稀疏矩阵指在矩阵中值为0的元素的数量远远多于非0值的矩阵 (非0元素总数/所有元素总数<=0.05) 稀疏矩阵的实现对象: csc_matrix() ...
Python学习-Scipy库信号处理signal(过滤、快速傅里叶变换、信号窗函数、卷积)
Python学习-Scipy库信号处理signal 目录 1.过滤:以某种方式修改输入信号 2.快速傅里叶变换 3.信号窗函数 4.卷积导入库 import matplotlib.pyplot as ...
python导入scipy库、sympy库遇到的问题及解决方式
首先从cmd中导入scipy库,输入代码: pip install scipy 注意: pip版本最好也要更新到最新版,否则容易发生版本冲突的问题. 但是出现异常:read time out 这时想到 ...
Python中scipy库对mat文件进行读写操作
mat文件是以字典的格式进行存储的,有时候Python中需要对字典进行读写,使用Python处理matlab的mat文件时,可以使用scipy库中的函数进行操作. 导入scipy库对mat文件的读写 ...
python安装scipy库_Numpy与Scipy的安装
Python下大多数工具包的安装都很简单,只需要执行 "python setup.py install"命令即可.然而,由于SciPy和numpy这两个科学计算包的依赖关系较多,安 ...
Python中scipy库中csr_matrix()函数和csc_matrix()函数的解释
在使用Python进行科学计算时经常需要用到稀疏矩阵的构造,而python的科学计算包scipy.sparse是很好的一个解决稀疏矩阵构造/计算的包. 构造稀疏矩阵常用的两个函数为:csr_matri ...
python概率论_概率论中常见分布总结以及python的scipy库使用
概率分布有两种类型:离散(discrete)概率分布和连续(continuous)概率分布. 离散概率分布也称为概率质量函数(probability mass function).离散概率分布的例子有 ...
scipy 概率泊松分布_概率论中常见分布总结以及python的scipy库使用：两点分布、二项分布、几何分布、泊松分布、均匀分布、指数分布、正态分布......
概率分布有两种类型:离散(discrete)概率分布和连续(continuous)概率分布. 离散概率分布也称为概率质量函数(probability mass function).离散概率分布的例子有 ...
概率论中常见分布总结以及python的scipy库使用：两点分布、二项分布、几何分布、泊松分布、均匀分布、指数分布、正态分布...
概率分布有两种类型:离散(discrete)概率分布和连续(continuous)概率分布. 离散概率分布也称为概率质量函数(probability mass function).离散概率分布的例子有 ...
python参数估计_python简单实现最大似然估计scipy库的使用详解
python简单实现最大似然估计 1.scipy库的安装 wim+R输入cmd,然后cd到python的pip路径,即安装:pip install scipy即可 2.导入scipy库 from sc ...

Python的Scipy库解微分方程

微分方程：

初始值：

问题：

代码实现：

拟合情况：

Python的Scipy库解微分方程相关推荐

最新文章

热门文章