Solving ODE just for one time step

Question

Ari 2025-5-6

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2176929-solving-ode-just-for-one-time-step

评论： Sam Chak 2025-5-16

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 个评论
显示 10更早的评论隐藏 10更早的评论

John D'Errico 2025-5-6

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 2025-5-7

编辑：Ari 2025-5-7

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 2025-5-7

在 MATLAB Online 中打开

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 2025-5-8

编辑：Sam Chak 2025-5-8

在 MATLAB Online 中打开

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 2025-5-8

编辑：Sam Chak 2025-5-8

在 MATLAB Online 中打开

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

Sam Chak 2025-5-16

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 2025-5-16

在 MATLAB Online 中打开

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 2025-5-7

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2176929-solving-ode-just-for-one-time-step#answer_1564834

编辑：John D'Errico 2025-5-7

在 MATLAB Online 中打开

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Solving ODE just for one time step

12 个评论
显示 10更早的评论隐藏 10更早的评论

回答（1 个）

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

Solving ODE just for one time step

12 个评论 显示 10更早的评论隐藏 10更早的评论

回答（1 个）

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

12 个评论
显示 10更早的评论隐藏 10更早的评论

0 个评论
显示 -2更早的评论隐藏 -2更早的评论