39 views (last 30 days)

Show older comments

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);

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);

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!