Home Articles

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

Asked
Viewed 1729 times
0

我正在尝试解决写入文件的问题,时间从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值更新初始条件,我该如何解决这个错误?

1 Answer

  • 0

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

    为了解决你的问题,我会这样做,留下循环时间步骤的任务 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'))

    再"I am met with this error",你永远不更新(至少在你的代码've shown) the values of the initial conditions, so that it'很难猜测出了什么问题,但我想你使用类似的 y0 = y ,而应该是 y0 = y[1] - 如果打印 y 你就会明白为什么 .

    最后,从我能理解,也许你的时间步骤有一点比什么是理想的时间越长,如在第一时间步骤结束的常住人口为0(更准确地说, 2.72457580e-13 ),并有可能要遵守过程而不是看最终结果 .


    Addendum

    作为your comment的后续行动,您可以循环搜索结果,检查结果并采取适当的措施,例如:

    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)
    

Related