How do I include a whileloop in the Ode45?

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.

[Star Strider]

I suggest you explore the ‘Events’ property of odeset. (See the section on ‘ODE Events Property’ about half-way down the page.) It may do what you want.

[Jan]

@Star Strider: Please post answers in the answer section. Otherwise they cannot be accepted.

[Star Strider]

@Jan — Noted. Thank you.

This was more a suggestion than anything else. I eventually intend to experiment with it and ‘Events’ and post that result as an answer.

In the context of ODEs, thank you for your discussions on numerics wrt the ode solvers. They sent me back to the books (literally), to my benefit.

[Star Strider]

@Hassan — What is your atmosisa function? I know it must be the Interantional Standard Atmosphere at y(2), but I cannot simulate your function without being able to call it.

Also, please put an end statement at the end of FS1, and please put some comments indicating what your code does as it runs. It is very difficult for me to interpret your code.

[Hassan]

@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]

@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.

Answer 1

[Azzi Abdelmalek]

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)

[Hassan]

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]

maby I did'nt understand your question try the corrected code above