fixed time step for ODE45

117 Ansichten (letzte 30 Tage)
aakash d
aakash d am 24 Okt. 2020
Kommentiert: aakash d am 27 Okt. 2020
Hi,
aakash here.
I have written the code below to solve the second order ODE using ODE45. There is an event condition in my code. Because of presence of this event condition, I am not (i dont know) able to use fixed time step in ODE45. Can somebody please help me with this. I request all to please help me to add the fixed time step in this code.
Waiting for the responses.
Thank you
~aakash
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc
clear all
close all
global k Th0 e m c t
%% Input parameters (Length (L), co-efficient of restitution (e), acceleration due to gravity (g))
e=0.9; m=1; c=0; k=246.46;
%% Initial condition
Thid = 8*40/(5*pi)^3; % position from where the mass is released
Thiv = 0; % inital velocity (at t=0)
Th0 = 0.05; % Position where the wall is there
tfi = 5; % time final time
%% ode45 (function call)
eta0=[Thid Thiv]; % Initial conditions
tspan=[0 tfi];
options = odeset('Events',@aakash_function);
t_vec = [];
eta_vec = [];
tcon=0;
ii=1;
while tcon<tfi %% to check whether final time is reached or not
tspan=[tcon tfi];
[t,eta] = ode45(@aakash_mass,tspan,eta0,options);
eta0=[eta(end,1) -(e*eta(end,2))]; %% modified inital conditions after impact
tcon=t(end,1);
t_vec = [t_vec;t];
eta_vec = [eta_vec;eta];
end
t_vec = [t_vec;t];
eta_vec = [eta_vec;eta];
%%
syms u5 x
u5 = sin(5*pi*x)*eta_vec(:,1)
%% Plot
% close all
for i = 1:length(u5)
fig = ezplot(u5(i))
set(fig,'Linewidth',1.5)
xlim([0 1])
ylim([-0.25 0.25])
grid on
title('only third mode')
pause(0.0000000000001)
end
  3 Kommentare
Ameer Hamza
Ameer Hamza am 24 Okt. 2020
ode45() does not use a fixed time-step for solving ODE. It uses an adaptive algorithm to adjust the step-size. Why do you want to use a fixed time-step?
aakash d
aakash d am 24 Okt. 2020
Thanks Ameer for response. I have to check the output at fixed time intervals (say 0.01seconds). Therefore, I want to run code at fixed time steps.
Is there any way to get output at known time intervals. Example: Time duration to run the code is 5 seconds So t_span is from 0 to 5. I want to get response/output at each step of 0.01 second. i.e., 0,. 0.01,0.02,0.03,....5. Please help me if there is any way to get this..

Melden Sie sich an, um zu kommentieren.

Akzeptierte Antwort

Ameer Hamza
Ameer Hamza am 24 Okt. 2020
If you are only concerned about output at fixed time-step, then you can pass tspan as a vector of time-values
tspan= tcon:0.01:tfi;
[t,eta] = ode45(@aakash_mass,tspan,eta0,options);
  1 Kommentar
aakash d
aakash d am 27 Okt. 2020

Thankyou so much Ameer Hamza. That really helped. Your "tspan" related comment worked. Thanks

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (1)

Walter Roberson
Walter Roberson am 24 Okt. 2020
https://www.mathworks.com/matlabcentral/answers/98293-is-there-a-fixed-step-ordinary-differential-equation-ode-solver-in-matlab-8-0-r2012b
has code for fixed step ode such as ode4.
ode45 is inherently variable time step and cannot be used in fixed step.
Event functions are triggered whenever they need to be according to the condition you program in. The ode45 solver will automatically take steps on both sides of the event to find the zero of the event function so as to give as accurate a crossing time as feasible.
If you were to use an event function with a fixed step solver then it would be unlikely that the event would happen to occur exactly at a multiple of the step size. If you were modelling a ball bouncing between two plates, should the fixed step solver fire the event before the ball reaches the plate or after it has passed through it?
  3 Kommentare
Walter Roberson
Walter Roberson am 24 Okt. 2020
No. Your requirements for using the event functions are incompatible with using a fixed time step. You can do what Ameer suggested to pass in a vector of time values, but if the event function fires and is terminal then the last entry will be at the event time rather than at the fixed time.
aakash d
aakash d am 27 Okt. 2020

Okay thanks for response.

Melden Sie sich an, um zu kommentieren.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by