using ode45 using while loop

Hi,

I am trying to solve a dynamics problem using ode45 and while loop. The code is below. It has two case. First a carriage is dropping by free fall and second it is dropping on a damper. The problem is I am getting the value of [at,ay] but when the second case starts the [at,ay] values overwrite the previous at and ay values. Is ther any way to avoid over writing?

clc;
clear;
g=9.81;
tspan = [0 1];
tstart = tspan(1);
tend = tspan(end);
y0=[-0.75 0];
c=6850;
m=450;
s=[];
f=2650;
e=535;
w=4410;
 t = tstart;
 y = y0;
 fcn = @pra;
at=[];
ay=[];
  options = odeset('Events', @freefall);
     
while t(end) < tend 
   
 
  [at,ay] = ode45(fcn, [t(end) tend], y(end,:), options);
   

  t = cat(1, t, at(2:end));
  y = cat(1, y, ay(2:end,:));
   
  if y(end,1)<=0
       
   

      fcn = @pra;
      options = odeset('Events', @freefall);
     
      a = g * ones(1, length(s));
  elseif y(end,1)>0
   
      fcn=@simulation;
      options = odeset('Events', @event);
      a=((c*(ay(:,2))+f+e-w)/m)/9.81;
     
  end       
 end
  
        


figure(1);                           
plot(t,y(:,1))                            
hold on
figure(2);
plot(t,y(:,2))
hold on
figure
plot (at,a)
%%%%% functions%%%%%%
function fval = simulation( t,y )

x=y(1);
v=y(2);

   c=6850;
        m=450;
        k=0;
        f=2650;
        e=535;
        w=4410;

fval(1,1)=v;
 %fval(2,1)= -(c*v+k*x)/m
 fval(2,1)=(-c*v-f-e+w+k*x)/m;


end
function val = pra( t,y )
g=9.8;
y(1);
y(2);


% c=6850;



val = [y(2); g];

end
%%event function%%%%
function [position,isterminal,direction] = freefall(t,y)
position = y(1); % The value that we want to be zero
isterminal = 1;  % Halt integration 
direction = 0;   % The zero can be approached from either direction

Answer 1

[Jan]

I am getting the value of [at,ay] but when the second case starts the [at,ay] values overwrite the previous at and ay values.

Why does this matter? You insert the values of at and ay to the arrays and y. So why do you want to access at and ay later? Use t and y instead.

[Muhammad Sarmad Zahid]

JI jan,

I actually need the previous values as well because i have to make a plot through those time values. is there any way to avoid overwriting ? I know i am saving all the time values but for the first case ineed seperate array of ay and at

[Jan]

But you do store all values of at inside t and all parts of ay in y already:

  [at,ay] = ode45(fcn, [t(end) tend], y(end,:), options);
   
  t = cat(1, t, at(2:end));
  y = cat(1, y, ay(2:end,:));  

After the loop all the data are stored together. Therefore I do not understand, why you ask for storing them another time.

But if you have a good reason to store the data twice:

tC = {};   % Before the loop
yC = {};
while t(end) < tend
  [at, ay] = ode45(fcn, [t(end) tend], y(end,:), options);
  
  % Collect the output
  t = cat(1, t, at(2:end));
  y = cat(1, y, ay(2:end,:));
  
  % Store the output another time:
  tC{end + 1} = at;
  yC{end + 1} = ay;
  
  ...
end

[Muhammad Sarmad Zahid]

Jan, I think you are storing all the value again. But i want to store the values of at and ay for each event seperately

[Jan]

But i want to store the values of at and ay for each event seperately

And this is exactly, what my code does: For the i.th interval the times tC{i} and the trajectory yC{i} is stored. If you need any other kind of "separately" post an example of what you mean.

[Muhammad Sarmad Zahid]

Jan, the problem now is resolved that i can save the ay at values of each event seperately. What i was trying to do is to calculate the second derivative. I still cant calculate the second derivative for first event. I want to combine the second derivatives of both first and second event in a plot. is it possible to do that ? just like what i am doing for displacement and velocity?

[Muhammad Sarmad Zahid]

The reason why I wanted to get t and y seperately beacause i wanted to draw a graph of second derivative. but it does not seems to work. If I ry to explain my goal in one sentence that is " i want to get acceleration plot over the whole time that is second derivative" . Your guidance will be valueable if you guide me how can I get seconf derivation for each event.

[Jan]

dy  = [];   % insert these lines to your code
tC  = {};
yC  = {};
dyC = {};
while t(end) < tend
  [at, ay] = ode45(fcn, [t(end) tend], y(end,:), options);
  
  dy  = simulation(tsol, ysol.').';
  ady = dy(:, 2);
  
  % Collect the output
  t  = cat(1, t, at(2:end));
  y  = cat(1, y, ay(2:end,:));
  dy = cat(1, dy, ady(2:end,:));
  
  % Store the output another time:
  tC{end + 1}  = at;
  yC{end + 1}  = ay;
  dyC{end + 1} = ady;
  ...
end


function fval = simulation(t, y)
x = y(1, :);
v = y(2, :);
c = 6850;
m = 450;
k = 0;
f = 2650;
e = 535;
w = 4410;
fval(1, :) = v;
fval(2, :) = (-c * v - f - e + w + k * x) / m;
end

Answer 2

[Muhammad Sarmad Zahid]

The dimensions of dy and ady were not same but i corrected it. the answer i getting is as follow. but dont you think that fot the first event it should be simply 9.8 m/s or 1 g? from 0s ro 0.3912s it should be constant but here we are getting a slope. I actually did it my self as well but i could not solve this problem. you can check the attached file for acceleration graph.

[Jan]

the answer i getting is as follow.

Your diagram does not contain a title or labels. Without seeing the code I can only guess, what was plotted. I do iknow know the final code you are using. So I cannot find out, if it contains a problem.