Problem using solve (Symbolic Toolbox)

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

[madhan ravi]

How would that result in L^3?

Answer 1

[Walter Roberson]

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.

[madhan ravi]

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]

Thank you very much, yes it works now.

Answer 2

[madhan ravi]

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 .