Interp1 returning NaN for ode45

Question

Rafael Félix Soriano 2018-11-10

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/429129-interp1-returning-nan-for-ode45

回答： Bruno Luong 2018-11-10

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 个评论
显示 -1更早的评论隐藏 -1更早的评论

Walter Roberson 2018-11-10

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.

请先登录，再进行评论。

请先登录，再回答此问题。