pythonnumerical-methodsdifferential-equations

One adaptive-time-step solution for ODE


I want to solve some ODE just for one time step. The time step is decided internally by the solver, not by me. Is there a python adaptive-time-step ode solver which could do this? For illustration the solver may look like this

solution = solve(ode_func, t_span, init_cond, step_num) 

If step_num = n then the solver calculates and returns solution for t1,t2,...,tn only (initial time is t0). Those t1,t2,...,tn are decided by the solver using its adaptive-time-step algorithm, I simply supply the value of n. In my case, I will set step_num = 1. If anyone is wondering, this one-time-step solver is required because I need to update the ode_func every time step.

Of course I could use any adaptive-time-step ode solver and then just take solution for t1 but I don't want to waste computational resources here.

Update: This is my attempt in matlab. It's very inefficient since matlab solvers don't have one-step-integration capability.

sol_t = 0; % initial time
sol_y = 1; % initial value
window = 2; % time window (2 sec)
tmax = 6;
ode = @(t,y) 0.25*(1+sin(t))-0.1*sqrt(y);
step = 1;

while sol_t(end) < tmax
    % solve ode for 2 sec time interval.
    % this will produce solution (t0,y0),(t1,y1),...,(tn,yn).
    sol = ode45(ode, [sol_t(end) sol_t(end)+window], sol_y(end));

    % take just (t1,y1) and discharge the rest
    sol_t = [sol_t sol.x(step+1)];
    sol_y = [sol_y sol.y(step+1)];
    
    % do something here
end

stem(sol_t,sol_y)

Solution

  • you can use scipy.integrate solvers directly, for example RK45

    from scipy.integrate import RK45
    
    def func(t, y):
        return -0.5 * y
    
    t = []
    y = []
    t_final = 10
    t.append(0)
    y.append([1])
    solver = RK45(func, t0=t[0], y0=y[0], t_bound=t_final)
    while solver.status == 'running':
        error = solver.step()
        if error:
            print(error)
            break
        t.append(solver.t)
        y.append(solver.y)
    
    # plot result with matplotlib
    import matplotlib.pyplot as plt
    
    plt.plot(t, y)
    plt.show()
    

    The curve looks jagged because the points are joined with a line, you can interpolate the result at each step with solver.dense_output()(some_time)


    I need to update the ode_func every time step.

    you shouldn't. the solvers assume the function does not change, solvers are very very sensitive to discontinuities and discontinuities in any derivatives. if you are okay with discontinuous derivatives then you can create a new solver on each step.

    solver = RK45(func, t0=t[-1], y0=y[-1], t_bound=t_final)
    error = solver.step()
    if error:
        raise ValueError(error)
    t.append(solver.t)
    y.append(solver.y)
    
    # create a new solver starting from the end time of the last one, and a different function
    solver = RK45(func2, t0=t[-1], y0=y[-1], t_bound=t_final)
    # use the new solver
    

    most solvers run code to set their internal states in the initializer, if the problem changes then this code needs to re-run by recreating the solver, for example RK45 needs to re-estimate the initial step size if it is not given.

    keep in mind that you will have discontinuous derivatives at each step, this only prevents the solver state from carrying over from the last step, so if whatever system you are solving should have continuous derivatives then you have the wrong solution here, or rather you are trying to solve the wrong problem.

    Full code:

    from scipy.integrate import RK45
    
    def func1(t, y):
        return -0.5 * y
    def func2(t, y):
        return -0.2 * y
    
    t = []
    y = []
    t_final = 10
    t.append(0)
    y.append([1])
    
    current_function = func1
    solver = RK45(current_function, t0=t[-1], y0=y[-1], t_bound=t_final)
    while solver.status == 'running':
        if current_function is func1:
            current_function = func2
        else:
            current_function = func1
        solver = RK45(current_function, t0=t[-1], y0=y[-1], t_bound=t_final)
        error = solver.step()
        if error:
            print(error)
            break
        t.append(solver.t)
        y.append(solver.y)
    
    # plot result with matplotlib
    import matplotlib.pyplot as plt
    
    plt.plot(t, y)
    plt.show()
    

    You can see how discontinous the output derivative looks. enter image description here

    Realistically speaking, i only needed to do this to implement a switch, you don't need to do this to solve ODEs .... at least you now know how to implement a switch, that's how you create accurate first order discontinuities.