222 - Topic 12 - Multivariate Extrema - Aluminum Trough

2800 days ago by Professor222

Applied Extrema for Multivariate Functions

John Travis

Mississippi College

Experimenting with the standard problem of determining the best way to crimp a sheet of aluminum so that the resulting trough can hold the largest possible amount of material.

# Folding up a sheet of aluminum to create a trough with maximum cross-sectional area L = 15 # width of the aluminum sheet var('x,theta') A = (L-2*x)*x*sin(theta) + x^2*sin(theta)*cos(theta) A(x,theta) = A @interact def _(crimp=slider(1,L/2,1/10,1),angle=slider(1,89,1,10)): global x0 x0 = crimp global a0 anglerad=angle*pi/180 a0 = anglerad h = crimp*sin(anglerad) b = crimp*cos(anglerad) G = line([(0,h),(b,0),(b+L-2*crimp,0),(2*b+L-2*crimp,h)],thickness=5) G += text('Area = %s'%str(A(crimp,anglerad).n()),(L+1,h/2),color='black') G.show(aspect_ratio=1) Surf = plot3d(A,(x,0,L/2),(theta,0,pi/2),color='yellow')+point3d((x0,a0,A(x0,a0)),size=10) Surf.show() C = contour_plot(A,(x,0,L/2),(theta,0,pi/2),contours=30,labels=True) C += point((x0,a0),zorder=2,color='red') C.show() 
       
crimp 
angle 

Click to the left again to hide and once more to show the dynamic interactive window

# Let's look for the critical point analytically fx = diff(A,x) ft = diff(A,theta) grad = [fx,ft] print 'The gradient vector' show(grad) 
       
The gradient vector
The gradient vector
 
       
# Since sin(theta)=0 does not yield a maximum for this problem, we can divide it out. # fx = fx/sin(theta) solns = solve([fx,ft],(x,theta),solution_dict=True) print(html('The collection of solutions to grad(f)=0')) # show((sol[x],sol[t]) for sol in solns) 
       
Traceback (click to the left of this block for traceback)
...
ValueError: self must be a relation
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "_sage_input_20.py", line 10, in <module>
    exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("IyAgU2luY2Ugc2luKHRoZXRhKT0wIGRvZXMgbm90IHlpZWxkIGEgbWF4aW11bSBmb3IgdGhpcyBwcm9ibGVtLCB3ZSBjYW4gZGl2aWRlIGl0IG91dC4KIyBmeCA9IGZ4L3Npbih0aGV0YSkKc29sbnMgPSBzb2x2ZShbZngsZnRdLCh4LHRoZXRhKSxzb2x1dGlvbl9kaWN0PVRydWUpCnByaW50KGh0bWwoJ1RoZSBjb2xsZWN0aW9uIG9mIHNvbHV0aW9ucyB0byBncmFkKGYpPTAnKSkKIyBzaG93KChzb2xbeF0sc29sW3RdKSBmb3Igc29sIGluIHNvbG5zKQ=="),globals())+"\\n"); execfile(os.path.abspath("___code___.py"))
  File "", line 1, in <module>
    
  File "/tmp/tmp2lSDGG/___code___.py", line 4, in <module>
    solns = solve([fx,ft],(x,theta),solution_dict=True)
  File "/home/sageserver/sage-7.5/local/lib/python2.7/site-packages/sage/symbolic/relation.py", line 875, in solve
    sol_dict=[{eq.left():eq.right()} for eq in sol_list]
  File "sage/symbolic/expression.pyx", line 2379, in sage.symbolic.expression.Expression.right_hand_side (/home/sageserver/sage-7.5/src/build/cythonized/sage/symbolic/expression.cpp:16619)
ValueError: self must be a relation
# Even after simplifying by dividing out sin(theta) in the first partial, things still don't solve well. d1 = (diff(A,x)/sin(theta)).simplify_trig() d2 = (diff(A,theta)/x).simplify_trig() show(d1) show(d2) 
       

                                
                            

                                
# After trying to evaluate this cell, you may have to use the "Interupt" option in the "Action" drop down # box at the top of this page. You will not likely succeed symbolically by trying to solve these equations directly. solve([d1,d2],x,theta) 
       
[2*x*(cos(theta) - 2) + 15, (2*cos(theta)^2 - 2*cos(theta) - 1)*x +
15*cos(theta)]
[2*x*(cos(theta) - 2) + 15, (2*cos(theta)^2 - 2*cos(theta) - 1)*x + 15*cos(theta)]
# However, by doing a bit of substitution and eliminating the cos(theta) terms we get the following # expression in terms of x only. That one can be solved for x and then plugged back in to get theta. eqn = (L*x-2*x^2)*(4*x-L)/(2*x)+x^2*(2*(4*x-L)^2/(2*x)^2-1) eqn.show() solns = solve(eqn,x,solution_dict=true) for soln in solns: x0 = soln[x] if x0<>0: t0 = acos((4*soln[x]-L)/(2*soln[x])) html('x = '+str(soln[x])+' and $\\theta$ = '+str(t0)) 
       
x = 5 and  = 1/3*pi
                                
                            
x = 5 and  = 1/3*pi