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

Question

我正在尝试解决这个写入文件的问题，时间运行s 从 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.

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

Answer 1

你只有最后一个值，因为你在循环执行期间没有输出任何东西，只是在它结束时。

为了解决你的问题，我会这样做，留给 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'))

Re "I am met with this error"，你永远不会更新（至少在你展示的代码中）初始条件的值，所以很难猜出是什么出错了，但我想你使用了 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)

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

Writing data from equation solver to file over loop only outputs last 2 elements

python

python-3.x

scipy

differential-equations

odeint