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

4 次查看(过去 30 天)
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,

回答(2 个)

Walter Roberson
Walter Roberson 2019-7-5
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;
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)

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
  6 个评论
Akshay Pratap Singh
Linking the HF and MF to newton raphson code. I am not finding the way to link it. Please suggest.

请先登录,再进行评论。

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by