Why am I getting this interp1 error?
53 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
%% AEM 588 HW1 Question 2
clear all;
close all;
clc;
global A_const;
global B_const;
global C_const;
mu = 398600; %km^3/s^2
% % nu0 = 1e-3;
Isp = 100;
r0 = 2378 + 12800; %km
m0 = 8; %kg
T = 12/1000;
v0 = sqrt(mu/r0);
rho0 =1;
tau0 = 0;
nu0 = T/m0/9.81;
T0 = nu0*(mu*m0)/r0^2*1000.;
mdot = T0/Isp/9.81;
A_const = T0*r0^2/1000;
B_const = mu*m0;
C_const = mdot*sqrt(mu*r0^3);
tspan = [0 10000];
initialConditions = [1.;0.;1.];
%Solve ODE
[tau,rho] = ode45(@odefun,tspan,initialConditions);
t = tau;
r = rho(:,1);
rMoonNorm = ones(length(r),1).*(384400/r0);
b_line = polyfit(t,rMoonNorm,1);
y_line2 = polyval(b_line,t);
x_int = interp1((y_line2-r),t,0);
y_int = polyval(b_line,x_int);
figure(1); clf(1);
plot(t,r)
hold on
plot(t,rMoonNorm)
xlabel('Time, Non-Dimensional')
ylabel('Radius, Non-Dimensional')
t = tau./sqrt(mu/r0^3);
r = rho(:,1).*r0;
rMoon = ones(length(r),1).*384400;
time_intercept = x_int/sqrt(mu/r0^3);
hours_intercept = time_intercept/3600;
days_intercept = hours_intercept/24;
Prop_mass = mdot*time_intercept;
fun = @(x) T0./(m0-mdot*x)
dV = integral(fun,0,time_intercept);
figure(2); clf(2);
plot(t,r)
hold on
plot(t,rMoon)
plot(x_int,y_int,'rx')
xlabel('Time (s)')
ylabel('Radius (km)')
function dydt = odefun(t,y)
global A_const;
global B_const;
global C_const;
dydt = zeros(3,1);
dydt(1) = y(2);
dydt(2) = ((y(3)^2)-y(1))/(y(1)^3);
dydt(3) = (A_const/(B_const - C_const*t))*y(1);
end
2 commentaires
Dyuman Joshi
le 3 Fév 2023
Modifié(e) : Dyuman Joshi
le 3 Fév 2023
Because one of the values in your data is not finite (see below). And interpolation in between infinite set of points is rather un-defined (for the lack of a better word).
%% AEM 588 HW1 Question 2
clear all;
close all;
clc;
global A_const;
global B_const;
global C_const;
mu = 398600; %km^3/s^2
% % nu0 = 1e-3;
Isp = 100;
r0 = 2378 + 12800; %km
m0 = 8; %kg
T = 12/1000;
v0 = sqrt(mu/r0);
rho0 =1;
tau0 = 0;
nu0 = T/m0/9.81;
T0 = nu0*(mu*m0)/r0^2*1000.;
mdot = T0/Isp/9.81;
A_const = T0*r0^2/1000;
B_const = mu*m0;
C_const = mdot*sqrt(mu*r0^3);
tspan = [0 10000];
initialConditions = [1.;0.;1.];
%Solve ODE
[tau,rho] = ode45(@odefun,tspan,initialConditions);
t = tau;
r = rho(:,1);
min(r)
rMoonNorm = ones(length(r),1).*(384400/r0);
b_line = polyfit(t,rMoonNorm,1);
y_line2 = polyval(b_line,t);
min(y_line2-r)
x_int = interp1((y_line2-r),t,0)
y_int = polyval(b_line,x_int);
figure(1); clf(1);
plot(t,r)
hold on
plot(t,rMoonNorm)
xlabel('Time, Non-Dimensional')
ylabel('Radius, Non-Dimensional')
t = tau./sqrt(mu/r0^3);
r = rho(:,1).*r0;
rMoon = ones(length(r),1).*384400;
time_intercept = x_int/sqrt(mu/r0^3);
hours_intercept = time_intercept/3600;
days_intercept = hours_intercept/24;
Prop_mass = mdot*time_intercept;
fun = @(x) T0./(m0-mdot*x)
dV = integral(fun,0,time_intercept);
figure(2); clf(2);
plot(t,r)
hold on
plot(t,rMoon)
plot(x_int,y_int,'rx')
xlabel('Time (s)')
ylabel('Radius (km)')
function dydt = odefun(t,y)
global A_const;
global B_const;
global C_const;
dydt = zeros(3,1);
dydt(1) = y(2);
dydt(2) = ((y(3)^2)-y(1))/(y(1)^3);
dydt(3) = (A_const/(B_const - C_const*t))*y(1);
end
KSSV
le 3 Fév 2023
Check your tspan llimits.....You may reduce it? Also see are all the parameters used are in the same unit system.
Réponses (1)
Bhanu Prakash
le 13 Mar 2023
Hi Shawn,
As per my understanding, you are trying to interpolate a 1-D function using "interp1" function and are facing some errors in it.
If you look at the code, there is a matrix "r" which is passed as an argument to the "interp1" function. The matrix "r" contains integers values along with "NaN" and "Inf (infinity)" values. These values are the reason for the error.
One workaround for this error is to replace the "NaN" and "Inf" values with finite values. I have tried this at my end and the code worked fine.
Hope this answer helps you.
Thanks,
Bhanu Prakash.
0 commentaires
Voir également
Catégories
En savoir plus sur Creating and Concatenating Matrices dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!