Keep getting the "Array indices must be positive integers or logical values." error

Trying to simulate a water rocket using matlab, however I can't get my for loop to work properly.

This is what I have so far:

close all
clear all

P0 = 6e5; %Pa
Pout = 1e5; %Pa
V0 = 2; %l
V = 1.5; %l
gamma = 1.4;
rho_w = 998;  %kg/m^3

%Pin = P0*(V/V0).^-gamma;

dt = 0.01;
t = 1:dt:100;

Ve(1) = 0; %m/s
Pin(1) = 6e5; %Pa
v0 = 0; %m/s
v = 0; %m/s
M0 = 0.1; %kg
M = M0 + (V0-V)*rho_w;
for i = 1:100
    Pin(i+1) = Pin(i) - P0(V/V0).^-gamma;
    Ve(i+1) = sqrt(2*(Pin(i) - P0)/rho_w);
    %v = Ve*log(M0/M) + v0;
end

Really hoping one of you can enlighten me, because I'm completely lost here.

[Guillaume]

"however I can't get my for loop to work properly"

Rather than leaving it up to us to guess what the problem is. What don't you tell us? If you get an error, give us the full text of the error message. If you don't get the result you expected tell us what result you get and what you expected instead (and why).

The only comments in your code are just to indicate units. There is no comment explaining at all what the code is supposed to be doing or what the variables represent.

Certainly one puzzling thing is that P0(V/V0). What is that supposed to do ? (what it does is try to index the P0 array with the invalid index V/V0 = 0/0 = NaN)

[Image Analyst]

V/V0 = 1.5/2, not 0/0

madhan's Answer below lets the code run without error. No guarantee about the equations correctness though - I didn't try to verify them.

[Guillaume]

Oh! I didn't see that there were both V and V0 variables (uppercase) as well as v and v0 (lowercase).

That's really begging for trouble! You can be sure at some point the wrong case will be used. Never distinguish variable names just by casing.

Answer 1

[madhan ravi]

  Pin(i+1) = Pin(i) - P0*(V/V0).^-gamma;
                        ^----missed it

[Tom Oud]

Whoops! thanks for that! However I still get this error:

"Array indices must be positive integers or logical values.

Error in Water_Rocket (line 23)

Pin(i+1) = Pin(i) - P0*(V/V0).^-gamma;"

Any ideas?

[madhan ravi]

Ofcourse you need to preallocate the arrays before loop: (Wierd I didn't get any errors)

>> clear all
>> close all
clear all
P0 = 6e5; %Pa
Pout = 1e5; %Pa
V0 = 2; %l
V = 1.5; %l
gamma = 1.4;
rho_w = 998;  %kg/m^3
%Pin = P0*(V/V0).^-gamma;
dt = 0.01;
t = 1:dt:100;
Ve(1) = 0; %m/s
Pin(1) = 6e5; %Pa
v0 = 0; %m/s
v = 0; %m/s
M0 = 0.1; %kg
M = M0 + (V0-V)*rho_w;
for i = 1:100
    Pin(i+1) = Pin(i) - P0*(V/V0).^-gamma;
    Ve(i+1) = sqrt(2*(Pin(i) - P0)/rho_w);
    %v = Ve*log(M0/M) + v0;
end
>> Pin
Pin =
   1.0e+07 *
  Columns 1 through 3
   0.060000000000000  -0.029756411635696  -0.119512823271392
  Columns 4 through 6
  -0.209269234907088  -0.299025646542784  -0.388782058178480
  Columns 7 through 9
  -0.478538469814176  -0.568294881449872  -0.658051293085568
  Columns 10 through 12
  -0.747807704721264  -0.837564116356960  -0.927320527992656
  Columns 13 through 15
  -1.017076939628351  -1.106833351264047  -1.196589762899744
  Columns 16 through 18
  -1.286346174535439  -1.376102586171135  -1.465858997806831
  Columns 19 through 21
  -1.555615409442527  -1.645371821078223  -1.735128232713919
  Columns 22 through 24
  -1.824884644349615  -1.914641055985311  -2.004397467621008
  Columns 25 through 27
  -2.094153879256704  -2.183910290892400  -2.273666702528096
  Columns 28 through 30
  -2.363423114163792  -2.453179525799488  -2.542935937435185
  Columns 31 through 33
  -2.632692349070881  -2.722448760706577  -2.812205172342273
  Columns 34 through 36
  -2.901961583977969  -2.991717995613665  -3.081474407249361
  Columns 37 through 39
  -3.171230818885058  -3.260987230520754  -3.350743642156450
  Columns 40 through 42
  -3.440500053792146  -3.530256465427842  -3.620012877063538
  Columns 43 through 45
  -3.709769288699234  -3.799525700334931  -3.889282111970627
  Columns 46 through 48
  -3.979038523606323  -4.068794935242019  -4.158551346877715
  Columns 49 through 51
  -4.248307758513411  -4.338064170149107  -4.427820581784803
  Columns 52 through 54
  -4.517576993420500  -4.607333405056195  -4.697089816691892
  Columns 55 through 57
  -4.786846228327588  -4.876602639963284  -4.966359051598980
  Columns 58 through 60
  -5.056115463234677  -5.145871874870372  -5.235628286506069
  Columns 61 through 63
  -5.325384698141765  -5.415141109777461  -5.504897521413157
  Columns 64 through 66
  -5.594653933048853  -5.684410344684549  -5.774166756320246
  Columns 67 through 69
  -5.863923167955941  -5.953679579591638  -6.043435991227334
  Columns 70 through 72
  -6.133192402863030  -6.222948814498726  -6.312705226134423
  Columns 73 through 75
  -6.402461637770118  -6.492218049405815  -6.581974461041511
  Columns 76 through 78
  -6.671730872677207  -6.761487284312903  -6.851243695948598
  Columns 79 through 81
  -6.941000107584293  -7.030756519219989  -7.120512930855684
  Columns 82 through 84
  -7.210269342491380  -7.300025754127075  -7.389782165762770
  Columns 85 through 87
  -7.479538577398466  -7.569294989034161  -7.659051400669856
  Columns 88 through 90
  -7.748807812305552  -7.838564223941247  -7.928320635576942
  Columns 91 through 93
  -8.018077047212637  -8.107833458848333  -8.197589870484029
  Columns 94 through 96
  -8.287346282119724  -8.377102693755420  -8.466859105391116
  Columns 97 through 99
  -8.556615517026810  -8.646371928662505  -8.736128340298201
  Columns 100 through 101
  -8.825884751933897  -8.915641163569592
>> 

[Tom Oud]

Seems like something went wrong with the i = ... part. When I run it like this with i = 1:100 it works fine, however I would like to be able to set a specific timestep, so it would look something like i = 1:dt:10, with dt being 0,1 for example. When I try this though, I get the same error as before. Weird right? Or am I missing something?

[madhan ravi]

First of all you need to know in Matlab index should be an integer , try the below:

x=2;
x(1)
x(0.1) -> you will get error here

[Tom Oud]

Nvm, I'm an idiot. i has to be an integer, just like the error stated.

Thanks for the help!

[madhan ravi]

Nobody's perfect. If my answer helped you make sure to accept the answer.

Anytime :)

[Steven Lord]

Since from your code I'm guessing you want values of Pin and Ve corresponding to the elements of your t vector

t = 1:dt:100;

loop from 1 to the number of elements in t then use the corresponding element in t to create the corresponding element in Pin and Ve.

P = zeros(size(t));
for whichtime = 1:numel(t)
    thetime = t(whichtime);
    P(whichtime) = thetime.^2;
end
plot(t, P)

Yes, there's a simpler way to achieve that above example than a for loop. But I think this demonstrates the general technique clearly.

[Jan]

@Tom Oud: If you need a loop over non-integer values, care for integer indices:

v = 1:0.1:10;
for k = 1:numel(v)
    a = v(k);
    Pin(k) = a;
end

Now you have a fractional a, but an integer k.

[Tom Oud]

I tried to apply this in my script, however now when I change the timestep, the values I get from my loop just get smaller, instead of the timestep. So when I change timestep from 0,1 to 0,01 Pin for example changes from 6 to 0,6. Any clue what might cause that?

[Steven Lord]

Your original loop code doesn't show either Pin or Ve depending on either t or dt, so I'm not sure. Can you show your updated code?

[Tom Oud]

Of course!

P0 = 6e5; %Pa

Pout = 1e5; %Pa

V0 = 2; %l

V = 1.5; %l

gamma = 1.4;

rho_w = 998; %kg/m^3

dt = 0.001;

t = 1:dt:10;

Ve(1) = 0; %m/s

Pin(1) = 6e5; %Pa

v0 = 0; %m/s

v = 0; %m/s

M0 = 0.1; %kg

M = M0 + (V0-V)*rho_w;

for i = 1:numel(t)

k = t(i);

Pin(i+1) = Pin(i) - P0*(V/V0).^-gamma;

Ve(i+1) = sqrt(2*(Pin(i) - Pout)/rho_w);

if Pin(i) < P0

Ve(i+1) = Ve(i);

end

%v = Ve*log(M0/M) + v0;

end

[Tom Oud]

Turns out I interpreted my results wrong. Decreasing the timestep now just seems to add more steps in total, instead of making the increment smaller.