ode45 not working (giving non real values)

1 view (last 30 days)
Why am I getting non real values for M.
iFunction file
function dMdx=s(x,M,gamma_b,a1Vec,b1Vec,a11Vec,b11Vec,xVec,xMin)
a1=interp1(xVec,a1Vec,x);
b1=interp1(xVec,b1Vec,x);
a11=interp1(xVec,a11Vec,x);
b11=interp1(xVec,b11Vec,x);
if xVec==xMin
dMdx=(1/4)*(-gamma_b*a1+sqrt((gamma_b*a1)^2-4*((gamma_b-1)*a1^2-(1+gamma_b)*a11+(((1+gamma_b)^2)/2)*b11)));
else
dMdx= M*((1+((gamma_b-1)/2)*M^2)/(1-M^2))*(-a1+((1+(gamma_b*M^2))/2)*b1);
end
end
iScript
xVec=0:0.01:0.4; %axial combuster length (can be varied according to user input and units)
xi=xVec(1); %station 3 (combuster starting point)
x4=xVec(end);%station 4 (combuster end point)
gamma_b=1.36;%Average value of ratio of specific heat capacities during burning
gamma_c=1.4;%Average value of ratio of specific heat capacities during compression
X=((xVec-xi)/(x4-xi));
a3=0.0038;%station 3 cross secion area in m2.
A=a3*(1+X.^2);%Area Profile
a=gradient(A,xVec);%dA/dx in equation (6-90)
a1Vec=a./A; %(dA/dx)/A in equation (6-90)
d2a=gradient(a,xVec);
a11Vec=d2a./A;
tau_b=1.5; %Tt4/Tt2 in eq. (6-91) overall total temerature rise in burner
theta=2;%emperical constant in eq. (6-91)
Tt2=2200;%Total temperature at station 2 depends on inlet and compression conditions (USER SPECIFIED)
Ttx=Tt2*(1+(tau_b-1)*((theta*X)./(1+(theta-1)*X)));%Temperature profile eq. (6-91)
b=gradient(Ttx,xVec);%dTt/dx in equation (6-90)
b1Vec=b./Ttx;%(dTt/dx)/Tt in equation (6-90)
d2b=gradient(b,xVec);
b11Vec=d2b./Ttx;
fun = @(xVec)((a./A)-(((1+gamma_b)/2)*(b./Ttx))); % Function of x only
[fMin, idx] = min(abs(fun(xVec)));
xMin = xVec(idx);
initial_M=1;
[x,M]=ode45(@(x,M) s(x,M,gamma_b,a1Vec,b1Vec,a11Vec,b11Vec,xVec,xMin),[xMin xi],initial_M);%Solve ODE eq. (6-90)

Accepted Answer

Walter Roberson
Walter Roberson on 20 Jul 2015
Your test
if xVec == xMin
is testing a vector of values against a single value. Your xVec are distinct so at most once of those will compare equal, so you will get a logical vector that has at most one true. When you have an "if" that is being asked to check a vector, the result is considered true only if all of the values in the vector are true. That can never be the case for you so the "else" is always taken.
You probably want to test x == xMin instead of xVec == xMin. But if you do that you should probably keep in mind that it is not easy to establish conditions under which floating point numbers that are intended to be the same value will be bit-for-bit identical and so end up comparing true. Your code as it stands does use the right conditions so if you do not change it then if you compare x == xMin then it will compare true the first time through, but it is probably not a good idea to count on that.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!