# why do i get Nan as output?

39 views (last 30 days)
Sunetra CV on 4 Sep 2019
Commented: Sunetra CV on 16 Sep 2019
im trying to solve four such equations using an input vector containing 4 values.
fn = 1/(rho(1)*Y(1))*(1.113-2.14*Y(1)+2.81*Y(1)^2-1.76*Y(1)^3+0.396*Y(1)^4+0.08694*Y(1)^5-0.0579*Y(1)^6 +0.0077*Y(1)^7);
using ode45.
But the output vector consists of Nan values ; where as its returining numerical values on using the command window.
function file:
function fn = plume5(t,Y)
global rho cp
fn(1) = 1/(rho(1)*cp)*(760.91-1608.08*Y(1)+2449.22*Y(1)^2-2001.61*Y(1)^3+904.51*Y(1)^4-212.087*Y(1)^5+20.095*Y(1)^6);
fn(2) = 1/(rho(2)*cp)*(760.91-1608.08*Y(2)+2449.22*Y(2)^2-2001.61*Y(2)^3+904.51*Y(2)^4-212.087*Y(2)^5+20.095*Y(2)^6);
fn(3) = 1/(rho(3)*cp)*(760.91-1608.08*Y(3)+2449.22*Y(3)^2-2001.61*Y(3)^3+904.51*Y(3)^4-212.087*Y(3)^5+20.095*Y(3)^6);
fn(4) = 1/(rho(4)*cp)*(760.91-1608.08*Y(4)+2449.22*Y(4)^2-2001.61*Y(4)^3+904.51*Y(4)^4-212.087*Y(4)^5+20.095*Y(4)^6);
fn(5) = 1/(rho(5)*Y(5))*(1.113-2.14*Y(5)+2.81*Y(5)^2-1.76*Y(5)^3+0.396*Y(5)^4+0.08694*Y(5)^5-0.0579*Y(5)^6 +0.0077*Y(5)^7);
fn(6) = 1/(rho(6)*Y(6))*(1.113-2.14*Y(6)+2.81*Y(6)^2-1.76*Y(6)^3+0.396*Y(6)^4+0.08694*Y(6)^5-0.0579*Y(6)^6 +0.0077*Y(6)^7);
fn(7) = 1/(rho(7)*Y(7))*(1.113-2.14*Y(7)+2.81*Y(7)^2-1.76*Y(7)^3+0.396*Y(7)^4+0.08694*Y(7)^5-0.0579*Y(7)^6 +0.0077*Y(7)^7);
fn(8) = 1/(rho(8)*Y(8))*(1.113-2.14*Y(8)+2.81*Y(8)^2-1.76*Y(8)^3+0.396*Y(8)^4+0.08694*Y(8)^5-0.0579*Y(8)^6 +0.0077*Y(8)^7);
fn = fn';
main file
clear all
close all
global rho cp
rho = [0.4878;0.7862;1.1449;1.1611;0.4878;0.7862;1.1449;1.1611]
cp = 1.4857
t_span= 0:10
y = [0.0182;0.3105;1.2882;2.7341;0.0182;0.3105;1.2882;2.7341]
[A,P]= ode45('plume5',t_span,y);
##### 2 CommentsShowHide 1 older comment
Sunetra CV on 5 Sep 2019
Thank You

Jan on 4 Sep 2019
Edited: Jan on 4 Sep 2019
Use te debugger to determine the first occurrence of a NaN. Type this in the command window:
dbstop if naninf
Then run the code again. When Matlab stops, check the value of the trajectory. Most likely there is a pole and NaN or Inf is the correct numerical value. You can check this by:
plot(t, Y)
By the way, vectorization can simplify the code:
function fn = plume5(t,Y, rho, cp)
fn = zeros(8, 1);
fn(1:4) = 1 ./ (rho(1:4)*cp) .* (760.91 - 1608.08 * Y(1:4) + ...
2449.22 * Y(1:4) .^ 2 - 2001.61 * Y(1:4) .^ 3 + ...
904.51 * Y(1:4) .^ 4 - 212.087 * Y(1:4) .^ 5 + 20.095 * Y(1:4) .^6 );
fn(5:6) = 1 ./ (rho(5:8) .* Y(5:8)) * (1.113 - 2.14 * Y(5:8) + ...
2.81 * Y(5:8) .^ 2 - 1.76 * Y(5:8) .^ 3 + 0.396 * Y(5:8) .^ 4 + ...
0.08694 * Y(5:8) .^ 5 - 0.0579 * Y(5:8) .^ 6 + 0.0077 * Y(5:8) .^ 7);
Sunetra CV on 16 Sep 2019
As in making the coefficients of y(i) in both the equations of the same range.