代码之家  ›  专栏  ›  技术社区  ›  Drudox lebowsky

一种ode积分器步进器的研制

  •  0
  • Drudox lebowsky  · 技术社区  · 7 年前

    我正在创建一个类,以便对一个差异问题执行多个操作(例如solve、plot、save in to file the solution) 步进器是更重要的部分,我想把我的代码从简单的classmethod修改为生成器

    下面是我的代码片段:

    import numpy as np  
    
    first_startup , second_startup, third_startup = True, True, True
    
    
    um1,um2,um3 = 0.,0.,0.
    
    def step(f , t : np.float , u : np.float , dt):
        global first_startup , second_startup, third_startup
        global um1,um2,um3    
    
        if first_startup:
            um3 = u.copy()
            unext = u + dt*f(t,u)  #rungekutta.RK4.step(func,t,um2,dt)
            t += dt
            first_startup = False
        elif second_startup:
            um2 = u.copy()
            unext = unext = u + dt*f(t,u) #rungekutta.RK4.step(func,t,um2,dt) 
            t+= dt
            second_startup = False
        elif third_startup:  
            um1 = u.copy()
            unext = u + dt*f(t,u) #rungekutta.RK4.step(func,t,um1,dt)
            t += dt 
            third_step = False
        else:  # compute AB 4th order    
            unext = u + dt/24.* ( 55.*f(t,u) - 59.*f(t-dt,um1) + 37.*f(t-dt-dt,um2) \
                                                                   - 9.*f(t-dt-dt,um3 )) 
            um3 = um2.copy()
            um2 = um1.copy()
            um1 = u.copy() 
        return unext   
    
    def main():
    
        func = lambda t,u : -10*(t-1)*u
        t0 = 0. 
        tf = 2.      
        dt = 2/50 
        u = np.exp(-5)   
        t = t0
        with open('output.dat', 'w') as f:
            while True:
                f.write('%f %f \n' %(t, u) )
                u = step(func,t,u,dt)
                t += dt
                if t > tf:
                  break
    
    if __name__ == '__main__':
        main()
    

    可以为步进器上的3个调用创建一个生成器,该生成器将if、elif、else块中的值计算返回(这里我编写了一个简单的方法,但是在步进器中我对调用进行了注释)如果是一个好主意,我该怎么做?

    编辑 jdowner我已经注释掉了对runge-kutta的调用,因为我只需要提供一个最小的工作代码,对runge-kutta的调用是为启动多步骤而完成的(4步,因此`我需要计算3点才能启动该方法)

    编辑 但我也需要这3个变量,在运行多步方法的过程中在里面寻找

    编辑 我得到这个错误:

    Traceback (most recent call last):
      File "drive.py", line 377, in <module>
        main()
      File "drive.py", line 224, in main
        f.write('%f %f \n' %(t, u) ) # u[0]) )
    TypeError: must be real number, not generator
    
    1 回复  |  直到 7 年前
        1
  •  1
  •   jdowner    7 年前

    如果你想用生成器执行上面的计算,我会这样写,

    def step(f, t, u, dt):
    
        um3 = u
        yield u + dt * f(t, u)
    
        um2 = u
        yield u + dt * f(t, u)
    
        um1 = u
        yield u + dt * f(t, u)
    
        while True:
            k1 = 55.0 * f(t, u)
            k2 = 59.0 * f(t - dt, um1)
            k3 = 37.0 * f(t - 2 * dt, um2)
            k4 = 9.0 * f(t - 2 * dt, um3)
    
            yield u + dt * (k1 - k2 + k3 - k4) / 24.0
    
            um3 = um2
            um2 = um1
            um1 = u