清空筛选条件
清空筛选条件

Problem using solve (Symbolic Toolbox)

2 次查看（过去 30 天）
Angus Wong
Angus Wong 2020-7-10
编辑： Angus Wong 2020-7-10
I encountered a problem in which:
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

请先登录，再回答此问题。

采纳的回答

Walter Roberson
Walter Roberson 2020-7-10
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 个评论
madhan ravi
madhan ravi 2020-7-10
编辑：madhan ravi 2020-7-10
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?
Angus Wong
Angus Wong 2020-7-10
Thank you very much, yes it works now.

请先登录，再进行评论。

更多回答（1 个）

madhan ravi
madhan ravi 2020-7-10
编辑：madhan ravi 2020-7-10
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 .

类别

在 Help Center 和 File Exchange 中查找有关 Number Theory 的更多信息

标签

产品


版本

R2020a

Community Treasure Hunt

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

Start Hunting!

Translated by