Solving ODE just for one time step

I want to solve some ODE using matlab solver (ode45, ode15s, etc) just for one time step. The time step is decided internally by the solver. Is this possible to do?

[Sam Chak]

Do you mean "one time step" to refer to a single, discrete increment in time during which the solution is approximated by the MATLAB ODE solver from the initial condition?

Could you explain why a "one time step" is required?

[Ari]

Yes, "one time step" refers to a single, discrete increment in time during which the solution is approximated by the MATLAB ODE solver from the initial condition. This "one time step" is required because in every time step I need to change the input function of the ODE. The changes depend on the solution from the solver in every time step. Of course I can solve the ODE for several time steps and take the first solution, but I don't want to waste computational resources.

[Torsten]

Could you explain the problem that makes it necessary to update the ODE after each time step ?

[Steven Lord]

Are you trying to use a fixed-step solver?

Can you show us the mathematical form of the ODE (or perhaps delay-differential equation, DDE) that you're trying to solve? If it is a DDE MATLAB has solvers specifically designed for those types of problems: dde23, ddesd, and ddensd.

[John D'Errico]

The ODE solvers in MATLAB are adaptive things. They don't use a fixed time step. So even if you tell it to give you the solution at a fixed time in the "future", it may well have used multiple intermediate time steps before that point. You would expect exactly that from an adaptive tool.

This means if you ABSOLUTELY need to do as you want, then just choose your choice of ODE schemes. For example, you might use a Runge-Kutta scheme, of some order.

https://en.wikipedia.org/wiki/Runge–Kutta_methods

The code is trivial to write, easier by far than hoping to trick an adaptive tool to not be adaptive! Heck, you can surely find code on the file exchange written by a student where they use Runge-Kutta to solve an ODE.

But, if you look carefully at the idea of a Runge-Kutta algorithm, it forces you to evaluate the derivatives at intermediate points. Is that going against what you are trying to do? I'm not sure, since I have no idea why you are trying to do this.

If you CANNOT think about a solver using intermediate points, because the derivatives you must evaluate depend on the function at that point, then you will be forced to just use a simple Euler scheme, or perhaps a backwards Euler.

[Ari]

To all, this is a reinforcement learning (RL) problem using matlab environment (please don't suggest me to use simulink) by constructing custom step function. In my case the environment is an ODE. For example, the environtment could be a simple water tank model. The step function needs to solve the ODE for just one time step to give the RL agent a chance to receive reward and construct action signal (action signal = new input for the ODE). I am aware that matlab ode solver (such as ode45) has adaptive time step capability. It decides what time steps (the initial time is ) are best for the solution. I don't want to intervene with this capability, in fact I am counting on it. What I want to do is to tell the solver:

"Hey solver, solve this ODE. Use your adaptive time step algorithm to decide the best for the solution, but only calculate and return solution for "

Adaptive time step is very important here for all its merits. If I want to settle with constant time step, I will just use simulink. I hope this clears everything.

[Sam Chak]

Hi @Ari

Thank you for your clarifications. The input function of the ODE (a.k.a. the math model of a dynamical system) typically does not change at every time step, unless you are required to use a piecewise function with If–Else statements to model abrupt changes in dynamics, such as when a falling apple hits the ground. If this is not the case, then it is unnecessary to 'update' the ODE model at every time step.

I believe the confusion arises from the custom Step function, which I have included below. In that function, Euler integration (marked with %%) is used, as @John D'Errico suggested in his comment above. Therefore, you intended to use MATLAB's ODE solvers (such as ode45, ode15s, etc.), but were uncertain about how to implement that in the custom code, am I right?

function [NextObs,Reward,IsDone,NextState] = myStepFunction(Action,State)
    % Custom step function to construct cart-pole environment for the function
    % name case.
    %
    % This function applies the given action to the environment and evaluates
    % the system dynamics for one simulation step.
    
    % Define the environment constants.
    
    % Acceleration due to gravity in m/s^2
    Gravity = 9.8;
    % Mass of the cart
    CartMass = 1.0;
    % Mass of the pole
    PoleMass = 0.1;
    % Half the length of the pole
    HalfPoleLength = 0.5;
    % Max force the input can apply
    MaxForce = 10;
    % Sample time
    Ts = 0.02;
    % Pole angle at which to fail the episode
    AngleThreshold = 12 * pi/180;
    % Cart distance at which to fail the episode
    DisplacementThreshold = 2.4;
    % Reward each time step the cart-pole is balanced
    RewardForNotFalling = 1;
    % Penalty when the cart-pole fails to balance
    PenaltyForFalling = -10;
    
    % Check if the given action is valid.
    if ~ismember(Action,[-MaxForce MaxForce])
        error('Action must be %g for going left and %g for going right.',...
            -MaxForce,MaxForce);
    end
    Force = Action;
    
    % Unpack the state vector from the logged signals.
    XDot = State(2);
    Theta = State(3);
    ThetaDot = State(4);
    
    % Cache to avoid recomputation.
    CosTheta = cos(Theta);
    SinTheta = sin(Theta);
    SystemMass = CartMass + PoleMass;
    temp = (Force + PoleMass*HalfPoleLength*ThetaDot*ThetaDot*SinTheta)/SystemMass;
    
    % Apply motion equations.
    ThetaDotDot = (Gravity*SinTheta - CosTheta*temp) / ...
        (HalfPoleLength*(4.0/3.0 - PoleMass*CosTheta*CosTheta/SystemMass));
    XDotDot  = temp - PoleMass*HalfPoleLength*ThetaDotDot*CosTheta/SystemMass;
    
    %% Perform Euler integration to calculate next state.
    NextState = State + Ts.*[XDot; XDotDot; ThetaDot; ThetaDotDot];
    
    % Copy next state to next observation.
    NextObs = NextState;
    
    % Check terminal condition.
    X = NextObs(1);
    Theta = NextObs(3);
    IsDone = abs(X) > DisplacementThreshold || abs(Theta) > AngleThreshold;
    
    % Calculate reward.
    if ~IsDone
        Reward = RewardForNotFalling;
    else
        Reward = PenaltyForFalling;
    end
end

[Ari]

Yes, @Sam Chak you're right. Here I include an overly simplified and reduced example of the problem. I hope it helps. A better approach is welcomed.

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); % here 1+sin(t) is input function

while sol_t(end) < tmax
    % solve ode for 2 sec time interval.
    % this will produce solution (t0,y0),(t1,y1),...,(tn,yn).
    % for this particular case n is 10 to 12.
    sol = ode15s(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(2)];
    sol_y = [sol_y sol.y(2)];
    
    % we just wasted computational resources in the previous steps.
    % this iteration could be repeated million times.
    % imagine how much computational resources will be wasted.
    % if we could make the solver stop at (t1,y1), or few steps after it, that would be great.
    
    % I tried to use event function to stop the integration,
    % but the solver could only monitor "continuously changing" quantity.
    % should I write my own solver?

    % the RL agent gives new input function here.
    % in this example we just use the old input function.
    
    % for the next time step the initial point will be
    % the previously obtained (t1,y1).
end

% stem(sol_t,sol_y)                 % original by the OP
plot(sol_t, sol_y, 'o'), grid on    % added by Sam

[Sam Chak]

Hi @Ari

The looping appears inefficient because you instructed the ODE solver to integrate over a relatively long time window of 2 seconds (out of a maximum of 6 seconds) when you only want to extract a single solution point right after the initial value and discard the rest at every iteration. Perhaps you can follow @John D'Errico's advice to take the last point. This approach requires you to define or instruct RL to inject or update a new control input signal value to the ODE model at the sampling time.

This approach does not 'waste' resources because the adaptive ODE solver will internally determine the number of intermediate points necessary to accurately perform numerical integration between each sampling time. This allows your RL system to use the sampled data to maximize the cumulative rewards.

sol_t   = 0;        % initial time
sol_y   = 1;        % initial value
window  = 6/60;     % time window (0.1 sec)
tmax    = 6;        % intended simulation time interval
ode     = @(t,y) 0.25*(1 + sin(t)) - 0.1*sqrt(y);     % here 1+sin(t) is input function

while sol_t(end) < tmax
    sol     = ode15s(ode, [sol_t(end) sol_t(end)+window], sol_y(end));
    sol_t   = [sol_t sol.x(end)];
    sol_y   = [sol_y sol.y(end)];
end

figure
disp(sol_t)

  Columns 1 through 18

         0    0.1000    0.2000    0.3000    0.4000    0.5000    0.6000    0.7000    0.8000    0.9000    1.0000    1.1000    1.2000    1.3000    1.4000    1.5000    1.6000    1.7000

  Columns 19 through 36

    1.8000    1.9000    2.0000    2.1000    2.2000    2.3000    2.4000    2.5000    2.6000    2.7000    2.8000    2.9000    3.0000    3.1000    3.2000    3.3000    3.4000    3.5000

  Columns 37 through 54

    3.6000    3.7000    3.8000    3.9000    4.0000    4.1000    4.2000    4.3000    4.4000    4.5000    4.6000    4.7000    4.8000    4.9000    5.0000    5.1000    5.2000    5.3000

  Columns 55 through 62

    5.4000    5.5000    5.6000    5.7000    5.8000    5.9000    6.0000    6.1000

% discard the final solution point because it exceeds tmax
plot(sol_t(1:end-1), sol_y(1:end-1), 'o'), grid on
xlabel('t'), ylabel('y(t)')

[Ari]

Sorry @Sam Chak, I just responded. This is comparison between my approach and your suggestion.

I take (,) as the next point because this is decided by the solver (using its internal adaptive time step mechanism), so the overall solution is in adaptive time step framework. As we can see on the upper plot (I used stem on purpose, instead of plot, to show the differences among time steps) we have adaptive-time-step solution. If we take (), as you and @John D'Errico suggested, then we basically tell the solver to use uniform time step. Also, your solution gives up to (,) for every step, so there isn't much computational saving here.

[Sam Chak]

Hi @Ari

Could you describe the meaning of the desired "computational savings" in @John D'Errico's approach? I do not fully understand it. However, from my perspective, John's "last point" approach does not waste computational effort, as it allows the ODE solver to adaptively determine the internal time step needed for accurate numerical integration between the sampling points.

In contrast, your MPC-like "look-ahead-but-take-only-the-first-computation" approach computes one-third of the maximum time in each iteration. If possible, please provide a desired performance index so that both your approach and John’s can be fairly compared.

[Sam Chak]

For example, despite the sampling time now being set to 1 second, the solution provided by the ode45 solver remains very accurate at the sampling points. The forward Euler method, however, cannot achieve the same level of accuracy with a time step of 1 second.

sol_t   = 0;        % initial time
sol_y   = 1;        % initial value
window  = 1;        % time window (1 sec)
tmax    = 6;        % intended simulation time interval
ode     = @(t,y) 0.25*(1 + sin(t)) - 0.1*sqrt(y);     % here 1+sin(t) is input function

while sol_t(end) < tmax
    sol     = ode15s(ode, [sol_t(end) sol_t(end)+window], sol_y(end));
    sol_t   = [sol_t sol.x(end)];
    sol_y   = [sol_y sol.y(end)];
end

figure
disp(sol_t)

     0     1     2     3     4     5     6

% discard the final solution point because it exceeds tmax
plot(sol_t(1:end), sol_y(1:end), ':o'), grid on
title('Solution provided by the ode45 solver')
xlabel('t'), ylabel('y(t)')

Answer 1

[John D'Errico]

Then I fail to see the problem.

You want only the solution at the end point.

% A simple ODE. with y(0) == 1, the analytical solution is just 2*exp(t)-t-1
odefun = @(t,y) t + y;
[t,y] = ode45(odefun,[0,2],1)

t = 41×1

0 0.0500 0.1000 0.1500 0.2000 0.2500 0.3000 0.3500 0.4000 0.4500

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

y = 41×1

1.0000 1.0525 1.1103 1.1737 1.2428 1.3181 1.3997 1.4881 1.5836 1.6866

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Multiple time steps were generated. So what? TAKE THE LAST ONE.

y = y(end)

y = 11.7781

Is it correct?

2*exp(2) - 2 - 1

ans = 11.7781

What is the problem? This works, basically always. It gave you more information than you wanted, wringing your hands with worry. Take the result you want to see, and ignore the rest.

If it really, really, really, desperately, upsets you that you need to do that extra step, then wrap the ODE solver in a caller function, that does exactly that, returning only the final element of y. Again, what is the problem?

function yfinal = myode45(odefun,tspan,y0)
  [t,y] = ode45(odefun,tspan,y0);
  yfinal = y(end);
end

myode45(odefun,[0,2],1)

ans = 11.7781