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,
0 个评论
回答(2 个)
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)
0 个评论
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!