# Advanced Mathematics and Mechanics Applications Using MATLAB, 3rd Edition

### Howard Wilson (view profile)

14 Oct 2002 (Updated )

Companion Software (amamhlib)

[vn,rcn,irrn,vs,rcs,irrs,times,nt,np]=...
```function [vn,rcn,irrn,vs,rcs,irrs,times,nt,np]=...
ropesymu(A,B,M,X0,Y0,Z0,nt,np)
%
% [vn,rcn,irrn,vs,rcs,irrs,times,nt,np]=ropesymu(...
%                                A,B,M,X0,Y0,Z0,nt,np)
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%
% This program computes geometrical properties of a
% twisted rope having a cross section which is two
% circles of diameter B touching tangentially. The
% tangency point is at distance A from the rotation
% axis z. As the area is rotated, it is also twisted
% in a helical fashion. For a complete revolution
% about the z axis, the area is twisted through m
% turns. The resulting surface resembles a rope
% composed of two strands. Two results are obtained
% 1) by a numerical method where the surface is
% modeled with triangular surface patches and
% 2) by symbolic math. See functions twistrope and
% twistprop for descriptions of the problem parameters.
% Numerical results and computation times for the two
% methods are compared, and the related surface
% geometry is plotted
%
% User functions called: twistrope twistprop ropedraw
%----------------------------------------------------

if nargin==0 % Default data case
A=3; B=1; m=6; np=201; nt=25;
X0=0; Y0=0; Z0=-3*pi/2; M=6;
end

disp(' ')
disp('   COMPARISON OF NUMERICAL AND SYMBOLIC')
disp('GEOMETRICAL PROPERTIES FOR A TWISTED ROPE')

% Run the first time to get a crude grid for plotting
[vn,rcn,irrn,x,y,z,c]=twistrope(A,B,M,X0,Y0,Z0,nt,np);

% Numerical solution using a dense point grid to get
% close comparison with exact results. Calculations
% are run repeatedly for accurate timing.
Nt=4*nt; Np=4*np; n=50; tic;
for i=1:n
[vn,rcn,irrn]=twistrope(A,B,M,X0,Y0,Z0,Nt,Np);
end
timn=toc/n;

% Perform the symbolic analysis. This takes a long
% time.
tic;
[v,rc,vrr,vs,rcs,irrs]=twistprop(A,B,M,X0,Y0,Z0);
tims=toc; times=[timn,tims];

disp(' ')
disp('FOR THE TRIANGULAR SURFACE PATCH MODEL')
disp(['Volume = ',num2str(vn)])
disp(['Rg = [',num2str(rcn(:)'),']'])
disp('Irr = '), disp(irrn)
disp(['Computation Time = ',num2str(timn),' Secs.'])

% Print numerical comparisons of results
disp(' ')
disp('FOR THE SYMBOLIC MODEL')
disp(['Volume = ',num2str(vs)])
disp(['Rg = [',num2str(rcs(:)'),']'])
disp('Irr = '), disp(irrs)
disp(['Computation Time = ',num2str(tims),' Secs.'])

disp(' ')
disp(' NUMERICAL APPROXIMATION ERROR USING TRIANGULAR')
disp('SURFACE PATCHES. THE ERROR VALUES ARE DEFINED AS')
disp('        NORM(APPROX.-EXACT)/NORM(EXACT)')
einert=norm(irrn-irrs)/norm(irrs);
disp(['Volume Error = ',num2str(evol)])
disp(['Inertia Tensor Error = ',num2str(einert)])

disp(' ')
disp('COMPARISON OF SOLUTION TIMES')
disp(['(Symbolic Time)/(Numerical Time) = ',...
num2str(tims/timn)])
disp(' ')

% Draw the surface using a crude grid to avoid
% crowded grid lines
ropedraw(A,B,np,nt,M,X0,Y0,Z0);

%===========================================

function [x,y,z,t]=ropedraw(a,b,np,nt,m,x0,y0,z0)
%
% [x,y,z,t]=ropedraw(a,b,np,mp,m,x0,y0,z0)
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% This function draws the twisted rope.
if nargin==0
a=3; b=1; np=200; nt=25; m=6;
x0=0; y0=0; z0=-3*pi/2;
end

% Draw the surface
t=linspace(0,2*pi,nt); p=linspace(0,3*pi,np)';
t=repmat(t,np,1); p=repmat(p,1,nt);
xi=b*cos(t).*abs(cos(t)); eta=b*sin(t).*abs(cos(t));
rho=a+xi.*cos(m*p)+eta.*sin(m*p);
x=rho.*cos(p)+x0; y=rho.*sin(p)+y0;
z=-xi.*sin(m*p)+eta.*cos(m*p)+p+z0;
close; surf(x,y,z,t), title('TWISTED ROPE')
xlabel('x axis'), ylabel('y axis'), zlabel('z axis')
colormap('prism(4)'), axis equal, hold on

% Fill the ends
fill3(x(1,:),y(1,:),z(1,:),'w')
fill3(x(end,:),y(end,:),z(end,:),'w')
view([-40,10]), hold off, shg

%===========================================

function [v,rc,vrr,V,Rc,Irr]=twistprop(A,B,M,X0,Y0,Z0)
%
% [v,rc,vrr,V,Rc,Irr]=twistprop(A,B,M,X0,Y0,Z0)
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% This function computes geometrical properties of
% a twisted rope. Exact results are obtained using
% symbolic math to evaluate three surface integrals
% for the volume, centroidal radius, and inertia
% tensor. The symbolic calculations take about five
% minutes to run.
%
% A,B,N    - parameters defining the twisted rope
% X0,Y0,Z0 - center coordinates for the centroid of
%            the twisted rope
% v,rc     - symbolic formulas for the volume and
% vrr      - symbolic formula for integral of
%            r*r'*d(vol)
% V,Rc     - numerical values for volume and
% Irr      - numerical value for the inertia tensor

if nargin==0
A=6; B=1; M=6; X0=1; Y0=2; Z0=3;
end

syms a b m t p xi eta rho x y z r rt rp real
syms x0 y0 z0 real
syms n dv dv1 v vr1 vr rg vrr1 vrr real
a=sym(A); b=sym(B); Pi=sym('pi');
x0=sym(X0); y0=sym(Y0); z0=sym(Z0);

% Surface equation for the twisted rope
xi=b*cos(t)*abs(cos(t));
eta=b*sin(t)*abs(cos(t));
rho=a+xi*cos(m*p)+eta*sin(m*p);
x=rho*cos(p)+x0; y=rho*sin(p)+y0;
z=-xi*sin(m*p)+eta*cos(m*p)+p+z0;
Pi=sym('pi');

% Tangent vectors
r=[x;y;z]; rt=diff(r,t); rp=diff(r,p);

% Integrate to get the volume
dv=det([r,rp,rt]); dv1=int(dv,t,0,2*Pi);
v=simple(int(dv1,p,0,3*Pi)/3);

% First moment of volume
vr1=int(r*dv,t,0,2*Pi);
vr=simple(int(vr1,p,0,3*Pi)/4);

% Radius to the centroid
rc=simple(vr/v);

% Integral of r*r'*d(vol)
vrr1=int(r*r'*dv,t,0,2*Pi);
vrr=simple(int(vrr1,p,0,3*Pi)/5);

% Obtain numerical values
V=double(subs(v,{a,b,m,x0,y0,z0},...
{A,B,M,X0,Y0,Z0}));
Rc=double(subs(rc,{a,b,m,x0,y0,z0},...
{A,B,M,X0,Y0,Z0}));
Irr=double(subs(vrr,{a,b,m,x0,y0,z0},...
{A,B,M,X0,Y0,Z0}));

% Rigid body inertia tensor for a
% body of unit mass density
Irr=eye(3,3)*sum(diag(Irr))-Irr;

%===========================================

function [v,rc,vrr,x,y,z,t]=twistrope(...
a,b,m,x0,y0,z0,nt,np)
%
% [v,rc,vrr,x,y,z,t]=twistrope(...
%                  a,b,m,x0,y0,z0,nt,nm)
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Geometrical properties of a twisted rope.
% This example takes 1.3 seconds to run
if nargin<8, np=321; end; if nargin<7, nt=161; end
if nargin==0
a=6; b=1; m=6; x0=1; y0=2; z0=3;
end
t=linspace(0,2*pi,nt); p=linspace(0,3*pi,np)';
t=repmat(t,np,1); p=repmat(p,1,nt);

% Surface equation for the twisted rope

xi=b*cos(t).*abs(cos(t));
eta=b*sin(t).*abs(cos(t));
rho=a+xi.*cos(m*p)+eta.*sin(m*p);
x=rho.*cos(p)+x0; y=rho.*sin(p)+y0;
z=-xi.*sin(m*p)+eta.*cos(m*p)+p+z0;

[v,rc,vrr]=srfv(x,y,z);

%===========================================

function [v,rc,vrr]=srfv(x,y,z)
%
% [v,rc,vrr]=srfv(x,y,z)
% ~~~~~~~~~~~~~~~~~~~~~~
%
% This function computes the volume, centroidal
% coordinates, and inertial tensor for a volume
% covered by surface coordinates contained in
% arrays x,y,z
%
% x,y,z   - matrices containing the coordinates
%           of a grid of points covering the
%           surface of the solid
% v       - volume of the solid
% rc      - centroidal coordinate vector of the
%           solid
% vrr     - inertial tensor for the solid with the
%           mass density taken as unity
%
% User functions called: scatripl proptet
%-----------------------------------------------

% p=inline(...
%  'v*(eye(3)*(r(:)''*r(:))-r(:)*r(:)'')','v','r');

%d=mean([x(:),y(:),z(:)]);
%x=x-d(1); y=y-d(2); z=z-d(3);

[n,m]=size(x); i=1:n-1; I=i+1; j=1:m-1; J=j+1;
xij=x(i,j); yij=y(i,j); zij=z(i,j);
xIj=x(I,j); yIj=y(I,j); zIj=z(I,j);
xIJ=x(I,J); yIJ=y(I,J); zIJ=z(I,J);
xiJ=x(i,J); yiJ=y(i,J); ziJ=z(i,J);

% Tetrahedron volumes
v1=scatripl(xij,yij,zij,xIj,yIj,zIj,xIJ,yIJ,zIJ);
v2=scatripl(xij,yij,zij,xIJ,yIJ,zIJ,xiJ,yiJ,ziJ);
v=sum(sum(v1+v2));

% First moments of volume
X1=xij+xIj+xIJ; X2=xij+xIJ+xiJ;
Y1=yij+yIj+yIJ; Y2=yij+yIJ+yiJ;
Z1=zij+zIj+zIJ; Z2=zij+zIJ+ziJ;
vx=sum(sum(v1.*X1+v2.*X2));
vy=sum(sum(v1.*Y1+v2.*Y2));
vz=sum(sum(v1.*Z1+v2.*Z2));

% Second moments of volume
vrr=proptet(v1,xij,yij,zij,xIj,yIj,zIj,...
xIJ,yIJ,zIJ,X1,Y1,Z1)+...
proptet(v2,xij,yij,zij,xIJ,yIJ,zIJ,...
xiJ,yiJ,ziJ,X2,Y2,Z2);
rc=[vx,vy,vz]/v/4; vs=sign(v);
v=abs(v)/6; vrr=vs*vrr/120;
vrr=[vrr([1 4 5]), vrr([4 2 6]), vrr([5 6 3])]';
vrr=eye(3,3)*sum(diag(vrr))-vrr;

%vrr=vrr-p(v,rc)+p(v,rc+d); rc=rc+d;

%===========================================

function v=scatripl(ax,ay,az,bx,by,bz,cx,cy,cz)
%
% v=scatripl(ax,ay,az,bx,by,bz,cx,cy,cz)
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Scalar triple product dot(cross(a,b),c) where
% the cartesian components of vectors a,b,and c
% are given in arrays of the same size.
v=ax.*(by.*cz-bz.*cy)+ay.*(bz.*cx-bx.*cz)...
+az.*(bx.*cy-by.*cx);

% =========================================

function vrr=tensprod(v,x,y,z)
%
% vrr=tensprod(v,x,y,z)
% ~~~~~~~~~~~~~~~~~~~~
% This function forms the various components
% of v*R*R'. The calculation is vectorized
% over arrays of points
vxx=sum(sum(v.*x.*x)); vyy=sum(sum(v.*y.*y));
vzz=sum(sum(v.*z.*z)); vxy=sum(sum(v.*x.*y));
vxz=sum(sum(v.*x.*z)); vyz=sum(sum(v.*y.*z));
vrr=[vxx; vyy; vzz; vxy; vxz; vyz];

% =========================================

function vrr=proptet(v,x1,y1,z1,x2,y2,z2,...
x3,y3,z3,xc,yc,zc)
%
% vrr=proptet(v,x1,y1,z1,x2,y2,z2,x3,y3,z3,...
%                                    xc,yc,zc)
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% This function computes tensor properties of a
% tetrahedron with its base being a triangular
% surface and its apex at the origin
vrr=tensprod(v,x1,y1,z1)+tensprod(v,x2,y2,z2)+...
tensprod(v,x3,y3,z3)+tensprod(v,xc,yc,zc);```