清空筛选条件
清空筛选条件

How do I include a whileloop in the Ode45?

1 次查看（过去 30 天）
Hassan
Hassan 2012-9-22
Hey!
I´m trying to solve a set of equation that describes the motion of a body in a ballistic trajectory. I'm interested of knowing at what downrange the body will hit the ground so I want to set a limit for the altitidue that H>0.
How can I include a whileloop in the ode45 function to set that limit.
Here is the code:
function [yp] = FS1(t,y)
ms1 = 698390; %structure mass for the first stage
[T, a, P, rho] = atmosisa(y(2));
g0=9.8; %Gravity at sea level [N/m^2]
Re=6371000; %Earth radius [m]
g=g0*(Re/(Re+y(2)))^2; %Gravitational variation
A=10; % Reference area
Cd=0; % Drag coefficient
D=0.5*rho*y(3)^2*Cd*A;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%The four first order differential equations yp=zeros(4,1);
yp(1)=y(3)*cos(y(4));
yp(2)=y(3)*sin(y(4));
yp(3)= -D/ms1 -g*sin(y(4)) +(((y(3)*(cos(y(4))))^2)*sin(y(4)))/(Re+y(2));
yp(4)= (-cos(y(4))*g)/y(3) +((y(3)*cos(y(4)))^2*cos(y(4)))/(y(3)*(Re+y(2)));
%%%%%%%%%%%%%%%%%%%%%%%%%%
y0 = zeros(4,1);
y0(1,1) = XR1(end); %Initial downrange, can be set to zero
y0(2,1) = HR1(end); %initial altitude, can be set to 7.6730e+004
y0(3,1) = VR1(end); %initial velocity, can be set to 2.7684e+003
y0(4,1) = gammaR1(end); %initial angel of attack, can be set to 0.5599
tspan = [0 500];
[t,y] = ode45(@FS1,tspan,y0);
XF1 = y(:,1);
HF1 = y(:,2);
VF1 = y(:,3);
gammaF1 = y(:,4);
figure(3)
plot(XF1,HF1)
The initial values is not that important, you can put the values you like. I will be very thankfull for any help.
  6 个评论
显示 4更早的评论
Hassan
Hassan 2012-9-23
@Star Strider - y(2) is the value of the altitude and it's changing value for ecry step of the solving process, for every new y(2) that is calculated the atmosisa function is calling that value. It should be called by itself in the FS1 function.
I'm sorry for the lack of comments, I have tried to put more explaning comments and I hope that should clarify a little.
Thank you
Star Strider
Star Strider 2012-9-23
@Hassan — Thank you for your clarification. The problem is that without the atmosisa function, I cannot run your code and try my ‘Events’ option with it.
In any event, if you have solved your problem (as you seem to have in your comment to Azzi's answer), then any solution I can come up with is likely unnecessary anyway.

请先登录，再进行评论。

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

回答（1 个）

Azzi Abdelmalek
Azzi Abdelmalek 2012-9-23
编辑：Azzi Abdelmalek 2012-9-23
correted code
y0 = zeros(4,1);
y0(1,1) = XR1(end);
y0(2,1) = HR1(end);
y0(3,1) = VR1(end);
y0(4,1) = gammaR1(end);
tf=500
test=1
while test>0
tf=tf+100
tspan = [0 tf];
[t,y] = ode45(@FS1,tspan,y0);
HF1 = y(:,2);
idx=find(HF1<0,1)
if ~isempty(idx)
test=0
end
end
XF1 = y(1:idx,1);
HF1 = y(1:idx,2);
VF1 = y(1:idx,3);
gammaF1 = y(1:idx,4);
figure(3)
plot(XF1,HF1)
  2 个评论
Hassan
Hassan 2012-9-23
编辑：Hassan 2012-9-23
Thank you for your answer Azzi but it's not working for me, I tried this solution but when I ran the code idx was equal to 1, that gave me XF1 and HF1 as scalar numbers and not vectors as I want them. Anyway I've solved the problem by using this if statement:
HF1 = [];
for i=1:length(y(:,2))
HF1 = [HF1 y(i,2)];
if y(i,2)<=0
break
end
end
if HF1(end)<0
HF1 = HF1(1:end-1);
elseif HF1(end)>0
HF1 = HF1;
end
I would be glad if you could comment on this because I would more like to use the whileloop to solve this problem.
Sincerely
Azzi Abdelmalek
Azzi Abdelmalek 2012-9-23
maby I did'nt understand your question try the corrected code above

请先登录，再进行评论。

类别

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

标签

Community Treasure Hunt

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

Start Hunting!

Translated by