Solving ODE just for one time step
Show older comments
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?
12 Comments
Torsten
on 6 May 2025
Could you explain the problem that makes it necessary to update the ODE after each time step ?
Steven Lord
on 6 May 2025
John D'Errico
on 6 May 2025
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.
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.
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:
(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 
"
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
on 7 May 2025
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
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
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)
% 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
on 16 May 2025
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.
,) 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
on 16 May 2025
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.
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)
% 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)')
Answers (1)
John D'Errico
on 7 May 2025
Edited: John D'Errico
on 7 May 2025
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)
Multiple time steps were generated. So what? TAKE THE LAST ONE.
y = y(end)
Is it correct?
2*exp(2) - 2 - 1
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)
Categories
Find more on General Applications in Help Center and File Exchange
on 16 May 2025
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
