# How do I interpolation into a ODE45 function for my trajectory code?

4 views (last 30 days)
Dylan Allen on 30 Jan 2023
Edited: Torsten on 31 Jan 2023
I'm a student and still learning MATLAB, I am trying to create a code to simulate the a falling payload using ODE 45 with two differnatial equations form velcoity and height. I have this section of the code working the trouble is when I try to also get Mach number as a function of time. This involves getting speed of sound which is a fucntion of temperature. I want to interpolate the temperature from the standard atmosphere and this is where I am running ito issue. Inside of the ODE function I have to create a temperature funciton to get the T value and I currently can get the velocity and height out of my ODE 45, but am unsure how to get the temperature and speed of sound into the function when it is needed to be interpolated so I'm unsure how to proceed. I have tried moving the entire interpolation into the fucntion, but it doesnt seem to work. The code is provided below.
Any feedback is helpful, thank you!!
This is my first time working with ODE 45 and simulation through MATLAB so I imagine there is a number of things that can be improved upon.
----------------------------------------------------------------------------------------------------------
clc
clear
%Defining temperature and speed of sound
T = 273.15+([15,8.5,2,-4.49,-10.98,-17.47,-23.96,-30.45,-36.94,-43.42,-49.90,-56.50,-56.50,-51.60,-46.64]);
a = sqrt(1.4*287.*T);
% define the initial conditions for the falling object
y0 = 30000; % initial height, in meters
v0 = 0; % initial velocity, in meters/second
% define a function that returns the derivative of the state vector
% the state vector contains the position and velocity of the object
% use ode45 to solve the differential equation for the falling object
tspan = linspace(0,150,30001); % time span for the simulation, in seconds
yinit = [y0; v0]; % initial state vector
[t, y] = ode45(@falling_object, tspan, yinit);
% unpack the position and velocity data from the solution
h = y(:,1);
V = y(:,2);
% plot the position and velocity of the falling object
figure;
plot(t, h, 'r-'); % position vs time
title('Position vs Time Cd = 0.1')
xlabel('Time (s)');
ylabel('Position (m)');
ylim([0 30000]) figure;
plot(V, h, 'r-',a, h1, 'k--'); % velocity vs position
Unrecognized function or variable 'h1'.
title('Velocity vs Position Cd = 0.1')
xlabel('Velocity (m/s)');
ylabel('Position (m)');
ylim([0 30000])
legend('Projected Velocity Profile', 'Speed of Sound')
% figure;
% plot(M, h, 'k-'); % Mach number vs Position
% xlabel('Mach Number (m/s)');
% ylabel('Position (m)');
function dydt = falling_object(t, y)
% unpack the state vector
h = y(1);
V = y(2);
% define the acceleration of the object due to gravity
g = 9.8; % in meters/second^2+
%Define the density as using logarithmic interpolation
% nrho = [1.225,1.112,1.007,.9093,.8194,.7364,.6601,.5900,.5258,.4671,.4135,.1948,.08891,.04008,.01841];
% N = linspace(0,30000,30001);
% lnrho = log(nrho);
% ilnrho = interp1(ypos,lnrho , N);
% rho = exp(ilnrho);
rho0 = 1.125; %Define density at sea level
R = 287;
if h < 25000
T = 273-131.21+ h.*0.00299;
elseif h >= 25000
T = 273-56;
end
%Define Ballistic Coefficenient
m = 9; %Mass inkg
W = m*g; %Weight in N
A = 0.09; %Cross sectional area in m^2
cd = .2; %Drag coefficent
B = W/(cd*A); %Ballistic Coefficient
% calculate the derivatives of the position and velocity
dydt = [-V; g-(.5*((rho0)*exp(-g*h/(R*T))*(V^2)*(g))/(B))];
%g-1/(2*W)*cd*A*rho0 * exp(-g*h/R*T)*V^2];
end
Dylan Allen on 30 Jan 2023
Yes, the temperature variable is wihin the dydt function, how would I be able to plot the mach number with respect to position as I have done with velocity?

Torsten on 31 Jan 2023
Edited: Torsten on 31 Jan 2023
I additionally told ode45 to stop integration at h = 0.
clc
clear
% define the initial conditions for the falling object
y0 = 30000; % initial height, in meters
v0 = 0; % initial velocity, in meters/second
% define a function that returns the derivative of the state vector
% the state vector contains the position and velocity of the object
% use ode45 to solve the differential equation for the falling object
tspan = linspace(0,150,30001); % time span for the simulation, in seconds
yinit = [y0; v0]; % initial state vector
options = odeset('Events',@myEvent);
[t, y] = ode45(@falling_object, tspan, yinit,options);
% unpack the position and velocity data from the solution
h = y(:,1);
V = y(:,2);
R = 287;
kappa = 1.4;
for i = 1:numel(h)
if h(i) < 25000
T = 273-131.21+ h(i)*0.00299;
elseif h(i) >= 25000
T = 273-56;
end
M(i) = V(i)/sqrt(kappa*R*T);
end
figure
plot(M, h, 'k-'); % Mach number vs Position
xlabel('Mach Number (1)');
ylabel('Position (m)'); function dydt = falling_object(t, y)
% unpack the state vector
h = y(1);
V = y(2);
% define the acceleration of the object due to gravity
g = 9.8; % in meters/second^2+
%Define the density as using logarithmic interpolation
% nrho = [1.225,1.112,1.007,.9093,.8194,.7364,.6601,.5900,.5258,.4671,.4135,.1948,.08891,.04008,.01841];
% N = linspace(0,30000,30001);
% lnrho = log(nrho);
% ilnrho = interp1(ypos,lnrho , N);
% rho = exp(ilnrho);
rho0 = 1.125; %Define density at sea level
R = 287;
if h < 25000
T = 273-131.21+ h.*0.00299;
elseif h >= 25000
T = 273-56;
end
%Define Ballistic Coefficenient
m = 9; %Mass inkg
W = m*g; %Weight in N
A = 0.09; %Cross sectional area in m^2
cd = .2; %Drag coefficent
B = W/(cd*A); %Ballistic Coefficient
% calculate the derivatives of the position and velocity
dydt = [-V; g-(.5*((rho0)*exp(-g*h/(R*T))*(V^2)*(g))/(B))];
%g-1/(2*W)*cd*A*rho0 * exp(-g*h/R*T)*V^2];
end
function [value,isterminal,direction] = myEvent(t,y)
value = y(1);
isterminal = 1;
direction = 0;
end

### Categories

Find more on Programming in Help Center and File Exchange

### Community Treasure Hunt

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

Start Hunting!