Index exceeds matrix dimensions.

I wrote a following program. Everything is ok to my knowledge. But it shows the error; "Index exceeds matrix dimensions."

clear all

clc

format longEng

syms x y z kh

a=[4;0.5;70];

% The Newton-Raphson iterations starts here

KH=linspace(0,0.55,12);

h=4;

q=100;

gma=14; nq=2*q/(gma*h);

delta=10;

phi=30;

lam=0.1;

A=lam*nq/(1+nq);

kv=0;

% kh=0;

da1=delta*(pi/180); da2=-delta*(pi/180); pha1=phi*(pi/180); pha2=phi*(pi/180);

dp1=delta*(pi/180); dp2=delta*(pi/180); php1=phi*(pi/180); php2=phi*(pi/180);

psi=atan(kh/(1-kv));

a=pha1+da1;

b=pha1-psi;

c=psi+da1;

alphac=atan((sin(a)*sin(b)+(sin(a)^2+sin(b)^2+sin(a)*cos(a)*sin(b)*cos(b)+A*cos(c)*cos(a)*sin(b))^0.5)/(A*cos(c)+sin(a)*cos(b)));

w=0.5*gma*h^2*(1/tan(alphac));

B=h*((1/tan(alphac))-lam);

ka1=((1+nq)*(1-A*tan(alphac))*(cos(b)-(sin(b)/tan(alphac))))/(cos(psi)*(cos(a)+tan(alphac)*sin(a)));

va2=asin(sin(da2)/sin(pha2))-asin(sin(psi)/sin(pha2))-da2-psi;

ka2=(1/cos(psi))*(cos(da2)*((cos(da2)-sqrt(sin(pha2)^2-sin(da2)^2)))/(cos(psi)+sqrt(sin(pha2)^2-sin(psi)^2)))*exp(-va2*tan(pha2));

vp1=asin(sin(dp1)/sin(php1))+asin(-sin(psi)/sin(php1))+dp1+psi;

kp1=(1/cos(psi))*(cos(dp1)*((cos(dp1)+sqrt(sin(php1)^2-sin(dp1)^2)))/(cos(psi)-sqrt(sin(php1)^2-sin(psi)^2)))*exp(vp1*tan(php1));

vp2=asin(sin(dp2)/sin(php2))+asin(-sin(psi)/sin(php2))+dp2+psi;

kp2=(1/cos(psi))*(cos(dp2)*((cos(dp2)+sqrt(sin(php2)^2-sin(dp2)^2)))/(cos(psi)-sqrt(sin(php2)^2-sin(psi)^2)))*exp(vp2*tan(php2));

% Qa1=q*ka1*(h+x); Qp2=q*kp2*y;

sinda1=sin(da1); sindp1=sin(dp1); sinda2=-sin(da2); sindp2=sin(dp2);

cosda1=cos(da1); cosdp1=cos(dp1); cosda2=cos(da2); cosdp2=cos(dp2);

% pa1=ka1*gma*0.5*(h+x)^2;

pp1=kp1*gma*0.5*(x^2);

pa1=0.5*gma*(1-kv)*ka1*h^2;

pa2=ka2*gma*(x*y+0.5*(y^2)); pp2=kp2*gma*(y*(h+x)+(0.5*(y^2)));

zp1=x/3;

za1=((((w*(1-kv)-pa1*sin(da1))*cos(pha1)*h)/(sin(alphac)*cos(alphac-pha1)))-(q*B*(lam*h+0.5*B)))/(pa1*cos(da1));

zp2=((0.5*(h+x)*(y^2))+((y^3)/3))/(((h+x)*y)+(0.5*(y^2)));

za2=((0.5*x*(y^2))+((y^3)/3))/((x*y)+(0.5*(y^2)));

e2=(pp1*cosdp1)+(pa2*cosda2)-(pa1*cosda1)-(pp2*cosdp2);

e3=(pp1*cosdp1*zp1)+(pp2*cosdp2*zp2)-(pa1*cosda1*za1)-(pa2*cosda2*za2);

g=[e2; e3]

J=jacobian([e2, e3], [x, y]);

A=zeros(2,numel(KH));

for i=1:numel(KH)

del=1;

indx=0;

while del>1e-6

gnum = vpa(subs(g,[x,y,kh],[a(1),a(2),KH(i)]));

Jnum = vpa(subs(J,[x,y,kh],[a(1),a(2),KH(i)]));

delx = -Jnum\gnum;

a = a + delx;

del = max(abs(gnum));

indx = indx + 1;

end

A(:,i)=double(a)

end

Answer 1

[KALYAN ACHARJYA]

Initially a defined as vectors

a=[4;0.5;70];  %where a(1)=4, a(2)=0.5 and a(3)=70

After that a is replaced in line 28, as scalars

a=pha1+da1;

But in line 61 & 62 you acess the a data as array as a(1) and a(2), where no a(2) data is avaliable, as a replaced by

a=pha1+da1; 

Thats why it shows error in line 61, next in line 62

gnum=vpa(subs(g,[x,y,kh],[a(1),a(2),KH(i)]));
Jnum=vpa(subs(J,[x,y,kh],[a(1),a(2),KH(i)]));

Command Window

>> a=pha1+da1
a =
    698.131700797732e-003
>> a(1)
ans =
    698.131700797732e-003
>> a(2)
Index exceeds matrix dimensions.

The problem way out is define this a=[4;0.5;70]; and this a=pha1+da1; in different variable name.

Hope you get the issue.

[Akshay Pratap Singh]

Thank you Kalyan. It resolved my issue.

Thanks a lot.