How do I interpolate to get the y axis values when I sub in the x axis values? My data is like a plot with no function associated to it.

I have a set of data points (an array, column 1 is for x values and column 2 is for y values) and I want to be able to sub in x values (not part of the initial data points) to get the y values based on interpolation of the current data points I have.

For example I have data for x values [1,3,5,6] and these x values map to y values of [2,4,6,8] but I want to be able to input x value 2 and get a y value based on the interpolation.

How can I do interpolation with a vector set of data like:

final_list(:,1) = x axis values
final_list(:,2) = y axis values

The cubic line from the tools is not able to give me a "best fit curve" hence I cannot use the equation as such. Any ideas how I can use interpolation function from matlab or am I looking at the wrong documentation?

Thanks in advance for any help/advise given!

[Mathieu NOE]

hello

sorry, your request is a bit unclear to me. Do you want to interpolate, curve fit or ???

seems the cubic polynomial fit is doing a good job so what is the issue ?

maybe a piece of code + a picture of what you want is appropriate

[Rachel Ong]

Hello, sorry for the unclear explanation.

The current data set I have is a bit long but let me just put it here for easier reference.

The list that I am referring to is final_list_2_5 and I need to use columns 6 and 2. So col6 is the altitude and col2 is for the elevator deflection. From this array, you will be able to see that the col6 does not have all altitudes eg it jumps from 0 to 90 to 181.

I will want to be able to lets say input altitude = 70 and find the associated value for col2 based on the set of array I have based on interpolation.

Attached is the code. So sorry for not explaining it well.

%% Maximum ROC vs V by taking the maximum ROC at each altitude 

clc
clear all

% ------ Constants ------ %
m = 1248.5;
g = 9.81;
W = m*g;
S = 17.1;

% ------ Functions to solve for FSOLVE (steady climbing flight) ------ %
% func1 = v' = 1/m *(T-D-mgsin(gamma)) = 0 
% func2 = gamma' = 1/mu * (L-mgcos(gamma)) = 0
% func3 = Cm = 0 

% ------ FSOLVE variables ------ %
% x(1) = velocity
% x(2) = deltae
% x(3) = gamma

% ------ Altitude range from flight envelope ------ %
hlist_2_5 = linspace(0,9000,100);
hlist_3 = linspace(0,17000,100);

% ------ Array for storing ------ %
% 1: V   2:deltae   3:gamma   4:aoa   5:roc   6. altitude
final_list_2_5 = zeros(length(hlist_2_5),6);
final_list_3 = zeros(length(hlist_3),6);

%% PS = 2.5
for j = 1:length(hlist_2_5)
    PS = 2.5 
    h = hlist_2_5(j)
    [rho, speedsound] = atmos(h)
    [aoa_start,aoa_end] = aoa_guess(PS)
    aoalist = linspace(aoa_start,aoa_end,150);
    store_v_de_gam_aoa_roc = zeros(length(aoalist),5); % 1: V   2:deltae   3:gamma   4:aoa   5:roc
    
    for i = 1:length(aoalist)
        qbar = @(x) 0.5*rho*x(1)^2*S;
        CL = @(x) 6.44*aoalist(i) + 0.355*x(2);
        L = @(x) CL(x)*qbar(x);
        CD = @(x) 0.03 + 0.05*(CL(x))^2;
        D = @(x) CD(x)*qbar(x);
        Cm = @(x) 0.05 - 0.683*aoalist(i) - 0.923*x(2);
        T = @(x) PS *( (7+x(1)/speedsound)*200/3 + h/1000*(2*x(1)/speedsound -11));

        f1 = @(x) (1/m)*(T(x)-D(x)-m*g*sin(x(3)));
        f2 = @(x) (1/(m*x(1)))*(L(x)-m*g*cos(x(3)));
        f3 = @(x) Cm(x);

        func = @(x) [f1(x);f2(x);f3(x)];

        x = fsolve(@(x) func(x), [40 0 0])
        
        store_v_de_gam_aoa_roc(i,1) = x(1); % V
        store_v_de_gam_aoa_roc(i,2) = x(2); % delta e
        store_v_de_gam_aoa_roc(i,3) = x(3); % gamma
        store_v_de_gam_aoa_roc(i,4) = aoalist(i); % aoa
        
        thrust = PS * ( 200/3*(7 + x(1)/speedsound) + h/1000*(2*x(1)/speedsound - 11)  );
        liftcoeff = 6.44*aoalist(i) + 0.355*x(2);
        drag = 0.5*rho*x(1)^2*S*(0.03 + 0.05*liftcoeff^2);
        ROC = (thrust-drag)*x(1)/W;
        store_v_de_gam_aoa_roc(i,5)= ROC; % ROC
    end    
     
    [max_roc, index] = max(store_v_de_gam_aoa_roc(:,5));
    final_list_2_5(j,5) = max_roc;  % ROC
    final_list_2_5(j,1) = store_v_de_gam_aoa_roc(index,1); % V
    final_list_2_5(j,2) = store_v_de_gam_aoa_roc(index,2); % deltae
    final_list_2_5(j,3) = store_v_de_gam_aoa_roc(index,3); % gamma
    final_list_2_5(j,4) = store_v_de_gam_aoa_roc(index,4) % aoa
    final_list_2_5(j,6) = hlist_2_5(j) % h
    
end



% ------ Plotting graph ------ %
figure     % Elevator deflection vs altitude for PS = 2.5
plot(final_list_2_5(:,6),final_list_2_5(:,2))
xlabel('Altitude')
ylabel('Elevator Deflection')
title('Elevator settings for each altitude at Power setting = 2.5')
legend('PS = 2.5')


%% Atmos function
% Calculating T, p ,rho and speedsound for every altitude in the ISA atmosphere
function [rho, speedsound] = atmos(h)
h1 = 11000;    % Height of tropopause
h2 = 20000;    % End height of table
g = 9.81;
R = 287; 
c = 6.51e-3;       % temperature lapse dt/dh = - c = -6.51 degcelcius/km

T0 = 15+273.15;       % Temperature sea level
p0 = 101325;       % pressure sealevel
rho0 = 101325/R/T0;      % density sealevel = pressure / R*T, R=287, T = 15 degcelcius       

T1 = T0 - c*h1;    % Temperature at 11km
p1 = p0 * (T1/T0)^5.2506;   % Pressure at 11km
rho1 = rho0 * (T1/T0)^4.2506;   % Density at 11km
T2 = T1;    % Temperature at 20km
p2 = p1 * exp(-g/(R*T2)*(h2-h1));  % Pressure at 20km
rho2 = rho1 * exp(-g/(R*T2)*(h2-h1));  % Density at 20km


if h <= h1
%       disp('Troposphere');
    T = T0 - c*h;
    p = p0 * (T/T0)^5.2506;
    rho = rho0 * (T/T0)^4.2506;
    speedsound = (1.4*R*T)^0.5;
elseif h <= h2
%       disp('Tropopause');
    T = T1;
    p = p1 * exp(-g/(R*T)*(h-h1));
    rho = rho1 * exp(-g/(R*T)*(h-h1));
    speedsound = (1.4*R*T)^0.5;
end

return 

end

%% AOA guess function 
% To guess the range of AOA at each altitude
function [aoa_start,aoa_end] = aoa_guess(PS)
h = 0;
W = 1248.5*9.81;        % Weight of UAV
S = 17.1;       % Wing area        

%     CL = (6.44*aoa + 0.355*eledefl);
%     CD = (0.03 + 0.05*(6.44*aoa + 0.355*eledefl)^2);
%     Cm = (0.05 - 0.683*aoa - 0.923*eledefl);

%     thrust = (PS * ( (7 + V/speedsound )*200/3  + h*(2*(V/speedsound) - 11) ));

%%
% V = x(1);
% aoa = x(2);
% eledefl = x(3);
[rho, speedsound] = atmos(h)

qbar = @(x) 0.5*rho*x(1)^2*S
CL = @(x) (6.44*x(2) + 0.355*x(3));
CD = @(x) (0.03 + 0.05*(CL(x))^2);
Cm = @(x) 0.05 - 0.683*x(2) - 0.923*x(3)
thrust = @(x) ( PS * ( 200/3*( 7 + x(1)/speedsound ) + h/1000*( 2*x(1)/speedsound - 11) ) );

func_1 = @(x) qbar(x)*CL(x) + thrust(x)*sin(x(2)) - W;
func_2 = @(x) qbar(x)*CD(x) - thrust(x)*cos(x(2));
func_3 = @(x) Cm(x)

FUNC = @(x) [func_1(x); func_2(x); func_3(x)];

% ------ Guess --> Range of AOA ------ %
Xlow = fsolve(@(x) FUNC(x), [10 0 -30/180*pi]) % high AOA
aoa_end = Xlow(1,2)

Xhigh = fsolve(@(x) FUNC(x), [100 25/180*pi -30/180*pi]) % low AOA
aoa_start = Xhigh(1,2)
end

Answer 1

[Star Strider]

After first creating ‘FinalList’ as a table and writing it to a file:

FinalList = array2table(final_list_2_5, 'VariableNames',{'ROC','V','Delta_e','Gamma','AOA','H'});

pth = pwd;

writetable(FinalList, fullfile(pwd,sprintf('FinalList PS=%.1f.txt',PS)))

fprintf('\nTable written to: %s\n',fullfile(pwd,sprintf('FinalList PS=%.1f.txt',PS)))

Try this —

FinalList = readtable('FinalList PS=2.5.txt');

% ------ Plotting graph ------ %
figure     % Elevator deflection vs altitude for PS = 2.5
plot(FinalList{:,6},FinalList{:,2})
hold on
Alt = [70, 1111, 4225, 6789, 8030];
ElevDefl = @(Alt) interp1(FinalList{:,6}, FinalList{:,2}, Alt);             % Interpolate
plot(Alt, ElevDefl(Alt), 'sr')
hold off
xlabel('Altitude')
ylabel('Elevator Deflection')
title('Elevator settings for each altitude at Power setting = 2.5')
legend('PS = 2.5')

.

[Rachel Ong]

I have done the compilation into a text file and I want to bring it over to another script but it seems like it is giving me an error at a different part of the code after putting the interoplation part into the code. Is there something wrong with what I am doing?

%% ODE45_estimation
clear all
clc
% ------ Constants ------ %
m = 1248.5;
g = 9.81;
W = m*g;
S = 17.1;

% ------ Time interval for simulation ------ %
t_span = [0,3000];         % [start time, end time]

% ------ Initial conditions ------ %
alpha0 = 0.0584957579231716;      % Using trim conditions at sealevel for PS = 2.5
u0 = 55.2682238485092;
q0 = 0;
theta0 = 0.0584957579231716;
mass0 = 1248.5;
h0 = 0;
displacement0 = 0;

init_cond = [alpha0;u0;q0;theta0;mass0;h0;displacement0];

% ------ Solution ------ %
% eledefl_table = zeros(length(t),1)
[t,x] = ode45(@(t,x) odefcn(t,x), t_span, init_cond)
height = x(:,6);
deltae_store = zeros(length(height),1)
for i=1:length(height)
	deltae_store(i,1) = 1.387e-15*height(i)^3 - 1.314e-11*height(i)^2 - 9.405e-07*height(i) - 0.02487;   
end

% ------ Plotting graph ------ %
figure               % alpha
subplot(3,2,1)
plot(t, x(:,1))
xlabel('Time [s]')
ylabel('Angle of attack, alpha [rad]')
title('Angle of attack vs Time')
hold on

subplot(3,2,2)               % velocity                              
plot(t, x(:,2))
xlabel('Time [s]')
ylabel('Velocity [m/s]')
title('Velocity vs Time')
hold on

subplot(3,2,3)               % q
plot(t, x(:,3))
xlabel('Time [s]')
ylabel('Pitch angular velocity, q [rad/s]')
title('Pitch angular velocity vs Time')
hold on

subplot(3,2,4)               % theta
plot(t, x(:,4))
xlabel('Time [s]')
ylabel('Theta [rad]')
title('Theta vs Time')
hold on

subplot(3,2,5)               % mass
plot(t, x(:,5))
xlabel('Time [s]')
ylabel('Mass [kg]')
title('Mass vs Time')
hold on

subplot(3,2,6)               % altitude
plot(t, x(:,6))
xlabel('Time [s]')
ylabel('Altitude [m]')
title('Altitude vs Time')

figure               % trajectory
plot(x(:,7), x(:,6))
xlabel('Horizontal displacement [m]')
ylabel('Altitude [m]')
title('Trajectory of flight')

figure               % delta e
plot(t, deltae_store(:,1))
xlabel('Time [s]')
ylabel('Elevator deflection [rad]')
title('Control - Elevator Deflection')
%% ODE equations function
function dxdt = odefcn(t,x)
    g = 9.81;
    b = 10.2;
    c = 1.74;
    S = 17.1;
    initial_mass = 1248.5;
    Ix = 1421;
    Iy = 4067.5;
    Iz = 4786;
    Ixz = 200;
    csfc = 200/60/60/1000;
    
    PS = 2.5;
    
    alpha = x(1);
    u = x(2); 
    q = x(3);
    theta = x(4);
    mass = x(5);
    height = x(6);
    displacement = x(7);
    
    [rho,speedsound] = atmos(x(6));
    
    % HERE IS THE PART WHERE I AM PUTTING THE INTERPOLATION PART
    % I need to find the deltae at the altitude of x(6) from the ode 
    % Finding delta e from the data of RC max estimation 
    FinalList_2_5 = readtable('FinalList PS=2.5.txt');
    Alt = x(6)
    deltae = @(Alt) interp1(FinalList_2_5{:,6}, FinalList_2_5{:,2}, Alt)
    
%     deltae = 1.387e-15*height^3 - 1.314e-11*height^2 - 9.405e-07*height - 0.02487;    
    
    % EPSILON FIXED AT 0 FIRST 
    epsilon = 0;    
    
    % AERODYNAMICS
    qbar = 0.5*rho*u^2*S;
    CL = 6.44*alpha + 3.8*c*q/(2*u) + 0.355*deltae;
    L = qbar*CL;
    CD = 0.03 + 0.05*CL^2;
    D = qbar*CD;
    Cm = 0.05 - 0.683*alpha - 9.96*c*q/(2*u) - 0.923*deltae;
    m = qbar*c*Cm;
    
    thrust = PS * ((7+u/speedsound)*200/3 + height/1000*(2*u/speedsound - 11) );
    
    % SYSTEM OF ODE
    % dxdt(1) = alpha, dxdt(2) = u, dxdt(3) = q, dxdt(4) = theta, 
    % dxdt(5) = mass, dxdt(6) = height, dxdt(7) = displacement 
    dxdt = [q - (thrust*sin(alpha+epsilon) + L)/(mass*u) + g/u*cos(theta-alpha);
                 (thrust*cos(alpha-epsilon) - D)/mass - g*sin(theta-alpha);
                 m/Iy;
                 q;
                 -csfc*thrust;
                 u*sin(theta-alpha);
                 u*cos(theta-alpha)];
             
end

%% Atmos function
% Calculating T, p ,rho and speedsound for every altitude in the ISA atmosphere
function [rho, speedsound] = atmos(h)
h1 = 11000;    % Height of tropopause
h2 = 20000;    % End height of table
g = 9.81;
R = 287; 
c = 6.51e-3;       % temperature lapse dt/dh = - c = -6.51 degcelcius/km

T0 = 15+273.15;       % Temperature sea level
p0 = 101325;       % pressure sealevel
rho0 = 101325/R/T0;      % density sealevel = pressure / R*T, R=287, T = 15 degcelcius       

T1 = T0 - c*h1;    % Temperature at 11km
p1 = p0 * (T1/T0)^5.2506;   % Pressure at 11km
rho1 = rho0 * (T1/T0)^4.2506;   % Density at 11km
T2 = T1;    % Temperature at 20km
p2 = p1 * exp(-g/(R*T2)*(h2-h1));  % Pressure at 20km
rho2 = rho1 * exp(-g/(R*T2)*(h2-h1));  % Density at 20km


if h <= h1
%       disp('Troposphere');
    T = T0 - c*h;
    p = p0 * (T/T0)^5.2506;
    rho = rho0 * (T/T0)^4.2506;
    speedsound = (1.4*R*T)^0.5;
elseif h <= h2
%       disp('Tropopause');
    T = T1;
    p = p1 * exp(-g/(R*T)*(h-h1));
    rho = rho1 * exp(-g/(R*T)*(h-h1));
    speedsound = (1.4*R*T)^0.5;
end

return 

end

[Star Strider]

I changed the code to read my ‘FinalList PS=2.5.txt’ table online, and corrected the interpolation call.

It needs to be:

Alt = x(6)

ElevDefl = @(Alt) interp1(FinalList_2_5{:,6}, FinalList_2_5{:,2}, Alt); % Interpolate

deltae = ElevDefl(Alt)

However now there are other errors that you need to correct. See the red text below.

(It took too long to run your earlier code online here, limited to 59 seconds, so I ran it offline and created the text file so I could use those results here.)

%% ODE45_estimation

% ------ Constants ------ %
m = 1248.5;
g = 9.81;
W = m*g;
S = 17.1;

% ------ Time interval for simulation ------ %
t_span = [0,3000];         % [start time, end time]

% ------ Initial conditions ------ %
alpha0 = 0.0584957579231716;      % Using trim conditions at sealevel for PS = 2.5
u0 = 55.2682238485092;
q0 = 0;
theta0 = 0.0584957579231716;
mass0 = 1248.5;
h0 = 0;
displacement0 = 0;

init_cond = [alpha0;u0;q0;theta0;mass0;h0;displacement0];

% ------ Solution ------ %
% eledefl_table = zeros(length(t),1)
[t,x] = ode45(@(t,x) odefcn(t,x), t_span, init_cond)

Alt = 0

deltae = -0.0248

Alt = 0

deltae = -0.0248

Alt = -1.9339e-13

deltae = NaN

Alt = -1.3752e-12

deltae = NaN

Output argument "rho" (and possibly others) not assigned a value in the execution with "solution>atmos" function.

Error in solution>odefcn (line 109)
[rho,speedsound] = atmos(x(6));

Error in solution (line 25)
[t,x] = ode45(@(t,x) odefcn(t,x), t_span, init_cond)

Error in ode45 (line 302)
f5 = odeFcn_main(t5, y5);

height = x(:,6);
deltae_store = zeros(length(height),1)
for i=1:length(height)
	deltae_store(i,1) = 1.387e-15*height(i)^3 - 1.314e-11*height(i)^2 - 9.405e-07*height(i) - 0.02487;   
end

% ------ Plotting graph ------ %
figure               % alpha
subplot(3,2,1)
plot(t, x(:,1))
xlabel('Time [s]')
ylabel('Angle of attack, alpha [rad]')
title('Angle of attack vs Time')
hold on

subplot(3,2,2)               % velocity                              
plot(t, x(:,2))
xlabel('Time [s]')
ylabel('Velocity [m/s]')
title('Velocity vs Time')
hold on

subplot(3,2,3)               % q
plot(t, x(:,3))
xlabel('Time [s]')
ylabel('Pitch angular velocity, q [rad/s]')
title('Pitch angular velocity vs Time')
hold on

subplot(3,2,4)               % theta
plot(t, x(:,4))
xlabel('Time [s]')
ylabel('Theta [rad]')
title('Theta vs Time')
hold on

subplot(3,2,5)               % mass
plot(t, x(:,5))
xlabel('Time [s]')
ylabel('Mass [kg]')
title('Mass vs Time')
hold on

subplot(3,2,6)               % altitude
plot(t, x(:,6))
xlabel('Time [s]')
ylabel('Altitude [m]')
title('Altitude vs Time')

figure               % trajectory
plot(x(:,7), x(:,6))
xlabel('Horizontal displacement [m]')
ylabel('Altitude [m]')
title('Trajectory of flight')

figure               % delta e
plot(t, deltae_store(:,1))
xlabel('Time [s]')
ylabel('Elevator deflection [rad]')
title('Control - Elevator Deflection')
%% ODE equations function
function dxdt = odefcn(t,x)
    g = 9.81;
    b = 10.2;
    c = 1.74;
    S = 17.1;
    initial_mass = 1248.5;
    Ix = 1421;
    Iy = 4067.5;
    Iz = 4786;
    Ixz = 200;
    csfc = 200/60/60/1000;
    
    PS = 2.5;
    
    alpha = x(1);
    u = x(2); 
    q = x(3);
    theta = x(4);
    mass = x(5);
    height = x(6);
    displacement = x(7);
    
    [rho,speedsound] = atmos(x(6));
    
    % HERE IS THE PART WHERE I AM PUTTING THE INTERPOLATION PART
    % I need to find the deltae at the altitude of x(6) from the ode 
    % Finding delta e from the data of RC max estimation 
%     FinalList_2_5 = readtable('FinalList PS=2.5.txt');
    FinalList_2_5 = readtable('https://www.mathworks.com/matlabcentral/answers/uploaded_files/913459/FinalList%20PS=2.5.txt');
    Alt = x(6)
    ElevDefl = @(Alt) interp1(FinalList_2_5{:,6}, FinalList_2_5{:,2}, Alt);             % Interpolate
    deltae = ElevDefl(Alt)
    
%     deltae = 1.387e-15*height^3 - 1.314e-11*height^2 - 9.405e-07*height - 0.02487;    
    
    % EPSILON FIXED AT 0 FIRST 
    epsilon = 0;    
    
    % AERODYNAMICS
    qbar = 0.5*rho*u^2*S;
    CL = 6.44*alpha + 3.8*c*q/(2*u) + 0.355*deltae;
    L = qbar*CL;
    CD = 0.03 + 0.05*CL^2;
    D = qbar*CD;
    Cm = 0.05 - 0.683*alpha - 9.96*c*q/(2*u) - 0.923*deltae;
    m = qbar*c*Cm;
    
    thrust = PS * ((7+u/speedsound)*200/3 + height/1000*(2*u/speedsound - 11) );
    
    % SYSTEM OF ODE
    % dxdt(1) = alpha, dxdt(2) = u, dxdt(3) = q, dxdt(4) = theta, 
    % dxdt(5) = mass, dxdt(6) = height, dxdt(7) = displacement 
    dxdt = [q - (thrust*sin(alpha+epsilon) + L)/(mass*u) + g/u*cos(theta-alpha);
                 (thrust*cos(alpha-epsilon) - D)/mass - g*sin(theta-alpha);
                 m/Iy;
                 q;
                 -csfc*thrust;
                 u*sin(theta-alpha);
                 u*cos(theta-alpha)];
             
end

%% Atmos function
% Calculating T, p ,rho and speedsound for every altitude in the ISA atmosphere
function [rho, speedsound] = atmos(h)
h1 = 11000;    % Height of tropopause
h2 = 20000;    % End height of table
g = 9.81;
R = 287; 
c = 6.51e-3;       % temperature lapse dt/dh = - c = -6.51 degcelcius/km

T0 = 15+273.15;       % Temperature sea level
p0 = 101325;       % pressure sealevel
rho0 = 101325/R/T0;      % density sealevel = pressure / R*T, R=287, T = 15 degcelcius       

T1 = T0 - c*h1;    % Temperature at 11km
p1 = p0 * (T1/T0)^5.2506;   % Pressure at 11km
rho1 = rho0 * (T1/T0)^4.2506;   % Density at 11km
T2 = T1;    % Temperature at 20km
p2 = p1 * exp(-g/(R*T2)*(h2-h1));  % Pressure at 20km
rho2 = rho1 * exp(-g/(R*T2)*(h2-h1));  % Density at 20km


if h <= h1
%       disp('Troposphere');
    T = T0 - c*h;
    p = p0 * (T/T0)^5.2506;
    rho = rho0 * (T/T0)^4.2506;
    speedsound = (1.4*R*T)^0.5;
elseif h <= h2
%       disp('Tropopause');
    T = T1;
    p = p1 * exp(-g/(R*T)*(h-h1));
    rho = rho1 * exp(-g/(R*T)*(h-h1));
    speedsound = (1.4*R*T)^0.5;
end

return 

end

.

[Rachel Ong]

Seems like my data has errors. I'll look at it again but thank you for teaching me how to do the interpolation with the data as well as importing it into a text file. Really appreciate! :)

[Star Strider]

As always, my pleasure!