和欧拉用 python 养鱼

小小詹同学

发布于 2018-09-25 09:55:59

78900

代码可运行

文章被收录于专栏：小詹同学小詹同学

运行总次数：0

代码可运行

宠物鱼店主戴夫最近引进了一批彩虹鱼，打算捞一笔，他引进了30条彩虹鱼，每天要卖出20条，未来会发生什么呢？让我们和欧拉预测一波

已知

20条鱼每天能生殖14条小鱼
每条鱼每天平均产出0.7条鱼
照料的很好，死亡率 = 0
简化模型，刚出生的小鱼就能生殖小鱼
鱼缸最大环境承受值750条
越接近环境承受值，小鱼的生育力越低（暗含死亡率）
每天卖出20条小鱼

建模

P(0) = 30 第0天的小鱼数量 P_init = 30
P(1) = P(0) + 0.7 * P(0) - 20 第1天的小鱼的数量
P(t+Δt) = P(t) +0.7 * P(t) * Δt -20 * Δt 第t+Δt天相对第t天的小鱼数量
P(t+Δt) = P(t) +0.7 *(1 - P(t)/750) * P(t) * Δt -20 * Δt 加入环境承受值对于生育率的影响(简单模型)

简单的分析一下P的变化-微分和导数

这是一个P的导数，相关与P函数本身的一个微分方程，Autonomous differential equations 自控微分方程。看上去是不是很复杂，这个时候我们就要呼唤欧拉了：欧拉方法，命名自它的发明者莱昂哈德·欧拉()，是一种一阶数值方法，用以对给定初值的常微分方程（即初值问题）求解。它是一种解决数值常微分方程的最基本的一类显型方法（Explicit method）。

python实现

函数和初始值

欧拉方法解微分方程的关键点在于Δt的选取，Δt越接近0，函数图像越准确
在这里我们将Δt作为预测函数的参数

def fish_predict(Dt):  #Δt
    t_init = 0 #第0天开始
    t_end = 30 #第30天结束
    P_init = 30 #初始数量30
    n_steps = int(round((t_end-t_init)/Dt))

在这里n_steps = (30-0)/Δt ，代表着一共做多少次欧拉方法来绘制函数图像

引入数学分析工具

import numpy as np #矩阵
import matplotlib.pyplot as plt #绘图

建立自变量和因变量矩阵

    t_arr = np.zeros(n_steps + 1) #创建一维矩阵t，记录自变量变化（初始为零） 
    P_arr = np.zeros(n_steps + 1) #创建一维矩阵P，记录因变量变化（初始为零） 
    t_arr[0] = t_init #默认值                      
    P_arr[0] = P_init #默认值

欧拉方法-步进，得到函数离散值

当前节点函数值 = 前一个节点函数值 + 前一个节点的导数 * 自变量变化量
P_arr[i] = P + Dt*dPdt

# Euler's method
    for i in range (1, n_steps + 1):
        P = P_arr[i-1]
        t = t_arr[i-1]
        dPdt = 0.7*P*(1-P/750)-20 #导数                       
        P_arr[i] = P + Dt*dPdt #求当前节点函数值              
        t_arr[i] = t + Dt #自变量步进变化
    return P_arr,t_arr

这样的循环下来，我们就将欧拉方法融入python中，返回两个离散的P_arr、t_arr矩阵，帮助我们描述函数了

在不同变化量下调用函数

为了更加深刻的理解欧拉法求解微分方程，我在这里使用三个不同的变化量使用欧拉方法

p1,t1 = fish_predict(1)
p2,t2 = fish_predict(0.5)
p3,t3 = fish_predict(0.25)

绘图 - 曲线

fig = plt.figure() #创建图像                      
plt.plot(t1, p1, linewidth = 4) #绘制曲线  
plt.plot(t2, p2, linewidth = 4)
plt.plot(t3, p3, linewidth = 4)

绘图 - 标签、轴、网格

plt.title('fish', fontsize = 25) #标签
plt.xlabel('t', fontsize = 20)
plt.ylabel('P(t)', fontsize = 20)
plt.legend(['Dt = 1','Dt = 2','Dt = 3']) #图例
plt.xticks(fontsize = 15) #字体
plt.yticks(fontsize = 15) 
plt.grid(True) #网格                          
plt.axis([0, 25, 0, 800]) #坐标轴范围