using ode45 using while loop

1 次查看（过去 30 天）
Muhammad Sarmad Zahid
评论： Jan 2021-3-19
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

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

采纳的回答

Jan
Jan 2021-3-17
编辑：Jan 2021-3-17
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.
  7 个评论
显示 5更早的评论
Muhammad Sarmad Zahid
Muhammad Sarmad Zahid 2021-3-19
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
Jan 2021-3-19
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

请先登录，再进行评论。

更多回答（1 个）

Muhammad Sarmad Zahid
编辑：Muhammad Sarmad Zahid 2021-3-19
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.
  1 个评论
Jan
Jan 2021-3-19
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.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by