Project Help Notebook by Keith Stroyan

Help NoteBook Copyright 1998 by Academic Press, Inc. - All rights reserved.
Modify this Notebook to prepare your project report.

Project 13: The Catenary

Initialization

This horrible stuff sets various defaults and defines some functions needed to draw the geometric side of the vectro lexicon.

Off[General::spell];
Off[General::spell1];
Off[ParametricPlot3D::ppcom]
Needs["Graphics`Colors`"]

SetOptions[Graphics,Axes->False,AxesOrigin->{0,0},AspectRatio->Automatic];
SetOptions[Graphics3D,ViewPoint->{5,2,1},Boxed->False,AspectRatio->Automatic];
SetOptions[ParametricPlot3D,ViewPoint->{5,2,1},AspectRatio->Automatic];
SetOptions[Plot3D,ViewPoint->{5,2,1},AspectRatio->Automatic];
SetOptions[ParametricPlot,AspectRatio->Automatic];
SetOptions[Plot,AspectRatio->Automatic];

colorQ[a_]:=(Head[a]===RGBColor);
grayQ[a_]:=(Head[a]===GrayLevel);
thickQ[a_]:=(Head[a]===Thickness);
absthickQ[a_]:=(Head[a]===AbsoluteThickness);
notruleQ[u_] := !(Head[u]===Rule);
ruleQ[u_] := (Head[u]===Rule);

Norm::usage = "Norm[expr] returns the Euclidean length of the
		vector, 'expr'.";

Norm[a_?VectorQ] := Sqrt[a.a];

Angle[a_?VectorQ,b_?VectorQ] := ArcCos[(a.b)/(Norm[a] Norm[b])]

Proj[a_?VectorQ,b_?VectorQ] := ((a.b)/(b.b))b

tol = .5*10^-6;

OrthUp[a_?VectorQ] :=
	Module[{b,c,d,f},
		d = N[a/Norm[a]];
		If[Chop[d[[3]]^2 - 1,tol] == 0, b={1,0,0} ,
		b = {0,0,1} - d[[3]] d; b = b/Norm[b];];
	Return[b]];
	
OrthSet[a_?VectorQ] :=
	Module[{b,c,d},
		d = N[a/Norm[a]];
		If[Chop[d[[1]],tol]==0,b={1,0,0},
			 If[Chop[d[[2]],tol]==0,b={0,1,0},
			  b={-d[[2]],d[[1]],0}]];
			b=b/Norm[b];c = Cross[d,b];
	Return[{b,c/Norm[c],d}]];

				
Rotate[a_?VectorQ,theta_] :=
	Module[{u,v,r},
		u = OrthSet[a];
		v = Transpose[u];
		r = N[{{Cos[theta],-Sin[theta],0},{Sin[theta],Cos[theta],0},{0,0,1}}];
	Return[v.r.u]];

Vector2D[tail_,tip_,opts___]:=
	Module[{ntail,ntip,sopts,Punit,up,arrow},
	ntail = N[tail];
	ntip = N[tip];
	
	sopts = Join[Select[{opts},colorQ],Select[{opts},
				thickQ],Select[{opts},absthickQ],Select[{opts},grayQ]];

If[Chop[ntail.ntail+ntip.ntip,tol] == 0,
			(*too small to plot or not number*)
			Return[Graphics[Point[{0,0}]]];
	,						
If[Chop[(ntip-ntail).(ntip-ntail),tol] == 0,
		Return[Graphics[{Join[sopts,{Line[{ntail,ntip}]}]}]];
	,
	Punit = 0.075 (ntip-ntail) ;
	up = 0.15 {-Punit[[2]],Punit[[1]]};
	arrow = Line[{ntail,ntip,ntip-Punit+up,ntip,ntip-Punit-up}];
	]];
Return[Graphics[{Join[{AbsoluteThickness[2]},Join[sopts,{arrow}]]}]];
];

Vector3D[tail_,tip_,opts___]:=
	Module[{ntail,ntip,sopts,Punit,up,u2,u3,u4,rotmat,arrow},
	ntail = N[tail];
	ntip = N[tip];
	
	sopts = Join[Select[{opts},colorQ],Select[{opts},
				thickQ],Select[{opts},absthickQ],Select[{opts},grayQ]];

If[Chop[ntail.ntail+ntip.ntip,tol] == 0,
			(*too small to plot or not number*)
			Return[Graphics3D[Point[{0,0,0}]]];					
	,		
If[Chop[(ntip-ntail).(ntip-ntail),tol] == 0,
		Return[Graphics3D[{Join[sopts,{Line[{ntail,ntip}]}]}]];
	,
	Punit = (ntip-ntail);
	rotmat = Rotate[Punit,N[Pi/2]];
	up = OrthUp[Punit];
	Punit = .075 Punit;
	up = .15 Norm[Punit] up;
	u2 = rotmat.up;
	u3 = rotmat.u2;
	u4 = rotmat.u3;
	arrow = Line[{ntail,ntip,
				ntip - Punit + up,ntip,
				ntip - Punit + u2,ntip,
				ntip - Punit + u3,ntip,
				ntip - Punit + u4}];
	]];
Return[Graphics3D[{Join[{AbsoluteThickness[2]},Join[sopts,{arrow}]]}]];
];

aRCtAN[y_,x_] := Module[{v},
	If[ Chop[N[x],tol] == 0 ,If[ N[y] < 0, v = -N[Pi/2], v = N[Pi/2]],
		v = N[ArcTan[y/x]];
		If[ x < 0 && y < 0, v = N[Pi -  v],v]];
	Return[v]];

Arc2D[X_,Y_,opts___] := Graphics[{opts,
		Disk[{0,0},1,{Min[
							aRCtAN[
						N[X[[2]]],N[X[[1]]]
								]
							,
							aRCtAN[
							N[Y[[2]]],N[Y[[1]]]
								]
						]
						,
					Max[
							aRCtAN[
						N[X[[2]]],N[X[[1]]]
								]
							,
							aRCtAN[
							N[Y[[2]]],N[Y[[1]]]
								]
						]}
				]}
				];
Arc3D[X_,Y_] := 
	Module[{t,U,V},
		U[V_]:= N[(1/Sqrt[V.V]) V];
		Return[Graphics3D[
		Polygon[Append[Table[U[t X + (1-t) Y],{t,0,1,0.1}],{0,0,0}]]
				]
		]
	];

(*Zero Arrow causes trouble. Put in Vector?*)

Vector[tail_,tip_,opts___]:=
	If[Length[N[tail]]==2,Return[Vector2D[tail,tip,opts]],
		Return[Vector3D[tail,tip,opts]]];
		
Arrow[vec_,opts___]:= 
	Module[{a,b},
				If[Length[N[vec]]==2,
				a=ToExpression["Tail"]/.Select[{opts},ruleQ]/.
					{ToExpression["Tail"]->{0,0}};
				Return[Vector2D[a,a+vec,opts]],
				a=ToExpression["Tail"]/.Select[{opts},ruleQ]/.
					{ToExpression["Tail"]->{0,0,0}};
				Return[Vector3D[a,a+vec,opts]]]
				];

Frame2D[delta_] := 
      Graphics[{
            {Green,Line[{{-delta/4,0},{delta,0}}]},
                    Text["x",{1.01 delta,0}],
            {Green,Line[{{0,-delta/4},{0,delta}}]},
                    Text["y",{0, 1.01 delta}]}];

Frame3D[delta_] := 
      Graphics3D[{
            {Green,Line[{{-delta/4,0,0},{delta,0,0}}]},
                    Text["x",{1.01 delta,0,0}],
            {Green,Line[{{0,-delta/4,0},{0,delta,0}}]},
                    Text["y",{0, 1.01 delta,0}],
            {Green,Line[{{0,0,-delta/4},{0,0,delta}}]},
                    Text["z",{0,0,1.01 delta}]}];

View2D[stuff_] := Show[stuff,Frame2D[1]];
View3D[stuff_] := Show[stuff,Frame3D[1]];

by A. Student and B. Pupil

with technical help from Da LabMonitor and Help Notebook

Hanging the Clothesline

Write this last.  Begin with the basic question and state your main finding clearly, but with few technical details.  You might even modify the title of this section to make it more interesting.

Beginning Explorations or Computer Experiments

The derivation of the equation of the catenary is a paper and pencil project, but following are some drawings to help in your write-up.

A General Diagram of Forces

The following program draws a segment of a catenary and shows the tension and weight forces acting on the segment.  (The large forces are scaled for drawing.)

\(Clear[x, y, H, w, Y0, a, b]; 
V[x_] := H\ Sinh[\(w\ x\)\/H]; 

    y[x_] := \(\(H\ \)\/w\) \((Cosh[\(w\ x\)\/H] - 1)\) + Y0; 

    T[x_] := {H, V[x]}; 
H = 700; w = 10; Y0 = 0; 
a = 10; 
b = 50; 

    curve = Plot[y[x], {x, a, b}, Axes -> None, 
        PlotStyle -> {AbsoluteThickness[4], Red}]; 
scale = 10; 

    arrows = Show[{Arrow[T[b]/scale, Tail -> {b, y[b]}, Brown], 
          Arrow[\(-T[a]\)/scale, Tail -> {a, y[a]}, Brown], 
          Arrow[{H, 0}/scale, Tail -> {b, y[b]}, Green], 
          Arrow[{\(-H\), 0}/scale, Tail -> {a, y[a]}, Green], 
		
          Arrow[{0, V[b]}/scale, Tail -> {b, y[b]}, Blue], 
          Arrow[{0, \(-V[a]\)}/scale, Tail -> {a, y[a]}, Blue], 
		
          Arrow[{0, \(-V[b]\) + V[a]}/scale, 
            Tail -> {\((a + b)\)/2, y[\((a + b)\)/2]}, SteelBlue]}]; 

    Show[curve, arrows, PlotRange -> All]; \)

Notice that the sum of the short blue and light blue arrows is a downward arrow the same size as the long upward blue arrow.

A Microscopic Diagram of Forces

The following program draws a tiny segment of a catenary and shows the tension  forces acting on the segment.  (Forces are scaled for drawing.)

\(Clear[x, y, H, V, T, w, Y0, a, b]; 
V[x_] := H\ Sinh[\(w\ x\)\/H]; 

    y[x_] := \(\(H\ \)\/w\) \((Cosh[\(w\ x\)\/H] - 1)\) + Y0; 

    T[x_] := {H, V[x]}; 
H = 700; w = 10; Y0 = 0; 
a = 20; 
b = 20.1; 

    curve = Plot[y[x], {x, a, b}, Axes -> None, 
        PlotStyle -> {AbsoluteThickness[4], Red}]; 
scale = 2000; 

    arrows = Show[{Arrow[T[b]/scale, Tail -> {b, y[b]}, Brown], 
			
          Arrow[{H, 0}/scale, Tail -> {b, y[b]}, Green], 
		
          Arrow[{0, V[b]}/scale, Tail -> {b + H/scale, y[b]}, Blue]}]; 

    Show[curve, arrows, 
	
      Graphics[{Line[{{a, y[a]}, {b, y[a]}, {b, y[b]}}], 
          Text["dx", {\((a + b)\)/2, y[a] -  .02}], 
          Text["dy", {b +  .02, \((y[a] + y[b])\)/2}], 
		
          Text["T",  .5  T[b]/scale + {b, y[b] +  .02}], 
		
          Text["H", { .5  H/scale, 0} + {b, y[b] -  .02}], 
		
          Text["V", {H/scale,  .5  V[b]/scale} + {b +  .02, y[b]}]}], 
      PlotRange -> All]; \)

Why is the aparant similarity

(dy/dx) ≈  (V/H)

important?

Solution of Differential Equations

Mathematica needs help to solve the differential equations in the general case:

Clear[x, y, V, H, w, Y0];
DSolve[{V[0] == 0,V'[x]==\(w\ \@\(1 + V[x]\^2\/H\^2\)\)},{V[x]},x]

Clear[x, y, V, H, w, Y0];
Y0 = 1;w = 1;H = 1;
DSolve[{y[0] == Y0,V[0] == 0,y'[x]==V[x]/H,V'[x]==\(w\ \@\(1 + V[x]\^2\/H\^2\)\)},{y[x], V[x]},x]

Hint 13.4

Mathematica can help check the answer:

Clear[V, x, H, w];
V[x_] = H Sinh[((w x)/H)];
(((((V'[x])/w)))^2)==\(Simplify[1 + V[x]\^2\/H\^2]\)

((V'[x])/w)

\(1 + V[x]\^2\/H\^2\)

\(Simplify[1 + V[x]\^2\/H\^2]\)

Hint 13.5

(dy/(dx )) =  Sinh[(w/H) x]

so

\(TextForm
    \`y[x]\  = \ \(H\/w\) \((Cosh[\(w\/H\) x]\  - \ 1)\)\  + \ Y0\)

because ...

\(Clear[x, y, w, H, Y0]; 

    y[x_] = \ ∫Sinh[\(w\/H\) x] \[DifferentialD]x\  + \ constant
        \)

y[0]

Solve[y[0] == Y0, constant]

Or check

\(Clear[x, y, w, H, Y0]; 

    y[x_]\  = \ \(H\/w\) \((Cosh[\(w\/H\) x]\  - \ 1)\)\  + \ Y0\)y[0]y'[x]\(w = 5; 
H = 1000; 
Y0 = 20; 
Plot[y[x], {x, \(-75\), 75}]; \)

Hint 13.7

Read the print project and use the following help.

\(Clear[x, y, w, H, Y0]; 

    y[x_]\  = \ \(H\/w\) \((Cosh[\(w\/H\) x]\  - \ 1)\)\  + \ Y0\)y[75] - y[0]

\(w = 5; 
FindRoot[y[75] == y[0] + 10, {H, 1000}]\)

Hint 13.8

Read the print project.  Compare your ideas with the computation:

\(Clear[x, y, H, w, T, V]; 
V[x_] := H\ Sinh[\(w\ x\)\/H]; 

    T[x_] := {H, V[x]}; 
Norm[T[x]]\)Simplify[%]

Where is the scalar tension  t[x] = Norm[T[x]] largest?

\(Clear[t, x, w, H]; 
t[x_] := \ H\ Cosh[\(w\ x\)\/H]\)t[x]t'[x]

Why is t[x] increasing for x > 0?  What happens to the derivative of an increasing function?

\(w = 5; 
H = 1441; 
Plot[t[x], {x, 0, 75}, PlotStyle -> Brown]; 

    Plot[\(t'\)[x], {x, 0, 75}, PlotStyle -> Orange, AspectRatio -> 1]; \)

Section 13.5 Optimal Strength

Here are some computational hints on the hints...  Read the print project section.

\(Clear[t, r, w, H]; 
tMAX[H_] := \ H\ Cosh[\(w\ r\)\/H]\)tMAX'[H]

Let  ξ = ((w r)/H) and look for a critical value to the general problem:

\(Clear[ξ]; 
Solve[ξ\ Sinh[ξ]\/Cosh[ξ] == 1, ξ]\)

Try for a numerical solution:

FindRoot[ξ Sinh[ξ]/Cosh[ξ] == 1, {ξ, 1}]

How does the maximal tension vary with  H?

\(w = 5; r = 75; 
Plot[tMAX[H], {H, 100, 5000}]; \)

FindRoot[tMAX'[H] == 0, {H, 500}]

The solution with  ξ:

H = r w/1.19968

tMAX[312.584]

tMAX[1441. ]tMAX[144. ]