Search code examples
pythondifferential-equations

Develop a ode integrator stepper


I'm creating a class in order to perform several operations to a differential problems (like solve, plot, save into file the solution) The stepper is the more important part and I want to modified my code from simple classmethod to a generator

Here a snap of my code usable:

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()

It is possible to create a generator that for the 3 calls at the stepper give back the value compute inside the if,elif,else block (here I wrote a simple method but in the stepper I make the call commented) if yoes Is a good idea and how can I do ?

EDIT jdowner I have comment out the calling to runge-kutta because I would just present a minimal working code the calls to runge-kutta are done for start-up the multistep (4step so` I need to calculate 3 point in order to start the method)

EDIT Ilia L but I need this 3 variable also during the run of the multistep method .. look inside the for

EDIT i got this error :

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

Solution

  • If you want to perform the above computation using a generator, I would write it like,

    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