Using MATLAB to do Numeric Calculus

So I am culminating all of my newfound knowledge to solve and show a simple calculus problem. However I am running into issues...

Given: A piecewise function where v(t)=70t @ t<=20s & 1596.2-9.81t @ 20<t<t_max <---(eq 1)

Also given: s(t)=35t^2+7620 @ t<=20s & 1596.2t-4.905t^2-8342 @ 20<t<=t_max <---(eq 2)

Find: If aircraft ignites its engines at t=0s and accelerates vertically, there is only enough fuel for an engine burn of 20 seconds. At which time will the ship become a projectile? Find value of t_max. Plot the crafts velocity vs. t, (eq 1). Plot the crafts altitude vs. t, (eq 2).

My Solution:

%% Finding Value of tmax
% Finding tmax and rounding it to the nearest whole number
%%

tmaxr=roots([-9.81, 1596.2]);
tmax=ceil(tmaxr);

%% Plotting Equation 1
% Creating Figure 1 from Equation 1
%%

t=0:0.1:tmax;
v=zeros(size(t));
for N=1:length(t)
   if t(N)<=20
       v(N)=70*t(N);
   else
       v(N)=1596.2-9.81*t(N); % Only valid for t > 20 seconds
   end 
end
figure(1)
plot(t,v)
xlabel('Time (s)')
ylabel('Velocity (m/s)')
title('Velocity vs. Time')

%% Plotting Equation 2
% Creating Figure 2 from Equation 2
%%

t=0:0.1:tmax;
s=zeros(size(t));
for N=1:length(t)
   if t(N)<=20
       s(N)=35.*t(N).^2+7620;
   else
       s(N)=1596.2.*t(N)-4.905.*t(N).^2-8342;
   end 
end
figure(2)
plot(t,s)
xlabel('Time (s)')
ylabel('Altitude (m)')
title('Position vs. Time')

Issue: So I have found t_max, plotted the velocity but it looks funky, and also tried to plot my altitide but not sure if it's faultless.

P.s. Sorry for the formatting, there's no easy way to type piecewise functions that I know of.

Answer 1

[Torsten]

Your code is correct.

Here is a simpler way to define s and v. Just substitute tcrash by tmax.

tcrash = 1596.2/(2*4.905)+sqrt((1596.2/(2*4.905))^2-8342/4.905);
s = @(t)(35*t.^2+7620).*(t<=20)+(1596.2*t-4.905*t.^2-8342).*(t>20).*(t<=tcrash);
v = @(t)(70*t).*(t<=20)+ (1596.2-9.81*t).*(t>20).*(t<=tcrash);
t = 0:0.01:tcrash;
figure(1)
plot(t,s(t))
grid on

figure(2)
plot(t,v(t))
grid on

[Spaceman]

I think I am missing something in my position code as it does not look like a parabola as it should... Any thoughts? Yours came out good, but I am trying to use a for-loop.

[Torsten]

Your curves are the same as mine - only plotted up to the point of maximum altitude.

[Spaceman]

I see. I have also found the derivative approximation of velocity and plotted it to compare it to the equation:

ds_dt=diff(s)./diff(t); % Meters per second
figure(3)
plot(t,v)
hold on
plot(t(1:end-1),ds_dt,'o','MarkerIndices',1:10:length(ds_dt))
xlabel('Time (s)')
ylabel('Velocity (m/s)')
title('Velocity vs. Time')

I think everything looks right...

[Torsten]

syms t
s1 = 35*t^2+7620 ;
s2 = 1596.2*t-4.905*t^2-8342 ;
diff(s1,t)

ans = 

diff(s2,t)

ans = 

[Sam Chak]

I suppose @Kyle intended to estimate the rate of change based on the simulated data points, rather than directly from the kinematic equation. Observing the curve, it appears that the aircraft (modeled as a point mass) is descending rapidly towards the ground. Is this behavior intentional and meant to be displayed?

[Spaceman]

I was examining the kinematic equations to plot the information and compare it to the derivative approximation of velocity, etc. I was then plotting the integral approximation and comparing it to the past plots. I got it all figured out, though. :)