while funtion, matrix and vectors behaving oddly

Below is a code I made that behaves oddly. I am using Matlab R2014B on a Mac with OSX Yosemite V. 10.10.1. It is a pretty straight-forward while function that saves a number of answers in vector format and then sums all answers up in a matrix (and plots some) at the end of the code. The code does not produce any error messages. Here are the odd parts: when the program writes the matrix at the end ('answers'), it ends up being pretty much all zeroes. Meanwhile, writing out individual vectors that make up the matrix (for example 'Avector') show that most do indeed contain numbers. In fact, if you remove the vectors Cevector and Ctvector (where the error probably lies) from the answer matrix, it stops showing the first two columns/vectors of the matrix as zeroes. Also, the minimum value of Ctvector that the code writes out at the end is clearly not 0, even though the vector only contains zeroes when you write it out. I am trying to grasp why matlab changes the values of these vectors to zeroes under these conditions...

clear all;

%assignment is to find the width b that minimizes total costs Ct
b=0; % start width value [m]
y=2.5; % [m]
L=1; % [m]
M=33; % [m^(1/3)/s]
Q=35; % [m^3/s] mean discharge

T=5000; % [hours/year] operation time of power plant
excost=50; % [USD/m^3] excavation costs
valenergy=0.04; % [USD/kWh] value of produced energy
rates=0.08; % [percent/year] interest rates, repayment, maintenance

% Ct=Ce(value of energy losses) + Cc(capital costs)
%value of energy losses Ce=W * Value of produced energy
%Capital cost Cc=interest rate * excavated volume * excavation costs

% energy loss is W = P*T = 8.2*Q*hfr*T [kWh]
%hfr=((U^2)*L)/((M^2)*(R^(4/3)))

bvector=[];
Avector=[];
Ccvector=[];
Cevector=[];
Ctvector=[];
Pvector=[];
Rvector=[];
Uvector=[];
hfrvector=[];
Wvector=[];

while b<=50;
    A=b*y; % cross section area [m^2]
    P=b+(2*y); % wetted perimeter [m]
    R=A/P; % hydraulic radius [m]
    U=Q/A; % [m/s]
    hfr=((U^2)*L)/((M^2)*(R^(4/3))); % [m] friction losses

      W=8.2*Q*hfr*T; % [kWh] energy loss
      Cc=rates*A*L*excost; % [USD] capital costs, this is per m length!!! not whole channel
      Ce=W*valenergy; % [USD] value of energy losses
      Ct=Cc+Ce; % [USD] total costs

      bvector=[bvector; b];
      Avector=[Avector; A];
      Pvector=[Pvector; P];
      Rvector=[Rvector; R];
      Uvector=[Uvector; U];
      hfrvector=[hfrvector; hfr];
      Wvector=[Wvector; W];
      Ccvector=[Ccvector; Cc];
      Cevector=[Cevector; Ce];
      Ctvector=[Ctvector; Ct];

      b=b+0.01; 

end

answers=[Avector, Ccvector, Cevector, Ctvector]
disp('A, Cc, Ce, Ct')

Ctmin=min(Ctvector)

plot(Avector, Ccvector, 'r', Avector, Cevector, 'g', Avector, Ctvector, 'b') %three plots in one graph
title('Costs and Cross-sectional area')
xlabel('Area [m^2]')
ylabel('Cc(r), Ce(g), Ctotal (b) [USD]')

Answer 1

[Star Strider]

See the documentation for the format function. I believe that’s the problem. (I use ‘format short eng’.)

This is what I get for the first few rows of ‘answers’:

answers =

       0.0000e+000     0.0000e+000             Inf             Inf
      25.0000e-003   100.0000e-003   121.1538e+009   121.1538e+009
      50.0000e-003   200.0000e-003    12.0520e+009    12.0520e+009
      75.0000e-003   300.0000e-003     3.1278e+009     3.1278e+009
     100.0000e-003   400.0000e-003     1.2021e+009     1.2021e+009
     125.0000e-003   500.0000e-003   572.8513e+006   572.8513e+006
     150.0000e-003   600.0000e-003   312.7879e+006   312.7879e+006
     175.0000e-003   700.0000e-003   187.6019e+006   187.6019e+006

The range of the content is on the order of 1E+12 (ignoring the Inf entries). You’re likely not seeing everything that’s in your matrix because of the scale range.

EDIT —

I would also change the plot to:

semilogy(Avector(2:end), Ccvector(2:end), 'r', Avector(2:end), Cevector(2:end), 'g', Avector(2:end), Ctvector(2:end), 'b') %three plots in one graph
title('Costs and Cross-sectional area')
xlabel('Area [m^2]')
ylabel('Cc(r), Ce(g), Ctotal (b) [USD]')
grid
legend('Ccvector', 'Cevector', 'Ctvector', 'Location', 'NE')