ploting second derivative from second order differential equation using for loop

Question

Muhammad Sarmad Zahid 2021-3-20

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/778502-ploting-second-derivative-from-second-order-differential-equation-using-for-loop

评论： Jan 2021-3-20

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

5 个评论
显示 3更早的评论隐藏 3更早的评论

Jan 2021-3-20

在 MATLAB Online 中打开

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

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.
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 2021-3-20

编辑：darova 2021-3-20

在 MATLAB Online 中打开

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 2021-3-20

在 MATLAB Online 中打开

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.

请先登录，再进行评论。

请先登录，再回答此问题。

ploting second derivative from second order differential equation using for loop

5 个评论
显示 3更早的评论隐藏 3更早的评论

回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

ploting second derivative from second order differential equation using for loop

5 个评论 显示 3更早的评论隐藏 3更早的评论

回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

5 个评论
显示 3更早的评论隐藏 3更早的评论