Problem using solve (Symbolic Toolbox)
2 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
Angus Wong
il 10 Lug 2020
Modificato: Angus Wong
il 10 Lug 2020
I encountered a problem in which:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/329795/image.jpeg)
I expect it to return -2*t*L^3/15.
The unfinished code is as follows, with the input of
syms L t positive real;
shearcentre([0 0;0 L;L*sqrt(2) L/2],sym([3 1 3/4*t;1 2 t;2 3 3/4*t]))
on the command line.
Many thanks.
function [ec,final,qvec]=shearcentre(coor,barprop)
% finding shear centre
arguments
coor(:,2,1)sym % [x y]
barprop(:,3,1)sym % [n1 n2 t]
end
close all;
[final,coor,barprop]=beamprop(coor,barprop); % barprop=[n1 n2 t length]
try
double(coor);
double(final);
double(barprop);
catch
syms t L;
assume([t L],{'real','positive'});
coor=simplify(coor);
final=simplify(final);
barprop=simplify(barprop);
end
if isAlways(final(4)~=final(3))
theta=atan(2*final(5)/(final(4)-final(3)))/2+pi;
else
theta=sym(pi);
end
coor=([cos(theta) sin(theta);-sin(theta) cos(theta)]*coor')';
IX=final(3)*cos(theta)^2+final(4)*sin(theta)^2-2*final(5)*sin(theta)*cos(theta);
IY=final(3)*sin(theta)^2+final(4)*cos(theta)^2+2*final(5)*sin(theta)*cos(theta);
occur=sym(zeros(size(coor,1),1));
for ct1=1:size(coor,1)
for ct2=1:size(barprop,1)
if barprop(ct2,1)==ct1 || barprop(ct2,2)==ct1
occur(ct1)=occur(ct1)+1;
end
end
end
if sum(isAlways(occur==0))>0
error('Error! Invalid inputs!');
end
qvec=sym(zeros(size(barprop,1),1));
F=[qvec qvec qvec qvec];
q=sym(zeros(size(coor,1),4));
syms lv positive real;
while size(find(occur==1,1),1)==1
start=find(occur==1,1);
[mem,~]=find(barprop(:,1:2)==start);
bar=barprop(mem,:);
F(start(1),4)=-sign(sum([coor(bar(2),1)-coor(bar(1),1);coor(bar(2),2)-coor(bar(1),2)] ...
.*[-coor(bar(1),2);coor(bar(1),1)]));
if isAlways(start==bar(2))==1
bar(1:2)=fliplr(bar(1:2));
end
% Vy
l1=coor(bar(1),2);
l2=(coor(bar(1),2)+coor(bar(2),2))/2;
lbar(lv)=((l2-l1)/bar(4))*lv+l1;
A(lv)=bar(3)*lv;
F(start,1)=int(q(start,1)+lbar*A,[0 bar(4)])/IX;
q(bar(2),1)=q(start,1)+q(bar(2),1)+lbar(bar(4))*A(bar(4))/IX;
qvec(start,1)=q(start,1)+lbar*A/IX;
figure('NumberTitle','off');
fplot(qvec(start,1),[0 double(bar(4))]);
title(sprintf('Vy, from node %g to node %g',start,bar(2)));
grid minor;
% Vx
l1=coor(bar(1),1);
l2=(coor(bar(1),1)+coor(bar(2),1))/2;
lbar(lv)=((l2-l1)/bar(4))*lv+l1;
F(start,2)=int(q(start,3)+lbar*A,[0 bar(4)])/IY;
q(bar(2),3)=q(start,3)+q(bar(2),3)+lbar(bar(4))*A(bar(4))/IY;
qvec(start,2)=q(start,3)+lbar*A/IY;
figure('NumberTitle','off');
fplot(qvec(start,2),[0 double(bar(4))]);
title(sprintf('Vx, from node %g to node %g',start,bar(2)));
grid minor;
% Lever Arm (F(:,3))
m=(coor(bar(2),2)-coor(bar(1),2))/(coor(bar(2),1)-coor(bar(1),1));
if isAlways(m==0)==1
F(start,3)=abs(coor(bar(1),2));
elseif isAlways(isinf(m))==1 || isAlways(isinf(-m))==1
F(start,3)=abs(coor(bar(1),1));
else
F(start,3)=abs((coor(bar(1),2)-m*coor(bar(1),2))/(m+1/m)*sqrt(1+m^-2));
end
% finishing step
barprop=barprop([1:(mem-1) (mem+1):end],:);
occur(start)=occur(start)-1;
occur(bar(2))=occur(bar(2))-1;
end
if sum(occur)~=0
% box section
switch size(find(occur==3,1),1)
case 0
syms qv positive real;
qt=sym(zeros(size(barprop,1),5));
node=find(occur~=0);
qt(:,1)=node';
start=min(node);
qt(start,2)=-qv;
qt(start,3)=-qv;
[mem,~]=find(barprop(:,1:2)==start,1);
barpropt=barprop;
qti=sym([0;0]);
while sum(occur)>0 % assume a cut in the min node
bar=barpropt(mem,:);
if isAlways(start==bar(2))==1
bar(1:2)=fliplr(bar(1:2));
end
% Vy
l1=coor(bar(1),2);
l2=(coor(bar(1),2)+coor(bar(2),2))/2;
lbar(lv)=((l2-l1)/bar(4))*lv+l1;
A(lv)=bar(3)*lv;
qt(bar(2),2)=qt(start,2)+qt(bar(2),2)*(size(barpropt,1)>1)+lbar(bar(4))*A(bar(4));
qti(1)=qti(1)+int((qt(start,2)+lbar*A)/bar(3),lv,[0 bar(4)]);
% Vx
l1=coor(bar(1),1);
l2=(coor(bar(1),1)+coor(bar(2),1))/2;
lbar(lv)=((l2-l1)/bar(4))*lv+l1;
qt(bar(2),3)=qt(start,3)+qt(bar(2),3)*(size(barpropt,1)>1)+lbar(bar(4))*A(bar(4));
qti(2)=qti(2)+int((qt(start,3)+lbar*A)/bar(3),lv,[0 bar(4)]);
% finishing step
occur(start)=occur(start)-1;
occur(bar(2))=occur(bar(2))-1;
barpropt(mem,:)=sym(zeros(1,4));
start=bar(2);
if size(barpropt,1)>=1
[mem,~]=find(barpropt(:,1:2)==start,1);
end
end
qti(1)=solve(qti(1)==0,qv); % problem encountered here.
qti(2)=solve(qti(2)==0,qv);
case 2
syms q [2,1] positive real;
% not finished
otherwise
error('Error! Too complicated!');
end
end
ex=sum(F(:,1).*F(:,3).*F(:,4))/IX;
ey=sum(F(:,2).*F(:,3).*F(:,4))/IY;
ec=vpa(([cos(theta) -sin(theta);sin(theta) cos(theta)]*[ex;ey]));
end
function [final,coor,barprop]=beamprop(coor,barprop)
arguments
coor(:,2,1)sym % [x y]
barprop(:,3,1)sym % [n1 n2 t]
end
prop=sym(zeros(size(barprop,1),8));
for ct=1:size(barprop,1)
prop(ct,1)=sqrt((coor(barprop(ct,2),1)-coor(barprop(ct,1),1))^2+(coor(barprop(ct,2),2)-coor(barprop(ct,1),2))^2); % length
prop(ct,2)=prop(ct,1)*barprop(ct,3); % Area
prop(ct,3)=(coor(barprop(ct,1),1)+coor(barprop(ct,2),1))/2; % centre x-coor
prop(ct,4)=(coor(barprop(ct,1),2)+coor(barprop(ct,2),2))/2; % centre y-coor
if coor(barprop(ct,2),1)==coor(barprop(ct,1),1)
prop(ct,5)=pi/2; % angle of bar
else
prop(ct,5)=atan((coor(barprop(ct,2),2)-coor(barprop(ct,1),2))/(coor(barprop(ct,2),1)-coor(barprop(ct,1),1))); % angle of bar
end
end
barprop=[barprop prop(:,1)];
final=sym(zeros(5,1));
final(1)=sum(prop(:,2).*prop(:,3))/sum(prop(:,2));
final(2)=sum(prop(:,2).*prop(:,4))/sum(prop(:,2));
prop(:,3)=prop(:,3)-final(1);
prop(:,4)=prop(:,4)-final(2);
coor(:,1)=coor(:,1)-final(1);
coor(:,2)=coor(:,2)-final(2);
for ct=1:size(barprop,1)
X=prop(ct,3)*cos(prop(ct,5))+prop(ct,4)*sin(prop(ct,5));
Y=-prop(ct,3)*sin(prop(ct,5))+prop(ct,4)*cos(prop(ct,5));
IX=prop(ct,1)*barprop(ct,3)*(barprop(ct,3)^2/12+Y^2);
IY=prop(ct,1)*barprop(ct,3)*(prop(ct,1)^2/12+X^2);
IXY=prop(ct,1)*barprop(ct,3)*X*Y;
prop(ct,6)=IX*cos(-prop(ct,5))^2+IY*sin(-prop(ct,5))^2-2*IXY*sin(-prop(ct,5))*cos(-prop(ct,5)); % Ix
prop(ct,7)=IX*sin(-prop(ct,5))^2+IY*cos(-prop(ct,5))^2+2*IXY*sin(-prop(ct,5))*cos(-prop(ct,5)); % Iy
prop(ct,8)=(IX-IY)*sin(-prop(ct,5))*cos(-prop(ct,5))+IXY*(cos(-prop(ct,5))^2-sin(-prop(ct,5))^2); % Ixy
end
final(3)=simplify(sum(prop(:,6)));
final(4)=simplify(sum(prop(:,7)));
final(5)=simplify(sum(prop(:,8)));
% for ct=1:size(barprop,1)
% barprop(ct,6)=(coor(barprop(ct,1),1)+coor(barprop(ct,2),1))/2;
% barprop(ct,7)=(coor(barprop(ct,1),2)+coor(barprop(ct,2),2))/2;
% end
end
Risposta accettata
Walter Roberson
il 10 Lug 2020
syms qv positive real;
and
syms L t positive real;
and when you solve for qv, you say you are expecting the answer
-2*t*L^3/15
but t and L are positive, so 2*t*L^3/15 would be positive, and negative of that would be negative. So if L and t are positive, qv would have to be negative. But you told it that qv is positive. So there is no solution for qv.
The -(2*L^2*t)/15 that madhan came up with has the same problem, that it would be negative.
2 Commenti
madhan ravi
il 10 Lug 2020
Modificato: madhan ravi
il 10 Lug 2020
Ah it makes sense sir Walter. But the OP edited/posted the relevant part of the code after I posted the answer. Also shouldn’t L^3 be L^2?
Più risposte (1)
madhan ravi
il 10 Lug 2020
Modificato: madhan ravi
il 10 Lug 2020
By the way posting a screenshot is useless , always make sure to upload the code as text so that others could run it.
>> syms L qv t
eqn = formula(-(2*L^3)/3 - (5*L*qv)/t )
solve(eqn(1) == 0, qv)
eqn =
- (2*L^3)/3 - (5*L*qv)/t
ans =
-(2*L^2*t)/15
I get the results in MATLAB mobile 2020a version . The only thing which differs from mine and yours is you’re running it in debug mode, I don’t know if that would be relevant here. I suspect you created a custom solve.m if so please remove it from the path or rename it .
0 Commenti
Vedere anche
Categorie
Scopri di più su Number Theory in Help Center e File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!