ploting second derivative from second order differential equation using for loop

Hi I am trying to solve a dynamics problem. I have tried to simplify my code and have successfully ploted zero and first derivatives of two cases. The code is below. How can i call the function so that I can get second derivative of both events. Right now i am only getting second derivativeof first case

here is my code

clc;
clear all;
g=9.81;

y0=[-0.75 0];

c=6850;
m=450;
k=0;
f=2650;
e=535;
w=4410;
tspan = [0 1];
tstart = tspan(1);
t = tstart;
y = y0;
at=[];
ay=[]; 


for i=1:100    
    tstart=(i-1)*0.01;
    tend=i*0.01;    
    [at,ay] = ode45(@(t,y) simulation(t,y), [tstart tend], y(end,:));     
    t = cat(1, t, at(2:end));
    y  = cat(1, y, ay(2:end,:));    
    dy = simulation(t, y.').';       
end

plot (t,y(:,1))
figure
plot (t,y(:,2))
figure
plot(t,dy(:,2))
     %%%%%%function%%%%%%%%
function fval = simulation( t,y )
    x = y(1, :);
    v = y(2, :);    
    c=6850;
    m=450;
    k=0;
    f=2650;
    e=535;
    g=9.8;
    if y(1, :)>=0           
          F = -c * v - f - e ;          
    else    
           F=0;     
    end
    fval(1, :) = v;
    fval(2, :) = (F + (m*g) + k * x) / m;
end

[Jan]

I've asked you in your other questions to use the tools to format the code in the forum. Please, Muhammad, do not ignore my advice, because this makes your code more readable. Use a proper indentation also: Simply press Ctrl-a Ctrl-i in the editor. Removing the number of empty lines would be useful also.

I've showed you how to obtain the 2nd derivative twice already in your other questions. You have modified my suggestions and get the value already:

dy = simulation(t, y.').';

Doing this inside the loop is less useful, because this evaluates the 2nd derivative for all t and y repeatedly. Better call this after the loop.

I do not see, how this question differes from your previous questions.

[Muhammad Sarmad Zahid]

Jan thank you for letting me know the commands. I will keep that in mind next time. I apologize.

I did both things (putting inside and outside loop) but i am getiing a constant graph at value 9.81 through all time from 0 to 1s. although it should be constant from 0 to 0.3912 s and then the acceleration graph should change after the change of condition when F is not equal to zero.

[Jan]

There is no reasons to apologize. If your code is hard to read, and this will impede the formulation of answers, you are the only person how suffers. But what about editing your question and making the code more readable now? This takes some seconds only.

You have inserted this in the function to be integrated:

if y(1, :) >= 0
   F = -c * v - f - e ;
else
   F = 0;
end

1. Now the function is not differentiable anymore. Matlab's ODE45 integrator is designed to handle differentiable functions only. Otherwise the step size control can be confused such massively, that the resulting trajectory can be dominated by rounding errors. So don't do this, but let en event function change the function to be integrated.
2. Remember that Matlab's IF command requires a scalar condition. You provide a vector here. Then Matlab inserts an all(...) automatically. During ODE45 runs, this funtion is called with a scalar only, but for obtaining the 2nd derivative, y(1, :) is a vector and F is set to0 if any element of y is not >= 0.

[Muhammad Sarmad Zahid]

Jan, but the problem is same while using event functions. as well there i no step change in the velocity. the code is below. This is my previous code using while loop. My velocity and displacement curves are ideal but i am trying to calling a function.

clc;
clear;
dy  = [];
tC  = {};
yC  = {};
dyC = {};
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);

tC = {};   % Before the loop
yC = {};
% figure
% plot(tsol,ysol(:,1), 'r-');
% hold on;
% plot(t,y(:,1), 'r-');
% figure
% plot(tsol,ysol(:,2), 'g-');
% hold on;
% plot(t,y(:,2), 'g-');

while t(end) < tend        
    [at,ay] = ode45(fcn, [t(end) tend], y(end,:), options);
    Q1 = [t(end), tend; at(1), at(end); ay(1,:); ay(end,:)]
    % Collect the output
    t  = cat(1, t, at(2:end));
    y  = cat(1, y, ay(2:end,:));
    
    dy  = fcn(t.', y.').';
    d2y = dy(:, 2);
    
    % Store the output another time:
    tC{end + 1}  = at;
    yC{end + 1}  = ay;       
    if y(end,1)<=0                        
        fcn = @pra;
        options = odeset('Events', @freefall);                
    elseif y(end,1)>0        
        fcn=@simulation;
        options = odeset('Events', @event);                        
    end
end
figure(1);
plot(t,y(:,1))
hold on
figure;
plot(t,y(:,2))
hold on
figure
plot(t,d2y)

[Jan]

These lines of your code are meaningless:

dy  = [];
g=9.81;
c=6850;
m=450;
s=[];
f=2650;
e=535;
w=4410;
at=[];
ay=[];

   Q1 = [t(end), tend; at(1), at(end); ay(1,:); ay(end,:)]   

Filling your program with useless line will increase the confusion level. So I recomment to omit them.

I still have no idea, why you collect the output twice:

    % 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;

You cannot collect the value of d2y in exactly the same way, because you have used the method to start t and y with their initial value, but d2y has no initial value at startup. So you need this:

d2y = [];
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,:));
    
    dy  = fcn(t.', y.').';
    if t(end) < tend  % Crop last element:
       d2y = cat(1, d2y, dy(1:end-1, 2);
    else
       d2y = cat(1, d2y, dy(:, 2);
    end
    ...
        

Then plot(t,d2y) will work also.