Error using symengine. The dimensions of matrices or vectors are incompatible.

Why Matlab cannot perform such a calculation ? This is not the first time I perform these kind of calculations.

dbstop if error
clear all
clc
format longEng
syms x y h
phi=(pi/180)*39;
delta=(2*phi)/3;
gma=18.4;
a=[1.51;0.632];
% kh=0.3;
H=linspace(0.5,4,8);
% h=4;
lam=0;
q=50;
nq=2*q/(gma*(h+x));
A=lam*nq/(1+nq);
kh=0;
kv=0;
psi=atan(kh/(1-kv));
beta=1.2;
alfa=1.2;
R1=-1; 
R3=-1;
delm1=0.5*(1-R1)*delta;
delm3=-0.5*(1-R3)*delta;
m=phi+delm1;
b=phi-psi;
c=psi+delm1;

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

kg=(tan(alphac-phi)+(kh/(1-kv)))/(tan(alphac)*(cos(delm1)+sin(delm1)*tan(alphac-phi)));
r=1-lam*tan(alphac);
kq=r*kg;

pg=0.5*gma*(1-kv)*kg*(h+x)^2*cos(delm1);
pq=(1-kv)*q*kq*(h+x);

k3=(2*cos(phi-psi)^2)/(cos(phi-psi)^2*(1+R3)+cos(psi)*cos(delm3+psi)*(1-R3)*(1+sqrt((sin(phi+delm3)*sin(phi-psi))/cos(delm3+psi)))^2);


y1=0:0.01:1;
y2=0:0.01:1;
mmm=size(y1);
mm=mmm(1,2);
for i=1:mm
    if (y1(i)>=1-(1/beta)&& y1(i)<=1)
        R2(i)=3*(beta*(1-y1(i))).^0.5;
    else
        R2(i)=3;
    end
    if (R2(i)>=0 && R2(i)<1)
        delm2(i)=0.5*(1-R2(i)).*delta;
        k2(i)=(2*cos(phi-psi)^2)/(cos(phi-psi)^2*(1+R2(i))+cos(psi)*cos(delm2(i)+psi)...
            *(1-R2(i))*(1+sqrt((sin(phi+delm2(i))*sin(phi-psi))/cos(delm2(i)+psi)))^2);
    else 
        delm2(i)=0.5*(R2(i)-1)*delta;
        k2(i)=1+0.5*(R2(i)-1).*((cos(phi-psi).^2./(cos(psi).*(cos(delm2(i)+psi).*...
            (-sqrt((sin(phi+delm2(i)).*sin(phi-psi))./(cos(delm2(i)+psi)))+1).^2)))-1);
    end
    if (y2(i)>=0 && y2(i)<=(1/alfa))
        R4(i)=3*(alfa*(y2(i)))^0.5;
    else
        R4(i)=3;
    end
     if (R4(i)>=0 && R4(i)<=1)
         delm4(i)=0.5*(1-R4(i)).*delta;
         k4(i)=(2*cos(phi-psi)^2)/(cos(phi-psi)^2*(1+R4(i))+cos(psi)*cos(delm4(i)+psi)...
             *(1-R4(i))*(1+sqrt((sin(phi+delm4(i))*sin(phi-psi))/cos(delm4(i)+psi)))^2);
     else
         delm4(i)=0.5*(R4(i)-1)*delta;
         k4(i)=1+0.5*(R4(i)-1)*((cos(phi-psi)^2/(cos(psi)*(cos(delm4(i)+psi)*(-sqrt((sin(phi+delm4(i))...
             *sin(phi-psi))/(cos(delm4(i)+psi)))+1)^2)))-1);
     end
     
h2(i)=gma.*x.^2.*k2(i).*y1(i).*cos(delm2(i));
h4(i)=gma.*y.^2.*k4(i).*cos(delm4(i)).*y2(i)+gma.*(x+h).*y.*k4(i).*cos(delm4(i));
H3_a(i)=k3.*y2(i).*cos(delm3);
% H3_b=matlabFunction(k3*cos(delm3));
h3=gma.*y.^2.*H3_a+gma.*x.*y.*k3.*cos(delm3);%integral(H3_b,0,1);
HF=h2-h4+h3-pg-pq;
%
M2(i)=k2(i).*y1(i).*cos(delm2(i)).*(1-y1(i));
m2=gma.*x.^3.*M2;
M4_a(i)=y2(i).^2;
M4_b(i)=y2(i);
m4=gma.*y.^3.*k4(i).*cos(delm4(i)).*M4_a+gma.*(x+h).*y.^2.*k4(i).*cos(delm4(i)).*M4_b;
M3_a(i)=k3.*y2(i).^2.*cos(delm3);
M3_b(i)=k3.*y2(i).*cos(delm3);
m3=gma.*y.^3.*M3_a+gma.*x.*y.^2.*M3_b;
MF=m2+m4-m3-pg.*(h+x).*(1/3)-0.5.*pq.*(h+x);
end
% The Newton-Raphson iterations starts here
g=[HF; MF];
J=jacobian([HF, MF], [x, y]);
Z=zeros(2,numel(H));

for j=1:numel(H)
    del=1;
    indx=0;
while del>1e-6
   gnum = vpa(subs(g,[x,y,h],[a(1),a(2),H(j)]));
   Jnum = vpa(subs(J,[x,y,h],[a(1),a(2),H(j)]));
   delx = -Jnum\gnum;
   a = a + delx;
   del = max(abs(gnum));
   indx = indx + 1;
end
 Z(:,j)=double(a)
end
'NEWTON-RAPHSON SOLUTION CONVERGES IN ITERATIONS',indx,
'FINAL VALUES OF a ARE';a,

Answer 1

[Walter Roberson]

g=[HF; MF];
J=jacobian([HF, MF], [x, y]);

In the first of those, you arrange HF and MF into rows, getting a 2 x N result for N = size(HF,2)

In the second of those, you arrange HF and MF into columns, getting a 1 x (2*N) result, which you then take the jacobian of with respect to two variables, getting a (2*N) x 2 result.

   gnum = vpa(subs(g,[x,y,h],[a(1),a(2),H(j)]));
   Jnum = vpa(subs(J,[x,y,h],[a(1),a(2),H(j)]));

The vpa(subs()) does not change the shape when given those arguments, so gnum is going to be 2 x N, and Jnum is going to be (2*N) x 2.

      delx = -Jnum\gnum;

https://www.mathworks.com/help/matlab/ref/mldivide.html

If A is a rectangular m-by-n matrix with m ~= n, and B is a matrix with m rows, then A\B returns a least-squares solution to the system of equations A*x= B.

In this notation, A = Jnum, m = (2*N) and n = 2, and B = gnum which must be a matrix with (2*N) rows -- but it has 2 rows instead.

In context it appears to me that you want 2 outputs. You need to decide whether you want that as a 1 x 2 output or as a 2 x 1 output, in order to figure out how to arrange your input variables. The number of columns of output is equal to the number of columns of B (here, gnum)

Answer 2

[Akshay Pratap Singh]

Thank you Roberson for suggestion:

May you please help me out in solving these two nonlinear equations, HF and MF. The HF and MF have been calculated for various k1 and k2. The given newton raphson code is suitable for constant k2 and k4 values and here k2 and k4 are varying.

Thank you in advance

[Walter Roberson]

Is it okay if the solutions are complex valued?

[Akshay Pratap Singh]

No, the solution will not be complex valued.

[Walter Roberson]

In order to handle varying k2 and k4, the easiest way is to loop using each combination in turn. The approach you are taking of trying to solve for all of the values simultaneously will not work: it will complain that the equations are inconsistent.

[Akshay Pratap Singh]

True, the k2 and k4 are the values of a segment that I have to make a function. This thing I am not able to do.

[Walter Roberson]

What do you find is preventing you from doing that?

[Akshay Pratap Singh]

Linking the HF and MF to newton raphson code. I am not finding the way to link it. Please suggest.