时间戳

首页 > 解决方案 > 通过循环将数据从方程求解器写入文件仅输出最后 2 个元素

python - 通过循环将数据从方程求解器写入文件仅输出最后 2 个元素

问题描述

我正在尝试解决写入文件的问题，时间从 0 -> 1000 以 100 为增量运行，但是当我运行脚本时，我一直只得到输出的最后 2 个元素。我认为循环存在问题，但我不确定。我希望所有数据都存在于 t=0 到 t=1000 之间，而不仅仅是最后的 2 个元素。可以在文件中看到，只记录了t=900和t=1000。似乎看不到要改变什么。

import matplotlib.pyplot as plt
from scipy.integrate import odeint
import pandas as pd


plt.ion()
plt.rcParams['figure.figsize'] = 10, 8

P = 0      # birth rate
d = 0.0001  # natural death percent (per day)
B = 0.0095  # transmission percent  (per day)
G = 0.0001  # resurect percent (per day)
A = 0.0001  # destroy percent  (per day)

# solve the system dy/dt = f(y, t)
def f(y, t):
     Si = y[0]
     Zi = y[1]
     Ri = y[2]
     # the model equations (see Munz et al. 2009)
     f0 = P - B*Si*Zi - d*Si
     f1 = B*Si*Zi + G*Ri - A*Si*Zi
     f2 = d*Si + A*Si*Zi - G*Ri
     return [f0, f1, f2]

# initial conditions
S0 = 500.             # initial population
Z0 = 30                 # initial zombie population
R0 = 60                 # initial death population
y0 = [S0, Z0, R0]     # initial condition vector

#looping over some time instead of integrating in one go.
t_a = 0 
oput = 500
t_b = t_a + oput
delta_t = t_a + 100 
tend = 1000

dfs=[]
while t_a < tend: 

    t_c = t_a + delta_t 
    t=[t_a,t_c]
    y = odeint(f,y0,t,mxstep=10000) #Integrator
    t_a = t_c

    if(t_a > oput):
        t_b = t_b +oput

        S = y[:,0]
        R = y[:,1]
        Z = y[:,2]



dfs.append(pd.DataFrame({'t': t, 'Z': Z,'R': R}))
g=pd.concat(dfs,axis=0)
g.to_csv('example.csv',mode='w',index=False)

出现的另一个问题是我想将我的初始条件更改为求解器的最终元素：

        S0 = S
        R0 = R
        Z0 = Z
        y0 = [S0,R0,Z0]

但我遇到了这个错误：

ValueError: Initial condition y0 must be one-dimensional.

如果我想用循环中的弹出值更新初始条件，我该如何解决这个错误？

标签： pythonpython-3.xscipydifferential-equationsodeint

解决方案

您只有最后一个值，因为您在循环执行期间不输出任何内容，而只是在循环结束时输出。

为了解决你的问题，我会这样，留给odeint循环时间步骤的任务，

In [21]: from scipy.integrate import odeint
In [22]: P, D, B, G, A = 0, 0.0001, 0.0095, 0.0001, 0.0001
In [23]: def fy(y, t):
    ...:      Si, Zi, Ri = y
    ...:      # the model equations (see Munz et al. 2009)
    ...:      return P-B*Si*Zi-D*Si, B*Si*Zi+G*Ri-A*Si*Zi, D*Si+A*Si*Zi-G*Ri
In [24]: y0 = [500.0, 30.0, 60.0]
In [25]: res = odeint(fy, y0, [i*100.0 for i in range(11)], mxstep=1000)
In [26]: print('\n'.join(','.join('%+15.9f'%x for x in row) for row in res))
 +500.000000000,  +30.000000000,  +60.000000000
   +0.000000000, +525.356080701,  +64.643919299
   -0.000000000, +525.999298342,  +64.000701658
   +0.000000000, +526.636115893,  +63.363884107
   +0.000000000, +527.266597017,  +62.733402983
   +0.000000000, +527.890804714,  +62.109195286
   +0.000000000, +528.508801154,  +61.491198846
   +0.000000000, +529.120648555,  +60.879351445
   -0.000000000, +529.726408164,  +60.273591837
   -0.000000000, +530.326140479,  +59.673859521
   -0.000000000, +530.919905407,  +59.080094593

如果要将 CSV 格式保存到文件，print支持重定向到打开的文件，例如print(..., file=open('res.txt'))

重新“我遇到了这个错误”，你永远不会更新（至少在你显示的代码中）初始条件的值，所以很难猜测出了什么问题，但我想你使用了类似的y0 = y，虽然它应该是y0 = y[1]——如果你打印y你就会明白为什么。

最后，据我所知，您的时间步长可能比期望的要长一些，因为第一个时间步长结束时的活人口为 0（更准确地说，2.72457580e-13）并且您可能想要遵循过程而不是看它的最终结果。

附录

作为您评论的后续行动，您可以遍历结果、检查它们并采取适当的措施，例如

res = odeint(...)

for S, R, Z in res:
    output_style = 0
    if S<40 and R> 30 : output_style = 1 
    elif ............ : output_style = 2
    ....
    output(S, R, Z, output_style)

推荐阅读