Interp1 returning NaN for ode45
4 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Rafael Félix Soriano
el 10 de Nov. de 2018
Respondida: Bruno Luong
el 10 de Nov. de 2018
The interp1 function returns a NaN value for the time-dependent functions that are later introduced in the ode45 function. Using pchip or spline instead of interp1 returns a value other than NaN, but the solutions are not valid.
My code is the one that follows, thank you in advance.
function [T_V , W_V , x_V , y_V , h_V , Vfin] = viraje(x0 , y0 , h0 , W0 ,...
V0 , Vf , gamma0 , gammaf , chi0 , chif , chip , CD0 , k , S , cE)
clc
g = 9.80665; %m/s2
tfin = (chif - chi0)/chip;
gammap = (gammaf-gamma0)/tfin;
Vp = (Vf-V0)/tfin;
t_interp = linspace(0 , tfin , 10);
V = V0 + Vp*t_interp;
gamma = gamma0 + gammap*t_interp;
mu = atan(chip/g*V.*cos(gamma)./(V/g*gammap + cos(gamma)));
h = Hviraje(t_interp , h0 , V0 , Vp , gamma0 , gammap);
rho = densidad(h);
function dWdt = dif(t,W,t_interp,V,gamma,mu,rho)
if Vp == 0
V = V0;
else
V = interp1(V,t_interp,t,'pchip');
end
if gammap == 0
gamma = gamma0;
else
gamma = interp1(gamma,t_interp,t,'pchip');
end
if Vp == 0 && gammap == 0
mu = mu(1);
else
mu = interp1(mu,t_interp,t,'pchip');
end
if gamma0 == 0 && gammap == 0
rho = densidad(h0);
else
rho = interp1(rho,t_interp,t,'pchip');
end
dWdt = -cE*(0.5*rho.*(V.^2)*S*CD0 + 2*k./(rho*S)*(chip*W.*cos(gamma)./...
(g*sin(mu))).^2 + W.*sin(gamma) + Vp/g*W)
end
tspan = [0 tfin];
ic = W0;
[t,W] = ode45(@(t,W) dif(t,W,t_interp,V,gamma,mu,rho) , tspan , ic);
W_V = W;
h_V = Hviraje(t , h0 , V0 , Vp , gamma0 , gammap);
x_V = Xviraje(t , x0 , V0 , Vp , gamma0 , gammap , chi0 , chip);
y_V = Yviraje(t , y0 , V0 , Vp , gamma0 , gammap , chi0 , chip);
rho = densidad(h_V)';
gamma = gamma0 + gammap*t;
V = V0 + Vp*t;
mu = atan(chip/g*V.*cos(gamma)./(V/g*gammap + cos(gamma)));
T_V = 0.5*rho.*V.^2*S*CD0 + 2*k./(rho*S).*(chip*W.*cos(gamma)./...
(g*sin(mu))).^2 + W.*sin(gamma) + Vp/g*W;
Vfin = Vf;
end
function H = Hviraje(T , h0 , V0 , Vp , gamma0 , gammap)
syms h(t)
eqn = diff(h,t) == (V0 + Vp*t)*sin(gamma0 + gammap*t);
cond = h(0) == h0;
ySol(t) = dsolve(eqn,cond);
H = ySol(T);
end
function X = Xviraje(T , x0 , V0 , Vp , gamma0 , gammap , chi0 , chip)
syms x(t)
eqn = diff(x,t) == (V0 + Vp*t)*cos(gamma0 + gammap*t)*cos(chi0 + chip*t);
cond = x(0) == x0;
ySol(t) = dsolve(eqn,cond);
X = ySol(T);
end
function Y = Yviraje(T , y0 , V0 , Vp , gamma0 , gammap , chi0 , chip)
syms y(t)
eqn = diff(y,t) == (V0 + Vp*t)*cos(gamma0 + gammap*t)*sin(chi0 + chip*t);
cond = y(0) == y0;
ySol(t) = dsolve(eqn,cond);
Y = ySol(T);
end
1 comentario
Walter Roberson
el 10 de Nov. de 2018
interp1() returns nan for most interpolation methods when you ask it to project at locations that are outside the range of the first parameter, unless you pass it the 'extrap' parameter.
Respuesta aceptada
Bruno Luong
el 10 de Nov. de 2018
V = interp1(V,t_interp,t,'pchip');
It looks arguments are swapped; I think it should be
V = interp1(t_interp,V,t,'pchip');
and similar for other INTERP1 calls.
0 comentarios
Más respuestas (0)
Ver también
Categorías
Más información sobre Dimensionality Reduction and Feature Extraction en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!