Index exceeds matrix dimensions.

clear;
%Input data
Q=22;
Z=2.5;
b=7.5;
L=3000;
So=0.0006;
n=0.015;
g=9.81;
%Set range and increments of y
y=0.01:.001:5;
%Define Functions
A=b*y+Z*y.^2;
B=b+2*Z*y;
P=b+2*y*sqrt(1+Z^2);
S=Q^2*n^2./(A.^3.3333.*P.^-1.3333);
Fr=sqrt(Q^2*B./(g*A.^3));
%Find values of y where (s-So) and (Fr-1) cross zero
[S,t]=min(abs(S-So));
yN=y(t)
[S,t]=min(abs(Fr-1));
yc=y(t)
%Reset range and increments of y to go from critical
%to normal 
clear y;
dy=0.1;
y=yc:dy:yN
%Calculate specific energy at new y values
A=b*y+Z*y.^2;
E=y+Q^2./(2*g*A.^2)
%Calculate x locations according to Eq.10.7.6
ym=(y(2:5)+y(1:4))/2;
A_ym=b*ym+Z*ym.^2;
P_ym=b+2*ym*sqrt(1+Z^2);
S_ym=Q^2*n^2./(A_ym.^(10/3).*P_ym.^(-4/3));
x(1)=L;
for i=2:5;
x(i)=x(i-1)+(E(i)-E(i-1))/(So-S_ym(i-1));
end;

[Walter Roberson]

Please post the complete error message including the line that the problem is happening on.

[Shifo]

Index exceeds matrix dimensions.

Error in Q3 (line 33) ym=(y(2:100)+y(1:99))/2;

Here is the error

Answer 1

[Walter Roberson]

Your yN and yc are extracted from the 101-element long y=0:.1:10 , so if you were to do yc:.1:yN the absolute best you could hope for would be 101 elements but because of the way you determine yc and yN, yc:.1:yN would almost certainly be shorter than 100. But you do not use .1 as the increment: you do yc:.2:yN so even if yc turned out to be y(1) and yN turned out to be y(end), the longest that yc:.2:yN could possibly be would be 51 elements. You then try to index that "at most 51 long" vector at element number 100.

Answer 2

[dpb]

...
yN =
  0.4000
yc =
  0.3000
y =
  0.3000
E =
  0.3844
Index exceeds matrix dimensions.

You've got all the debugging info you need; you defined

y=yc:dy:yN;

which is the same as

y=0.3:0.2:0.4;

so as the debugging info you printed shows, y turns out to be just yc because

yc+dy > yN

so the expansion is satisfied after the first element.

>> whos y
Name      Size            Bytes  Class     Attributes

y         1x1                 8  double

>>

You subsequently expect to have at least 200 elements in y.

Not sure just what you really are intending here, but that's why it fails. How to set the two limits to what really want or whether those are right but you're just expecting something different and dy should be much smaller or just not have as many elements is something you'll have to decide based on what it is you're actually trying to do.

[Shifo]

Thanks a lot, Okay what I need to do to make it even 50. what i really want is to guess y at different x say dx is 30 m so total we have 3000/30 is 100 location of x which should have 100 value of y and y is increasing with the increasing with distance.

[Shifo]

The code down is perfectly working, however i want get values for 50 values for y at least at different locations of x, I have tried to change i but its not working and giving matrix error as posted. What i really want is to increase the matrix to give values at the different locations

clear;
%Input data
Q=22;
Z=2.5;
b=7.5;
L=3000;
So=0.0006;
n=0.015;
g=9.81;
%Set range and increments of y
y=0.01:.001:5;
%Define Functions
A=b*y+Z*y.^2;
B=b+2*Z*y;
P=b+2*y*sqrt(1+Z^2);
S=Q^2*n^2./(A.^3.3333.*P.^-1.3333);
Fr=sqrt(Q^2*B./(g*A.^3));
%Find values of y where (s-So) and (Fr-1) cross zero
[S,t]=min(abs(S-So));
yN=y(t)
[S,t]=min(abs(Fr-1));
yc=y(t)
%Reset range and increments of y to go from critical
%to normal 
clear y;
dy=0.1;
y=yc:dy:yN
%Calculate specific energy at new y values
A=b*y+Z*y.^2;
E=y+Q^2./(2*g*A.^2)

%Calculate x locations according to Eq.10.7.6
ym=(y(2:5)+y(1:4))/2;
A_ym=b*ym+Z*ym.^2;
P_ym=b+2*ym*sqrt(1+Z^2);
S_ym=Q^2*n^2./(A_ym.^(10/3).*P_ym.^(-4/3));
x(1)=L;
for i=2:5;
x(i)=x(i-1)+(E(i)-E(i-1))/(So-S_ym(i-1));
end;

[Shifo]

Any help guys?

[dpb]

The underlying real problem definition is still lacking, but the first step could be

function y=shifo
Q=22;
Z=2.5;
b=7.5;
L=3000;
So=0.0006;
n=0.015;
g=9.81;

function S=fnS(y)
  A=b*y+Z*y.^2;
  P=b+2*y*sqrt(1+Z^2);
  S=Q^2*n^2./(A.^3.3333.*P.^-1.3333);
end

function Fr=fnFr(y)
  A=b*y+Z*y.^2;
  B=b+2*Z*y;
  Fr=sqrt(Q^2*B./(g*A.^3));
end

%Find values of y where (s-So) and (Fr-1) cross zero
yS=fzero(@(y) fnS(y)-So,1);
yF=fzero(@(y) fnFr(y)-1,1);
y=[yS;yF];
end

where I've used internal functions so the two calculating functions have access to all the constants without passing them directly or duplicating them in code; in later releases you can also include functions in a script file; can't yet w/ R2014b which is what I have...anyway, those logistics aside, the above finds the two y locations for the two functions directly -- the question is precisely what to do next. I'll note part of the difficulty in your original code is the reuse of y for different variables; keep the two solutions separate until you're done with them.

The above yields--

>> shifo
ans =
  1.2916
  0.8648
>>

for yS, yFr respectively.

Is there supposed to be a minimum or something else somewhere between these two points I gather? If so, what is it specifically you're trying to solve for next as end result?

Answer 3

[Image Analyst]

Maybe you should learn about and use the linspace() function instead of using the colon operator.

[dpb]

Fixed typo of 'colon' for 'color', IA, but while linspace could get him past his "number of elements in a vector" problem, the actual solution to his problem is to use fzero to find the original solutions then do something similar probably to find the final result he's looking for. That may be, if the crystal ball is functioning correctly this morning, an inverse interpolation problem...

Unfortunately, he's so far not yet shown us the paper/text from which Eq. 10.7.6 comes so we don't know the end objective (and I don't recognize the correlation from the equations presented so far, anyway).

[Shifo]

Hi dpb, thank you so much for your help, here is equation

[Shifo]

Here is the another a code which is really what i want, I have found this online but cant understand it or even use it

                    Q=22;
                    Z=3;
                    b=25;
                    L=3000;
                    So=0.001;
                    n=0.02;
                    g=9.81;
% ---------------------------------------------------------------------------
                    yy =yd(e);
                    zz = eval(['z',int2str(e),'(1);']); 
                    zz1 = zz-s(e)*L(e);
                    mm = eval(['m',int2str(e),'(1);']); 
                    mm1 = mm+L(e);
                    eval(['y',int2str(e),'(p(e)+1)','= yy;']) 
                    eval(['m',int2str(e),'(p(e)+1)','= mm1;']); 
                    eval(['z',int2str(e),'(p(e)+1)','= zz1;']);                     
                    A = (b(e)+zt(e)*yy)*yy;
                    pe = b(e)+2*yy*(1+zt(e)^2)^(1/2);
                    v2 = (Q/A)^2/(2*9.8);
                    ar= 1/n(e)*A^(5/3)/pe^(2/3);
                    E2 = zz1+yy*ss(e)+v2;
                    sf2=(Q/ar).^2;
                    % intial h
                    h = 2;
                    while h <= p(e)+1
                    yy1 = eval(['y',int2str(e),'(p(e)+3-h);']);   
                    mm2 = mm+(p(e)+1-h)*dx(e);
                    zz2 = zz-s(e)*(mm2-mm);
                    eval(['m',int2str(e),'(p(e)+2-h)','= mm2;']); 
                    eval(['z',int2str(e),'(p(e)+2-h)','= zz2;']); 
   for i = 1:num_itera;
   A = (b(e)+zt(e)*yy1)*yy1;
   pe = b(e)+2*yy1*(1+zt(e)^2)^(1/2);
   R = A/pe;
   T = b(e)+2*zt(e)*yy1;
   v2 = (Q/A)^2/(2*9.8);
   ar= 1/n(e)*A^(5/3)/pe^(2/3);
   E1 = zz2+yy1*ss(e)+v2;
   sf1=(Q/ar).^2;                            
   f = E2-E1+0.5*(sf1+sf2)*dx(e); %Aritmetic mean friction slope
   df = -(1-Q^2/9.8*T/A^3)+dx(e)/2*Q^2*n(e)^2*(-10/3*pe^(4/3)*A^(-13/3)*T+(4/3)*A^(-10/3)*pe^(1/3)*(2*(1+zt(e)^2)^0.5));
          yy1= yy1-f/df;                             
          if (abs(f/df)/(yy1 - 0.5*f/df)<tol)
          evalc(['y',int2str(e),'(p(e)+2-h)','= yy1;']); 
                 break
                        end
                            end
           sf2=sf1;
           E2=E1;
           h=h+1;
           end                    
           yy5 = eval(['y',int2str(e),'(1);']);
           A = (b(e)+zt(e)*yy5)*yy5;
           v2 = (Q/A)^2/(2*9.8);
           zz5 = eval(['z',int2str(e),'(1);']); 
           EE5 = yy5*ss(e)+v2+zz5;
           eval(['EN',int2str(e),'(1)','= EE5;']);                

           yy6 = eval(['y',int2str(e),'(p(e)+1);']); 
           A = (b(e)+zt(e)*yy6)*yy6;
           v2 = (Q/A)^2/(2*9.8);
           zz6 = eval(['z',int2str(e),'(p(e)+1);']); 
           EE6 = yy6*ss(e)+v2+zz6;
           eval(['EN',int2str(e),'(2)','= EE6;']);

[dpb]

Well, you stopped at the intriguing part of "the evaluation of conditions a location i+1 proceeds as follows:"

Looks like there would be the key to how the text did it.

I'd think if knew just what/how your variables coincided with those in the text that one could use the diff eq solvers in Matlab to solve 10.7.4 directly and then compute the desired quantities from that solution but I admit without some more context I have no idea what this all is about or what you're trying to achieve as an end result so not sure where to go, precisely.

What are the two y values we worked so hard to compute above and how do they bear on any of the above???

[Stephen23]

"Here is the another a code which is really what i want, I have found this online but cant understand it or even use it"

That seems self-contradictory, unless what you want is code that you cannot use. In any case, see also this question:

https://www.mathworks.com/matlabcentral/answers/358148-any-help-how-i-can-use-this-code-to-produce-result

[Shifo]

I am trying to use the netwotn raphson method to do iterative steps to find values of y along x direction, consider it as you have a river with flow and you know y " flow depth" and Q the discharge downstream so as shown above this is used to integrated to find the flow depth upward, I have tried everything. the second code seems to do the job in other text book, the first one here is okay to integrate from 0 to 5 distant upward.

[Shifo]

Hi dpb, the two y above is the normal and critical depth

[Stephen23]

Someone wrote a Newton-Raphson root finder with fifteen calls to eval ?! As a joke, I presume.

@Shifo: your posts and questions are a huge waste of your own time. After several days of posting useless code and random textbook excerpts, we are still none-the-wiser to what you are actually trying to do. Only a few minutes ago did you post the first hint of an explanation... You should read this:

http://xyproblem.info/

[Shifo]

I completely agree with you Stephen, I dont have that much experience in matlab but I am learning specially from you guys.